做电力系统状态估计的同学一定绕不开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这对组合状态估计算法有很多种但WLSWeighted 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^(k1) 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(n1: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(n1:2*n); if norm([dP; dQ], inf) tol converged true; return; end end converged false; end这段代码的雅可比矩阵build_powerflow_jacobian是潮流计算的关键公式书中都有不在这里展开。需要提醒的是NR潮流对初值比较敏感对IEEE 14节点系统用平启动V1, θ0基本都能收敛如果一个算例不收敛先检查节点类型和参数是否有误不要急着调初值。4.4 结果输出与误差统计估计完成后我习惯用一张图把结果可视化包含四个子图节点电压幅值的真值、估计值、NR潮流值对比节点相角的真值、估计值、NR潮流值对比电压幅值估计误差的柱状图相角估计误差的柱状图。同时打印误差统计表格格式类似下面这样误差项电压幅值 (p.u.)相角 (deg)最大绝对误差0.00230.031平均绝对误差0.00080.012均方根误差0.00110.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(xepsilon)-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电压幅值量测”的最简场景开始逐步增加量测类型每加一类就验证一次收敛性和精度。这样一旦出现问题你永远知道是刚加入的那部分代码出了问题。最后分享一个小习惯所有中间变量都命名得足够直白注释里写清楚量测类型和单位。状态估计代码最容易在量测索引上出错比如第几个量测对应哪个节点、是幅值还是相角索引错一个数字结果就是天壤之别。宁可多写几行注释也不要省这点时间。
