电力系统状态估计这个话题，老早以前是SCADA的天下，调度员靠RTU传来的有功、无功和幅值，再用非线性加权最小二乘去迭代，一套下来动不动几十次迭代，碰上坏数据还得来回排查。这几年PMU（相量测量单元）铺开之后，局面变了不少——它直接给出带时标的电压相量，幅值有了，相角也有了，状态估计从非线性一下子被拉回到线性模型，算起来清爽太多。今天想聊的这个项目，就是用MATLAB搭一套基于PMU量测的电力系统状态估计，标题叫《基于matlab PMU相量测量单元电力系统状态估计》，带源码，编号14925期。我按自己做同类项目的经验，把整体设计思路、数学模型、MATLAB实现细节和踩过的坑全部拆开讲一遍，给正在做电力系统课设、毕设或者入门广域监测系统的朋友一个可以直接落地的参考。
这套东西能解决什么问题？简单说：在系统里装了一批PMU之后，如何在已知网络拓扑、线路参数和部分量测数据的情况下，把全网各母线电压的幅值和相角给估出来。有了这个东西，后续的潮流追踪、扰动定位、静态稳定分析才有基础数据。适合谁看？电气工程相关专业的学生、刚接触WAMS（广域测量系统）的工程师，以及想把MATLAB数值计算能力和电力系统分析结合起来练手的开发者。下面按项目从设计到代码再到排错的顺序，从头捋一遍。
1. 项目整体设计与思路拆解
1.1 核心需求：为什么有了PMU，状态估计变简单了
传统状态估计用的是SCADA量测，量测类型主要是节点注入有功、无功、支路潮流和电压幅值。注意，这里面没有相角量测，因为SCADA根本没有统一时标去测相角。于是状态量和量测量之间是非线性关系，比如支路有功潮流公式是P_ij = V_i²G_ij - V_iV_j(G_ij cosθ_ij + B_ij sinθ_ij)，求解必须靠高斯-牛顿迭代，每次迭代都要重新算雅可比矩阵，计算量大，而且在初值离真值远的时候还可能发散。
PMU把这个痛点直接戳掉了。它能以GPS/北斗授时同步采样，输出的就是带绝对时标的电压相量，幅值V_i和相角θ_i是直接测出来的。因此状态量和量测量之间变成了线性关系：z = Hx + e。这种情况下，状态估计本质上就是一个线性加权最小二乘问题，不用迭代，一次矩阵运算就出结果，稳定性和速度都有质的提升。
打个比方：SCADA就像你只知道一个人的身高，却要推断他的站姿；PMU则是直接给你一张带角度的照片，身高、倾斜角度全都标好了。不是不需要估计，但估计的工作量小了很多。
1.2 方案选型：为什么用线性WLS而不是卡尔曼滤波
这个项目走的是加权最小二乘（WLS）路线，而且是线性版本。选它有几个理由：
第一，静态状态估计是基础。WLS是经典框架，理论成熟、代码易懂、结果可以直接跟真值做误差对比。对课设和毕设来说，这个路线做出来有理有据，答辩也好讲。
第二，线性WLS不需要迭代矩阵的反复分解。传统WLS每轮迭代都要做一次LU分解，而线性模型只需要一次求解，对MATLAB这种解释型语言来说性能友好很多。
第三，卡尔曼滤波虽然能处理动态过程，但它需要模型噪声和量测噪声的先验协方差，调参难度大。这里先做静态版本，后续要扩展成动态估计，可以在此框架上加状态转移方程，那属于进阶玩法。
结论：以线性WLS为骨架，把PMU量测方程、权重矩阵、可观测性分析、坏数据检测这几个模块串起来，是性价比最高的方案。
1.3 数据流和模块划分
整个项目的核心流程可以切分成五个模块：
- 数据输入：系统节点数、支路参数、PMU安装位置及量测值。
- 可观测性分析：检查量测是否足够让H矩阵列满秩。
- 核心估计：构建z、H、R，求解WLS状态估计。
- 坏数据检测：计算残差，做检测与剔除。
- 结果输出：展示估计值、误差指标，可视化对比。
这个划分也直接映射到MATLAB的代码结构上，每个功能对应一个函数，main脚本负责串联。后面第三章会给出各模块的代码骨架和关键参数选取逻辑。
2. 核心模型与算法拆解
2.1 PMU量测方程与状态向量的数学形式
先定义状态向量。对于节点数为n的电力系统，若所有节点都装了PMU，每个节点的状态是电压幅值V_i和相角θ_i，那么理论上的状态量总数是2n。但因为全网相角需要一个参考基准，通常把某个参考节点的相角锚定为0，实际上待估计的状态数是2n-1。
PMU的量测方程可以写成：
z = Hx + e
其中x是状态向量，e是量测噪声向量。假设在母线i上装了一台PMU，它能同时量测该母线的电压幅值和相角，还能量测与该母线相连的支路电流相量。以电压量测为例，对应H矩阵的行是第i个幅值状态位为1，第i个相角状态位为0；对应相角量测的行是第i个相角状态位为1，幅值状态位为0。如果PMU还提供支路电流相量，那需要基于线路的π型等值电路，把电流相量表达式转换为状态量的线性组合，这部分关键是要精确填H矩阵的系数，稍有差错，估计结果就会偏掉。
这里面有个细节很多人第一次做容易漏：PMU给的是绝对相角（以UTC为基准的全球同步相角），而状态估计里的状态量是相对参考母线的相角。处理方式是把所有相角量测减去参考母线的相角量测，再进入估计。否则H矩阵会多一列全零或者出现秩亏。
2.2 加权矩阵R的选取：为什么要按量测类型分权重
WLS的目标函数是J = (z-Hx)ᵀR⁻¹(z-Hx)，R是量测误差协方差矩阵。PMU的相量测量单元精度比传统RTU高很多，幅值误差通常在0.1%量级，相角误差在0.01°~0.02°量级。但不同PMU通道、不同幅值和相角的方差差异还是存在，所以R不能简单设成单位阵。
实践中常用做法是查PMU的精度指标，把幅值标准差和相角标准差换算成方差：
- 电压幅值量测：标准差约0.001~0.002 p.u.，对应R对角线元素约为1e-6~4e-6。
- 电压相角量测：标准差约0.0002~0.0004 rad，对应R对角线元素约4e-8~1.6e-7。
- 电流幅值和相角量测：取决于CT/PT的精度等级，通常比重会比电压量测的方差大一些。
权重的本质是量测的信任度：方差越小，权重越大，在求解时对结果的贡献越大。这个道理和加权平均是一样的，实际项目里我习惯先用均匀权重跑一遍，看残差分布，再用残差方差反推R，做一次迭代定权。这个小技巧能明显改善结果。
2.3 可观测性分析：为什么H矩阵必须满秩
线性状态估计能解出唯一解的前提，是量测方程个数不小于状态数，而且H矩阵列满秩。列满秩意味着每个状态量都被足够的独立量测覆盖，不存在某个母线电压怎么测都测不到的情况。
在MATLAB里判断很简单：计算秩rank(H)，列满秩的判据是rank(H)等于状态数2n-1。同时还可以看条件数cond(H)，条件数太大说明H矩阵近似病态，哪怕满秩，数值上也可能解出漫天乱跳的结果。条件数控制在1e6以内比较好，超过这个量级就要警惕。
这就引出一个部署问题：PMU数量不够怎么办？工程上常见的是PMU只装在部分关键节点，剩下的节点靠SCADA量测补齐。这种混合量测场景下，模型重新变成非线性，得用传统WLS迭代。如果非要保持线性模型，可以假设SCADA区域的状态初始值已知，那其实就不叫状态估计了，属于扰动分析，逻辑上要分清。
2.4 坏数据检测：标准化残差怎么用
PMU数据也不是百分百干净，通信丢包、相量计算异常、GPS失步都会产生坏数据。线性模型下，残差r = z - Hx_est，理论上服从零均值高斯分布。采用的是基于标准化残差的检测：
r_i_normalized = r_i / sqrt(R_ii * (I - H(HᵀR⁻¹H)⁻¹HᵀR⁻¹)_ii)
分子是第i个量测的残差，分母是残差方差的开方。标准化之后，r_i_normalized近似服从标准正态分布，用阈值λ（一般取3.0，对应99.7%置信度）去卡：超过阈值就判坏数据。
有一个容易踩的坑：多个坏数据同时存在时，逐次剔除比一次性剔除更稳。因为坏数据可能互相掩盖，残差会被拉平，单次残差检验会漏掉。每次只剔除标准化残差最大的那个量测，重新做一遍估计，再检，直到没有超阈值的点为止。后面第四章会展开讲这个问题的具体表现。
3. MATLAB实操实现与核心环节
3.1 数据准备：我用IEEE 9节点系统作为测试床
我复现这个项目时用的算例是IEEE 9节点系统，节点数少、拓扑清楚、Matpower里有现成数据，适合验证算法。数据准备阶段要明确几样东西：
- 节点表：9个节点，包括基准电压、类型（PQ/PV/平衡）。
- 支路表：每条支路的电阻、电抗、对地电纳，以及变压器变比。
- PMU位置：试验时我假定节点1、3、6、9装了PMU，量测覆盖这四点的电压相量，以及相连支路的电流相量。
在MATLAB里我推荐用struct组织这些数据，别用一堆散变量：
data.n = 9; data.branch = [ 1 4 0.0000 0.0576 0.0000 0; 4 5 0.0170 0.0920 0.1580 0; 5 6 0.0390 0.1700 0.3580 0; ... ]; data.pmu = [1; 3; 6; 9]; data.z_meas = YOUR_MEAS_VECTOR;如果手上没有实测PMU数据，可以用潮流计算结果作为真值，再叠加高斯噪声生成量测值。这个做法对验证代码正确性特别有用——因为你知道真值，就能算误差。
3.2 核心求解函数：从测量向量到状态量
这一段是代码的枢纽。给出一个核心函数框架，读者可以直接改成自己的数据规模。
function x_est = pmu_wls_se(z, H, R) % z: 量测向量 m x 1 % H: 量测矩阵 m x (2n-1) % R: 量测误差协方差矩阵 m x m G = H' * (R \ H); % 信息矩阵 b = H' * (R \ z); % 右端项 x_est = G \ b; % 最小二乘解 end实际项目里H不是手工填的，而是根据PMU位置和网络拓扑动态生成。构建H矩阵的逻辑分三步：
第一步，建立状态索引。每个节点分配两个索引：幅值索引idxV_i = 2*(i-1)+1，相角索引idxTheta_i = 2*(i-1)+2。参考节点的相角索引要特殊处理，要么在H中删掉该列，要么在x中固定为0，并同步调整量测方程。
第二步，填电压量测行。对于母线i的PMU，幅值量测行在idxV_i位置填1，相角量测行在idxTheta_i位置填1。相角量测如果用的是全局相角，记得统一减去参考节点的全局相角后再进估计。
第三步，填支路电流量测行。以π型等值电路为准，先算线路导纳Y_ij = G + jB，再根据电流相量I_ij = Y_ij(V_i - V_j) + jBsh/2 * V_i，把实部和虚部对状态量的偏导算出来，填到对应位置。这一步最容易写错，强烈建议先用一个简单两节点系统验证H矩阵的正确性。
3.3 结果验证：误差分析和残差分析
状态估计算完不能直接交差，得验证。我的做法是把估计值跟潮流真值做对比，计算每个节点的幅值误差和相角误差。通常用RMSE（均方根误差）来评价整体精度，公式是：
RMSE = sqrt(mean((x_est - x_true).^2))
从我的测试结果看，在R矩阵设置合理、量测噪声标准差符合PMU实际水平的前提下，9节点系统的幅值估计误差大约在1e-4 p.u.量级，相角估计误差大约在1e-3 rad量级。如果误差偏大一个数量级以上，优先检查H矩阵有没有填错，其次是R矩阵是否给了不合理的权重。
另外要画残差分布图。把标准化残差画成条形图，直观能看出有没有异常量测。正常情况下残差密布在±3之间，且没有明显单点突出；如果有突出点，先别急着删，确认是数据问题还是H矩阵问题。
3.4 整体脚本结构与运行流程
main脚本的结构可以这样安排：
% 第1步：加载系统数据和PMU量测 system_data = load_pmu_system('ieee9'); % 第2步：构建H矩阵和R矩阵 [H, R, z, idx] = build_linear_model(system_data); % 第3步：可观测性检查 assert(rank(H) == size(H,2), 'H矩阵秩亏，系统不可观测'); % 第4步：状态估计 x_est = pmu_wls_se(z, H, R); % 第5步：坏数据检测 [r_norm, bad_idx] = bad_data_detect(z, H, R, x_est); % 第6步：结果输出 plot_result(x_est, system_data);这里每步调用一个函数，函数内部再细分，可读性和可维护性都比写一个两三百行的主脚本强得多。有一个小建议：运行前用tic/toc记录时间，线性模型下9节点系统的计算时间应该在毫秒级，G矩阵是17x17，MATLAB分解起来非常快。如果跑出来要好几秒，基本可以确定有冗余循环或者H矩阵构建逻辑低效，需要排查。
4. 常见问题与排查技巧实录
4.1 H矩阵奇异或条件数过大怎么办
这是我在项目里遇到最多的一个问题，几乎每个初做PMU状态估计的人都会撞上一次。表现就是rank(H) < size(H,2)，或者cond(H)在1e12以上，估计值乱跳。
可能原因有三个。一是PMU覆盖不足，某些母线完全没有量测覆盖，对应H行全是零，状态不可观。二是参考相角没有处理干净，H矩阵里含有一列全零或者两列线性相关。三是线路参数填错，导致支路电流量测对应行与电压量测行产生线性相关关系。
排查顺序：第一步，打印H矩阵的稀疏模式，用spy(H)看哪列全零，哪两列是成比例关系。第二步，逐个PMU检查量测方程个数和类型，确认覆盖范围。第三步，用一个只有两台PMU的两节点系统做单元测试，H矩阵规模小，一眼能看出问题。
提示：如果是可观测性不足，不要试图在代码层面打补丁。老老实实增加PMU量测，或者把部分SCADA量测补进模型，用混合量测的非线性WLS。强行求解只会得到数值上看起来正常、实际上完全没意义的解。
4.2 相角参考基准冲突
PMU量测给出的是全球同步的绝对相角，不同厂家的PMU在接入同一系统时，由于GPS信号处理延迟的差异，可能会有微小的角度偏移。如果直接把不同PMU的相角拿来拼成同一个z向量，容易出现系统性偏差。
我的处理办法是：在数据预处理阶段，先把所有PMU相角量测减去同一个参考PMU的相角量测，得到相对相角序列再进入估计。这样即使PMU本身有固定延迟误差，只要延迟在短时间内稳定，相减后误差会被抵消掉。这个操作在代码里就是一行：z_theta = z_theta - z_theta_ref。
4.3 坏数据检测的误检与漏检
误检通常是因为R矩阵的方差设得太小，量测噪声本来没那么高精度，标准化之后残差就会偏大，超过阈值。漏检则常见于两个坏数据互相抵消的场景。举个例子：某条支路两端的电流量测同时坏掉，它们的残差可能方向相反，平均下来标准化残差不大，就会被漏掉。
实操中我的建议是：别只依赖一次残差检验。做一个循环剔除法——每次只删标准化残差最大的那个量测，重新估计后再检。同时把检验阈值从3.0放宽到2.8，多捕获一些边缘可疑点，宁可多剔除一个可疑量测，也不能放坏数据进门。当然，剔除的量测数不能太多，一般超过总量测数的5%就要回头检查是不是数据质量整体不行。
4.4 量测噪声设置与实际不符
测试阶段想要模拟真实场景，可以在潮流真值上加高斯白噪声，但噪声标准差的选择要有依据。我见过有人直接用randn加噪声，标准差设成0.1，结果估计误差大得离谱，还以为是算法有问题，实际上是噪声水平跟PMU的真实指标差了百倍。
PMU的幅值测量精度典型值在0.1%，相角测量精度在0.01°左右，对应弧度约1.7e-4 rad。按这个量级设置噪声，估计结果才能反映算法本身的性能。做完之后可以统计残差的标准差，和设置的噪声标准差做对比，如果差太多，说明H矩阵或者加权有问题。
4.5 计算效率与内存小技巧
PMU量测点数一旦多了，H矩阵和R矩阵的维度会涨得很快。比如IEEE 118节点系统全装PMU，状态量就是235个，量测可能有上千行。这时候如果H还是稠密矩阵，求逆操作会越来越慢。
两个改进方向：一是用稀疏矩阵存储H和R，MATLAB里直接sparse(H)；二是用信息矩阵G = H'RH，然后对G做Cholesky分解，而不是直接求H的伪逆。这两步能把计算时间下降两个数量级。我在118节点系统上试过，从几十秒降到几百毫秒，效果非常明显。
另外一个细节：R矩阵是对角阵，R \ H这一步不要写成inv(R) * H，直接用左除效率更高，数值也稳定得多。MATLAB里对稀疏对角阵的除法有专门优化，一定要利用上。
最后再分享一个扩展思路。这个静态估计框架跑通之后，如果还想往深做，可以尝试两个方向：一是把PMU量测数据按时间序列连续输入，加入状态转移模型，升级成动态状态估计（这时就该上卡尔曼滤波了）；二是在现有框架中把H矩阵的构建推广到三相不平衡系统，就能用来处理配电网状态估计。我个人做下来最大的体会是，这套模型的代码骨架一旦搭好，往各种方向扩展都很快，关键是前期的数据结构和H矩阵构建逻辑要写干净，别为省几行代码把后续的扩展性毁了。
