做电力系统状态估计的同学，一定绕不开WLS、PMU和Newton-Raphson这三样东西。我前阵子刚把一个完整的对比实验跑通：用PMU量测加上加权最小二乘（WLS）估计出系统的电压幅值和相角，再拿这个估计结果和Newton-Raphson潮流算出来的状态做对比，全程用Matlab实现。整个过程比想象中曲折，但跑完之后，很多以前只停留在课本上的概念，比如量测冗余度、雅可比矩阵、可观测性，都变得非常具体。这篇博客就把这个项目从原理到代码，再到调试踩坑的过程完整记录下来，给准备入坑电力系统状态估计的同行一个可以直接参考的完整样例。
这个项目非常适合两类人：一类是刚接触状态估计、想把理论公式变成可运行代码的研究生，另一类是已经在做EMS/SCADA相关工程，想评估PMU数据接入后对状态估计精度影响的工程师。你不需要很深的Matlab功底，但最好懂一点误差分析和线性代数，否则迭代公式那块会有点懵。我会把每条公式都拆开讲，代码也会给关键片段，你照着抄基本能复现。
1. 项目背景与核心思路拆解
1.1 这个项目到底在做什么
先明确一下任务边界。常规电力系统状态估计，会用SCADA系统提供的量测数据（节点注入功率、支路潮流、电压幅值）作为输入，通过某种估计算法求出系统最可能的电压状态，也就是每个节点的电压幅值V和相角θ。而PMU（相量测量单元）是近年来大量部署的同步测量装置，它借助GPS授时能够直接测得带时标的电压相量，也就是V和θ一次同时给出，精度还很高。
这个项目做的事情可以拆成四步：
- 生成一个IEEE标准测试系统的数据，我这边用的是IEEE 14节点系统。
- 先跑一次Newton-Raphson潮流，得到精确的电压幅值和相角，作为“真值”。
- 在真值上叠加一定的测量噪声，模拟PMU实际输出的带误差量测数据。
- 用WLS算法基于这些带噪声的量测做状态估计，再把估计结果和第一步的潮流真值放在一起比较，画出误差曲线，计算出平均偏差和均方根误差。
逻辑上这是一条很清晰的验证链：NR潮流提供参考基准，PMU量测提供观测数据，WLS负责从观测中还原状态。整个过程不需要真实的PMU硬件，也不需要联网实时数据，只需要一台电脑和Matlab，非常适合做算法验证和教学演示。
1.2 为什么选WLS和PMU这对组合
状态估计算法有很多种，但WLS（Weighted Least Squares，加权最小二乘）一直是工业界和学术界使用最广泛的经典算法。原因很简单：在量测误差服从正态分布的前提下，WLS估计是线性无偏且方差最小的估计器之一，数学性质非常漂亮。它的目标函数是所有量测残差乘以权重后的平方和，形式直观，求解过程也比较稳定。
PMU和WLS搭配的原因更实用。传统SCADA量测是非同步的，而且功率量测和电压幅值量测属于不同类型，量测方程是非线性的。PMU则直接提供同步的电压和电流相量，量测方程在有公共参考相角的前提下是线性的，整个估计问题可以简化，甚至在某些配置下可以变成线性WLS，一步算出结果。即便仍然使用非线性迭代，PMU在直角坐标系下的量测方程也要比极坐标下的功率方程简单不少，雅可比矩阵的计算量会下降。
另外，项目的目标不是单纯做状态估计，而是要和Newton-Raphson潮流做对比。NR潮流本质上是求解一个确定性的非线性方程组，输入是负荷和发电机出力，输出是唯一的运行状态。而状态估计输入是量测数据，输出是一个统计意义上的“最优估计”。两者在数学上同宗同源，都离不开雅可比矩阵和迭代求解，但应用场景完全不同。把这两种方法放在同一个框架里比较，能非常直观地理解“精确求解”和“噪声中估计”的差别。
1.3 为什么要把估计结果和Newton-Raphson潮流比较
我最早看到这个课题的时候也有疑问：状态估计的结果应该和谁比？如果没有真值，那只能看残差；但残差小不代表估计准确，因为可能有坏数据污染。为了客观评价估计算法，我们需要一个参考标准。
在这个项目里，Newton-Raphson潮流就是那个标准。给定相同的网络参数和运行点，潮流计算解出的电压幅值和相角，就是在无噪声的确定性条件下系统唯一的状态。然后把潮流解叠加噪声当量测，再让WLS去估计，估计值和潮流解之间的差距就是估计算法在噪声环境下的真实误差。这种做法的好处是可控性强：你可以在不同噪声水平、不同量测配置下反复实验，量化WLS和PMU组合对状态估计精度的提升幅度。
从工程角度看，这种对比也有实际意义。在线状态估计系统往往会把潮流计算结果作为粗校验，比如比较调度员潮流和状态估计结果的偏差，偏差过大说明量测数据有问题。所以搞明白WLS估计结果和NR潮流结果之间的误差分布特征，对于实际系统的数据质量分析很有帮助。
2. 电力系统状态估计的数学模型
2.1 状态向量与量测方程
电力系统状态估计的状态向量通常取所有节点的电压幅值V和相角θ，也就是：
x = [θ_2, θ_3, ..., θ_n, V_1, V_2, ..., V_n]^T
注意一般把参考节点（平衡节点）的相角固定为0，因此状态向量中不包含θ_1。这是因为相角是相对量，必须有一个参考基准。实际网络中的功率量测只和相角差有关，PMU电压相量量测则必须经过坐标变换才能放进这个以节点1为参考的状态框架里，这也是后面代码里要特别处理的地方。
量测方程可以统一表示为：
z = h(x) + e
其中z是量测向量，h(x)是量测函数，e是量测误差向量，通常假设为零均值高斯白噪声，其协方差矩阵为R。量测向量里可以包含：
- 支路有功和无功潮流：P_ij、Q_ij；
- 节点注入有功和无功：P_i、Q_i；
- 节点电压幅值：V_i；
- 节点电压相角：θ_i；
- 支路电流幅值和相角：I_ij、φ_ij。
前两类是传统SCADA量测，后三类是PMU可以提供的量测。在这个项目里，为了让问题聚焦在PMU的贡献上，量测向量主要使用PMU量测，也就是节点电压幅值、节点电压相角和支路电流相量。其中支路电流相量又可以通过线路导纳转化为节点电压相量的函数，因此整体量测方程依然只和电压状态有关。
2.2 WLS估计的迭代公式与雅可比矩阵
WLS的目标函数是加权残差平方和最小：
J(x) = [z - h(x)]^T R^{-1} [z - h(x)]
其中R是量测误差协方差矩阵。如果量测互相独立，R就是一个对角阵，对角线元素是各量测的方差σ_i^2。权重越大，表示该量测越可信。WLS的求解思路和牛顿法一样，对目标函数求梯度并令其为零，然后通过高斯-牛顿法迭代：
Δx = (H^T R^{-1} H)^{-1} H^T R^{-1} [z - h(x)]
这里H是h(x)对状态x的雅可比矩阵，也就是量测函数对每个状态变量的偏导数组成的矩阵：
H = ∂h(x)/∂x
每次迭代更新状态：
x^(k+1) = x^(k) + Δx
直到Δx的范数小于某个阈值，比如1e-6。
整个计算的核心就在雅可比矩阵H的构造上。如果量测是PMU电压相量，那么在极坐标模型下，V_i量测对V_i的偏导为1，对θ_i的偏导为0；θ_i量测对θ_i的偏导为1，对其他相角偏导为0，形式非常简单。支路电流相量稍复杂，由于电流I_ij和电压V_i、V_j都有关系，需要用线路导纳的实部和虚部把电流相量分解成电压相量的线性函数，再求偏导。若使用直角坐标模型，整个量测方程会进一步线性化，这里为了和NR潮流统一，还是采用极坐标形式。
2.3 PMU量测方程的特殊处理
PMU量测与传统功率量测最大的不同，是它可以直接测量相角。但相角必须有一个参考。PMU自身的参考是GPS时间信号，所有PMU测得的相角都是相对于全球参考时间基准的绝对相角。但在状态估计里，我们把节点1的相角固定为0，所以PMU电压相角量测不能直接拿来用，必须先做参考坐标变换。
假设PMU测得的节点i电压相角为θ_i^PMU（相对于GPS），节点1的参考相角为θ_1^PMU（也由PMU给出，但通常认为等于0），那么在状态估计模型里实际的量测方程应该写成：
θ_i^cal = θ_i - θ_1
其中θ_i和θ_1是状态向量里待求的相角。如果θ_1没有被当作状态量，则量测方程就是θ_i^cal = θ_i，但在实际代码中，参考节点的具体值可能需要特殊处理。简单的做法是：把PMU电压相角量测统一减去PMU量测到的参考节点相角，变成相对参考节点的相角差，再进入状态估计。我在代码里就是这么处理的，否则信息矩阵很容易奇异，迭代直接失败。
另一个容易忽略的问题是支路电流相量量测。电流相量的相角也是相对于GPS参考的，同样需要去参考。而且电流相量的大小与线路两端电压相关，需要对每条支路单独建立量测方程。在实际工程中，为了简化，很多人只使用PMU的电压相量量测，而把电流相量量测作为附加冗余信息。这个项目里我会同时考虑电压和电流两类PMU量测，但电流量测在程序里可以通过开关控制是否启用，方便做对比实验。
3. Newton-Raphson潮流计算与对比逻辑
3.1 潮流计算方程与牛顿法原理
Newton-Raphson方法求解潮流时，把节点功率平衡方程作为核心约束。对每个PQ节点，需要满足注入有功和无功等于该节点所有支路潮流之和；对PV节点，只需要满足有功平衡，电压幅值已知，无功作为待定量由计算得到。对应的方程组可以写成：
ΔP_i = P_i^spec - P_i^cal(x) = 0 ΔQ_i = Q_i^spec - Q_i^cal(x) = 0
其中P_i^spec和Q_i^spec是给定的节点注入功率，P_i^cal和Q_i^cal是由当前电压状态计算出的注入功率。然后同样构造雅可比矩阵J，求解修正方程：
J Δx = -ΔP / ΔQ
迭代更新状态。本质上，潮流计算和WLS状态估计的迭代框架非常像，都是用牛顿法求解一组非线性方程，区别只在于潮流方程的个数等于未知状态量个数，方程组是“方阵”，不存在冗余；而状态估计的量测方程个数远大于未知状态个数，是一个超定方程组，所以要用最小二乘而不是直接求解。
3.2 潮流与状态估计的异同
我把两者的区别整理成一张表，方便对照：
| 对比项 | Newton-Raphson潮流 | WLS状态估计 |
|---|---|---|
| 输入数据 | 负荷和发电机出力 | 冗余的量测数据 |
| 方程数量 | 与未知数相同，方阵 | 大于未知数，超定 |
| 求解目标 | 精确满足功率平衡 | 最小化加权残差平方 |
| 量测噪声 | 不考虑 | 统计建模，考虑噪声 |
| 核心矩阵 | 潮流雅可比J | 信息矩阵H^T R^{-1} H |
| 输出状态 | 确定性的电压幅值/相角 | 最优估计的电压幅值/相角 |
从这个表能直观理解为什么不能用潮流直接代替状态估计：实际量测有误差，如果直接拿量测代入潮流方程，方程一般不成立；状态估计通过大量冗余量测，能把误差的影响摊薄，得到更可靠的运行状态。反过来，NR潮流适合做离线分析和规划，因为输入数据是确定的。
3.3 如何构造对比场景
在这个项目里，我采用的方法是先计算潮流，然后把潮流解视为“理想真值”。接着在真值上叠加高斯噪声，生成模拟PMU量测。噪声的标准差取PMU的典型精度，电压幅值量测的标准差一般取0.001~0.005 p.u.，电压相角量测的标准差取0.001~0.005 rad，电流相量的精度也类似。为了验证算法的统计特性，可以多次随机生成量测，重复估计，统计估计误差的均值和方差。
比较时有三组关键指标：
- 最大绝对误差：max |V_est - V_true|，max |θ_est - θ_true|；
- 平均绝对误差：mean |V_est - V_true|；
- 均方根误差：sqrt(mean((V_est - V_true)^2))。
通过这三组指标，可以量化不同噪声水平、不同PMU配置下WLS估计的精度。还有一个额外的比较维度是迭代次数和计算耗时，虽然这个项目规模小，看不出太大差别，但能训练你评估算法效率的思维。
4. Matlab代码实现全流程
4.1 算例选择与量测数据生成
我选用IEEE 14节点系统，原因是规模适中，节点数不算多，但足够体现状态估计的复杂性，又不至于让代码调试变成灾难。Matlab里可以直接用Matpower工具箱加载ieee14.m数据，也可以自己手写14节点的线路、变压器和负荷数据。为了减少依赖，我在代码里用了Matpower读取数据，但核心WLS和NR算法都是自己实现的，没有调用现成函数。
量测数据生成步骤如下：
- 用自写的NR潮流函数计算系统精确电压状态；
- 选定要配置PMU的节点和线路；
- 对每个PMU量测，在精确电压和电流相量上叠加高斯噪声；
- 把带噪声的量测构造成z向量，同时生成对应的误差协方差矩阵R；
- 保存真值和量测数据，供估计使用。
这里有一个经验：如果你在Matlab里用randn生成噪声，记得固定随机数种子，比如rng(2024);，否则每次实验结果都不一样，不方便排查问题。
4.2 WLS-PMU状态估计核心代码
下面这段是WLS状态估计的主迭代代码，我简化了雅可比矩阵的生成细节，但保留了完整流程。雅可比矩阵部分需要根据量测类型逐行填充，这里略过具体的偏导推导，重点展示牛顿迭代框架。
function [V_est, theta_est, iter] = wls_pmu_estimator(z, R_inv, Ybus, ref_idx) % WLS状态估计主函数 % 输入: z - 量测向量, R_inv - 权重矩阵(逆协方差), Ybus - 导纳矩阵 % ref_idx - 参考节点索引 % 输出: V_est, theta_est - 估计的电压幅值和相角 n = size(Ybus, 1); % 节点数 V = ones(n, 1); % 电压幅值初值 theta = zeros(n, 1); % 相角初值 theta(ref_idx) = 0; % 参考节点相角固定 max_iter = 20; tol = 1e-6; for iter = 1:max_iter % 计算量测函数 h(x) hx = compute_hx(V, theta, Ybus); % 计算残差 residual = z - hx; % 构建雅可比矩阵 H H = build_jacobian(V, theta, Ybus); % 求解信息方程 G = H' * R_inv * H; dx = (G) \ (H' * R_inv * residual); % 更新状态 theta = theta + dx(1:n); V = V + dx(n+1:2*n); % 修正参考节点相角 theta = theta - theta(ref_idx); % 判断收敛 if norm(dx, inf) < tol break; end end V_est = V; theta_est = theta; end注意信息矩阵G在量测配置不足时可能奇异。比如如果系统里完全没有功率量测，只有少数几个节点的电压幅值量测，那么某些不可观测量对应的信息矩阵项可能为零，\运算会报错。我的建议是先用rank(G)检查一下矩阵是否满秩，或者对参考节点对应的行和列做特殊处理。
4.3 Newton-Raphson潮流核心代码
NR潮流代码比较成熟，我直接写了一个标准实现，重点在雅可比矩阵稀疏填充和节点类型处理。下面是核心迭代段：
function [V, theta, converged] = nr_powerflow(Ybus, Sbus, V0, theta0, ref_idx) % Newton-Raphson潮流求解 % Sbus - 节点注入功率向量, V0/theta0 - 初值 V = V0; theta = theta0; max_iter = 30; tol = 1e-8; for iter = 1:max_iter % 由当前状态计算节点注入功率 S_cal = V .* conj(Ybus * (V .* exp(1j * theta))); % 计算失配量 dP = real(Sbus - S_cal); dQ = imag(Sbus - S_cal); % 去掉参考节点的失配量 dP(ref_idx) = 0; dQ(ref_idx) = 0; % 构建潮流雅可比 J = build_powerflow_jacobian(V, theta, Ybus); % 求解修正量 dx = -J \ [dP; dQ]; % 更新 theta = theta + dx(1:n); V = V + dx(n+1:2*n); if norm([dP; dQ], inf) < tol converged = true; return; end end converged = false; end这段代码的雅可比矩阵build_powerflow_jacobian是潮流计算的关键，公式书中都有，不在这里展开。需要提醒的是，NR潮流对初值比较敏感，对IEEE 14节点系统用平启动（V=1, θ=0）基本都能收敛；如果一个算例不收敛，先检查节点类型和参数是否有误，不要急着调初值。
4.4 结果输出与误差统计
估计完成后，我习惯用一张图把结果可视化，包含四个子图：
- 节点电压幅值的真值、估计值、NR潮流值对比；
- 节点相角的真值、估计值、NR潮流值对比；
- 电压幅值估计误差的柱状图；
- 相角估计误差的柱状图。
同时打印误差统计表格，格式类似下面这样：
| 误差项 | 电压幅值 (p.u.) | 相角 (deg) |
|---|---|---|
| 最大绝对误差 | 0.0023 | 0.031 |
| 平均绝对误差 | 0.0008 | 0.012 |
| 均方根误差 | 0.0011 | 0.016 |
这里的数值取决于你设置的噪声标准差。如果PMU精度较高，误差会明显小于SCADA量测下的估计误差，这也是PMU项目最直接的收益体现。
5. 仿真结果与影响参数分析
5.1 典型算例的收敛性与精度
在默认参数下，WLS估计一般迭代3到5次就可以收敛，NR潮流也差不多。精度的核心影响因素是量测噪声标准差。我做了几组实验：当电压幅值量测噪声标准差设为0.002 p.u.、相角噪声标准差设为0.002 rad时，估计值相对真值的最大幅值误差在0.001~0.003 p.u.之间，相角误差在0.01~0.03度之间。这个精度符合PMU作为高精度量测的预期。
如果你发现估计误差比设置的噪声标准差还大很多，那一定不是噪声问题，而是代码有bug。最常见的问题有两个：一是雅可比矩阵求错了，尤其是支路电流量测的偏导；二是参考节点相角没有处理好，导致全部相角估计值整体偏移。我的调试经验是：先跑一个纯电压量测的简单场景，把相角量测全部去掉，只估计幅值，看幅值误差是否正常；再加相角量测，分步定位问题。
5.2 PMU量测误差和权重对结果的影响
权重矩阵R_inv的设定直接决定估计结果。如果某个量测的权重设得过大，算法会不分青红皂白地“信”它，一旦该量测有较大偏差，整体估计就会被带偏；如果权重设得过小，该量测就几乎没有贡献，相当于没有使用这份数据。
我做过一个特别直观的实验：把某个PMU电压相量量测的权重设成其他量测的10000倍，结果估计出的整个系统相角都被拉到接近这个量测的值，哪怕附近节点的电压相量量测明显不一致。这说明权重不是越大越好，要和量测的真实精度匹配。所以我在代码里不是直接填权重，而是填R_inv，也就是误差协方差矩阵的逆。每个量测的方差需要根据PMU技术手册和实际标定结果来设定，不能一拍脑袋。
5.3 可观测性分析与PMU布点策略
状态估计的前提是系统可观测，也就是量测数据足以唯一确定所有状态。PMU配置如果太少，系统就可能不可观测，信息矩阵奇异。在这个项目里，我可以做一个小实验：只在一个节点配置PMU电压相量量测，其他量测全砍掉，结果WLS直接报错。
工程上解决可观测性的办法是配置足够的PMU，让每个节点至少被一个PMU覆盖。PMU布点优化是一个热门研究课题，常见思路是：在保证可观测的前提下，用最少的PMU覆盖所有节点。代码里可以通过遍历所有组合来验证，但节点数大了以后是典型的组合爆炸问题，需要用整数规划或者启发式算法。这个项目里我只做了简单的覆盖率检查，也就是枚举所有节点，看是否存在至少一个与之关联的PMU量测。
6. 实操中的常见问题与解决办法
6.1 信息矩阵奇异
如果H^T * R_inv * H的秩小于状态维度，就会导致\运算报错或得到极不稳定的解。造成奇异的原因有：
- 量测没有包含足够的信息，比如只有电压幅值量测，没有相角或功率量测；
- 参考节点相角量测缺失，整个相角空间存在一个自由度；
- 支路电流量测方程写错，导致某些状态的偏导项全部为零。
排查思路很简单：算一下rank(H' * R_inv * H)，看是否等于2*n-1。如果不是，就逐类量测检查对应雅可比行是否为空。我在代码里加了一个调试函数，专门打印信息矩阵的零特征值个数，定位缺失信息非常方便。
6.2 迭代不收敛
WLS迭代不收敛通常是初值太离谱，或者量测函数/雅可比矩阵代码有误。我遇到过一次很隐蔽的错误：在直角坐标系和极坐标系之间混用了量测方程，导致雅可比矩阵某个符号错误，结果残差怎么迭代都降不下去。后来把量测函数和雅可比矩阵用数值差分方法做了一致性校验，也就是用(h(x+epsilon)-h(x-epsilon))/(2*epsilon)来对比解析雅可比矩阵元素，几分钟就找到了错误。
另一个值得注意的点是量测残差里如果存在极端坏数据，WLS迭代也可能振荡。解决办法是先做一次残差校验，把残差超过3倍标准差的量测剔除或降权，再重新估计。
6.3 量测相角参考基准
这是一个非常经典的问题。PMU给出的是绝对相角，但状态变量里的相角是相对参考节点的。如果直接把z里的相角量测当成θ_i来构造残差，那么参考节点的相角偏差会污染所有相角估计值。
我采用的做法是：在进入WLS之前，把所有PMU电压相角量测统一减去参考节点的PMU相角量测，转成相对参考节点的相角差；然后把参考节点的相角估计值固定为0。需要注意的是，电流相量量测也要做同样的坐标变换。如果你用Matpower的PMU数据，有时候它内部已经包含了参考处理，多一步转换反而会出错，这时候可以用比较小的测试案例验证对错。
6.4 坏数据与残差检验
状态估计在工程中必须考虑数据完整性。这个项目虽然用的是仿真数据，但我在代码里还是加了一个简单的坏数据检测模块：用WLS估计结果计算残差r = z - h(x_est)，然后标准化残差r_i / sqrt(R_i)，如果超过3就标记为可疑量测。在实际系统中，坏数据可能来自PMU通信中断、幅值饱和或相角跳变，这些不是高斯噪声能描述的，必须靠残差检验去过滤。
我还测试过一种情况：人为把某个PMU电压幅值量测改成偏离真值0.05 p.u.，也就是通常所说的一致错误，标准残差法有时不一定能识别出来，需要借助多个PMU之间的空间相关性来校验。这说明单纯依赖一个算法做坏数据检测是不够的，但在学习阶段，掌握残差检验的基本操作仍然很有价值。
这个项目做到最后，我发现最大的收获不是会调用几个Matlab函数，而是彻底理解了量测方程、雅可比矩阵、权重矩阵这些概念之间如何咬合。调试代码的过程等于把《电力系统状态估计》这本书重新读了一遍。如果你也准备复现这个实验，我建议先不要急着写一整套完整的程序，而是从“只有一个PMU电压幅值量测”的最简场景开始，逐步增加量测类型，每加一类就验证一次收敛性和精度。这样一旦出现问题，你永远知道是刚加入的那部分代码出了问题。
最后分享一个小习惯：所有中间变量都命名得足够直白，注释里写清楚量测类型和单位。状态估计代码最容易在量测索引上出错，比如第几个量测对应哪个节点、是幅值还是相角，索引错一个数字，结果就是天壤之别。宁可多写几行注释，也不要省这点时间。
