做电网状态估计这个课题最折磨人的不是算法推导而是你明明把最小二乘法的公式背得滚瓜烂熟打开MATLAB却连量测矩阵都组织不明白。我第一次在MATLAB里实现加权最小二乘法WLS状态估计时光是一个雅可比矩阵的符号就排查了一下午。这篇博文围绕“基于最小二乘法和快速解耦法的电网状态估计”展开重点讲清楚两件事第一这两种算法的数学逻辑和它们在MATLAB环境下的落地方式第二所谓“节点电压无偏估计”在工程上到底意味着什么以及怎么验证。无论你是做课程设计、本科毕设还是刚接触能量管理系统算法开发的工程师这份实操笔记都值得你对照着调一遍代码。1. 状态估计到底在解决什么问题1.1 一堆SCADA量测数据为什么不能直接拿来用电力系统调度中心每天收到SCADA系统传回来的海量遥测包括节点注入有功、节点注入无功、母线电压幅值、线路两端潮流等。按理说这些数据直接反映了电网当前运行状态那为什么还要专门做一次状态估计答案很朴素这些数据不可信而且互相之间经常矛盾。先看第一个问题量测误差。CT、PT本身有精度等级二次回路有压降数据采集和模数转换也会引入误差所以任何一个遥测值都不是精确值。第二个问题是通信异常通道抖动、设备重启、主备切换都可能让某几个测点刷新出跳变数据。第三个问题是数据冗余引发的矛盾——量测是冗余配置的比如一条线路两端都有潮流测点多个测点之间还存在功率平衡关系一旦误差方向不一致数据之间就会“打架”。直接用这些矛盾数据去算潮流轻则迭代发散重则算出荒谬结果。状态估计做的事情就是在冗余量测的基础上做一次统计处理找出一组最接近“真实状态”的节点电压幅值和相角使得按这组状态算出的量测预测值和实际量测在误差允许范围内尽量接近。电压幅值和相角就是电力系统稳态分析最核心的状态量这也是“节点电压估计”这个标题的由来。1.2 加权最小二乘法凭什么成为行业默认选项加权最小二乘法WLS在状态估计里的地位基本等同于默认选项。它的逻辑一句话就能说清把所有量测误差的平方按权重累加求一组状态让这个累加值最小。权重大的量测相当于告诉算法“这个测点我很信任误差要压得低一点”权重小的量测则允许它偏离预测值远一些。那为什么偏偏选平方误差而不是绝对误差因为量测误差通常被建模为零均值高斯分布。在高斯噪声假设下WLS估计恰好等价于极大似然估计而且是最小方差无偏估计。换句话说在所有线性无偏估计量里WLS的估计误差方差最小。这个性质从理论上保证了只要量测系统没有系统性偏差WLS给出的节点电压估计就是无偏的——估计量的期望等于真值。这正是标题里“节点电压无偏估计”的底气来源。当然工程上也会用到其他准则比如最小绝对值LAV对坏数据更鲁棒但实现复杂度和调试难度都更高。WLS胜在成熟、稳定、好调配合残差检测程序足以应对绝大多数现场需求。我个人的习惯是第一版算法永远先上WLS跑通了再谈优化。1.3 快速解耦法在WLS基础上做减法WLS每次迭代都要重新计算量测函数、雅可比矩阵和增益矩阵其中增益矩阵的因子分解是最大开销。电网节点数动辄几百上千实时性要求高的时候WLS单步迭代的计算量就有点吃不消。快速解耦法就是冲着这个痛点来的。它的思路和潮流计算里的快速解耦法一脉相承。高压输电网中线路电抗远大于电阻有功功率流动主要受两端相角差影响无功功率主要受电压幅值差异影响。也就是说P和θ强相关Q和V强相关而P-V、Q-θ之间的耦合非常弱。基于这一点快速解耦法把WLS问题拆成两个子问题一个只管P和θ一个只管Q和V。然后再进一步把雅可比矩阵固定在平启动点各节点电压1.0标幺、相角差为0上从而把增益矩阵变成常数矩阵只要在迭代前分解一次后续直接复用。代价是迭代次数可能会略增工程上这完全能接受。实践表明正常工况下快速解耦法的收敛解和WLS非常接近节点电压估计值的差异通常在极小范围。但在需要高性能计算的场景里快速解耦法的计算量能比WLS降一个量级实用性很强。2. 两类算法的数学内核与推导逻辑很多人写代码卡壳根因不是循环写错而是公式和代码的对应关系没建立起来。这一章把两类算法的数学链路说透。2.1 量测方程与状态向量的定义方式先定义状态向量。电网稳态分析的状态量是节点电压幅值V和相角θ。由于全网电压相角需要一个参考基准通常把参考节点一般选平衡节点的相角固定为0所以状态向量写成x [θ_2, θ_3, …, θ_n, V_1, V_2, …, V_n]^T注意θ_1没有进状态向量因为它恒为0。这个细节极其关键很多新手把所有节点的θ都放进去增益矩阵必然奇异迭代直接崩掉。量测方程统一写成z h(x) e其中z是所有量测组成的向量h(x)是根据当前状态向量算出的各量测预测值e是量测误差向量。量测类型主要有四种节点注入有功P_i、节点注入无功Q_i、线路首端有功P_ij、线路首端无功Q_ij还有节点电压幅值V_i。每一种都对应一个关于x的非线性函数本质上就是潮流方程。比如节点注入有功的表达式是P_i V_i Σ V_j (G_ij cosθ_ij B_ij sinθ_ij)其中求和遍历所有与节点i相连的节点G_ij和B_ij是节点导纳矩阵对应元素的实部和虚部θ_ij θ_i - θ_j。线路潮流的表达式也类似只是多了一项线路自身导纳的处理。这些公式在构建h(x)时要一条条对应好一个萝卜一个坑。量测误差e的协方差矩阵记作R通常假设各量测相互独立所以R是对角矩阵对角线元素就是各量测方差σ_i^2。权重矩阵W取R的逆方差小、精度高的量测对应权重就大。2.2 WLS迭代求解的完整链路WLS的目标函数是J(x) [z - h(x)]^T W [z - h(x)]求极值令梯度为0得到-H^T W [z - h(x)] 0其中H ∂h/∂x是雅可比矩阵行数等于量测数列数等于状态数。由于h(x)是非线性的这个方程没有闭式解需要用高斯-牛顿法迭代求解。在第k次迭代点x^k处对h(x)做一阶泰勒展开代入梯度方程整理后得到核心迭代公式G Δx H^T W Δz其中G H^T W H是增益矩阵Δz z - h(x^k)Δx x^(k1) - x^k。每次迭代要重新计算H和G然后对G做Cholesky分解或LDL分解回代求解线性方程组。收敛后输出最优状态估计值。这里面有几个容易出错的点。一是雅可比矩阵的表达式尤其是线路潮流对两端电压幅值和相角的偏导正负号很容易搞混。二是状态向量的排列方式和雅可比矩阵的列索引如果状态是“先θ后V”的结构那么H的列也必须严格对齐否则矩阵相乘直接报错或者算出错误结果。2.3 快速解耦法的两个关键简化快速解耦法相比WLS核心就是两个简化。第一个简化是解耦。把雅可比矩阵H按量测类型和状态类型分块H [H_Pθ, H_PV; H_Qθ, H_QV]其中H_Pθ是有功量测对相角的偏导H_PV是有功量测对电压幅值的偏导H_Qθ是无功量测对相角的偏导H_QV是无功量测对电压幅值的偏导。快速解耦法直接忽略H_PV和H_Qθ两个耦合块把问题拆成两个独立子问题有功-相角子问题用有功量测估计相角θ电压幅值固定不动无功-电压子问题用无功量测和电压幅值量测估计V相角固定不动。两个子问题交替迭代直到状态修正量都小于阈值。第二个简化是常数化。在平启动点V_i 1.0θ_ij 0计算子雅可比矩阵然后得到两个常数增益矩阵B H_Pθ^T W_P H_Pθ B H_QV^T W_Q H_QV由于平启动点上的导纳关系非常规整B和B实际上只与电网导纳结构和量测权重配置有关和当前运行状态无关。所以迭代之前只需要做一次矩阵分解之后每一轮迭代就是反复求解两个固定系数的线性方程组B Δθ H_Pθ^T W_P [z_P - h_P(V, θ)] B ΔV H_QV^T W_Q [z_Q - h_Q(V, θ)]这里的H_Pθ和H_QV在快速解耦实现中也是常数矩阵右侧乘积可以预计算一部分单步迭代计算量大幅下降。本质上快速解耦法是用模型简化换计算速度收敛解和WLS非常接近但代码结构差别很大。3. MATLAB环境下的工程实现3.1 数据准备从IEEE算例到量测组织写代码之前先把数据准备好。我习惯用MATPOWER的IEEE 14节点系统当测试对象。MATPOWER提供了case14.m文件里面有完整的母线矩阵和支路矩阵用loadcase命令直接加载。如果不想依赖MATPOWER也可以手动定义bus和branch两个矩阵但自己写容易出错建议直接用成熟数据。数据准备分三个层次。第一层是网络参数也就是节点导纳矩阵Ybus。用MATPOWER自带的makeYbus函数可以一步生成。手写的话要遍历支路表把每条支路的串联阻抗转成导纳加上对地电纳再按节点汇聚。这里有个大坑变压器支路的变比处理变比非1的支路在形成导纳矩阵时自导纳和互导纳都要乘上变比系数漏掉这个系数后面全盘皆错。第二层是真值工况。做状态估计算法验证必须有一个“真实状态”当基准。我的做法是先用MATPOWER的runpf函数跑一次潮流把潮流解作为真值状态V_true和θ_true然后由真值状态计算各量测通道的精确值作为“无噪声量测”再叠加高斯白噪声模拟真实SCADA量测。第三层是量测结构。定义量测对象每个量测点需要记录四个信息量测类型P注入、Q注入、P线路、Q线路、电压幅值、位置注入量测记录节点号线路量测记录首末端节点号、量测数值、权重。权重一般取量测方差的倒数典型设置是注入功率方差0.01、线路潮流方差0.01、电压幅值方差0.0001。电压量测精度高权重就大很多。3.2 核心代码框架与关键函数拆解直接给可参考的核心代码框架。我按模块拆开讲。第一个模块是量测函数计算。输入当前状态V、θ和量测结构输出所有量测的预测值向量h。以节点注入有功为例我喜欢用复功率的写法代码简洁、不易出错function hp compute_P_injection(V, theta, Ybus, bus_list) n length(V); hp zeros(length(bus_list), 1); Vc V .* exp(1j * theta); I Ybus * Vc; S Vc .* conj(I); for k 1:length(bus_list) i bus_list(k); hp(k) real(S(i)); end end用复功率计算的好处是不用把三角函数的偏导公式一个个展开雅可比矩阵的解析计算也可以用类似技巧。不过更稳妥的做法是列出每个量测类型关于θ_j和V_j的偏导公式再填充矩阵。对状态估计这种精度敏感的算法解析雅可比比数值差分可靠速度也快得多。第二个模块是WLS主循环function [V, theta, iter] wls_estimator(Ybus, meas, W, max_iter, tol) n length(Ybus); V ones(n, 1); theta zeros(n, 1); state [theta(2:end); V]; % θ_1不参与估计 for iter 1:max_iter theta_full [0; state(1:n-1)]; V_full state(n:end); [h, H] calculate_measurements(V_full, theta_full, Ybus, meas); G H * W * H; dx (G 1e-8 * eye(size(G))) \ (H * W * (meas.z - h)); state state dx; if max(abs(dx)) tol break; end end theta [0; state(1:n-1)]; V state(n:end); end这里把参考节点相角从状态向量里排除了避免增益矩阵奇异同时给G加了对角微量扰动进一步提高数值稳定性。第三个模块是快速解耦法主循环function [V, theta, iter] fd_estimator(Ybus, meas, Wp, Wq, max_iter, tol) n length(Ybus); V ones(n, 1); theta zeros(n, 1); % 在平启动点计算常数雅可比和增益矩阵 [Hp_fl, Hq_fl] calculate_fd_jacobian(V, theta, Ybus, meas); Bp Hp_fl * Wp * Hp_fl; Bq Hq_fl * Wq * Hq_fl; % 预先做LU分解 [Lp, Up] lu(Bp); [Lq, Uq] lu(Bq); for iter 1:max_iter % 有功-相角迭代 hp compute_active_meas(V, theta, Ybus, meas); dp Hp_fl * Wp * (meas.zp - hp); dtheta Up \ (Lp \ dp); theta theta dtheta; % 无功-电压迭代 hq compute_reactive_meas(V, theta, Ybus, meas); dq Hq_fl * Wq * (meas.zq - hq); dV Uq \ (Lq \ dq); V V dV; if max(abs([dtheta; dV])) tol break; end end end快速解耦法里量测被拆成有功组和无功组分别组织成meas.zp和meas.zq权重矩阵Wp、Wq也分别定义。迭代先更新相角再更新电压交替进行。LU分解只做一次每次迭代就是解两个三角方程组计算量小很多。3.3 算例结果如何判断估计结果是否“无偏”代码跑通后要直面一个问题怎么判断估计结果是无偏的严格定义的无偏指估计量的期望等于真值。但单次运行只得到一个样本谈不上期望。所以验证无偏性必须做蒙特卡洛仿真。我的做法是固定真值工况生成1000组高斯噪声叠加到同一组量测上分别跑状态估计得到1000组V和θ的估计值。然后对每个状态量计算均值和真值对比。如果均值与真值的偏差显著小于单个量测噪声的标准差且偏差的正负号随机就可以认为无偏性成立。可以用偏差百分比来量化bias_i mean(V_est_i) - V_true_i在IEEE 14节点系统里比如节点3的电压幅值真值是1.0234标幺1000次蒙特卡洛仿真后估计均值如果是1.0233偏差0.0001基本就是无偏。而某次单次估计可能是1.0198或者1.0279这都属于正常波动不能说明有偏。还可以把所有节点的估计值画成箱线图。箱体中心线应该贴近真值线箱体宽度反映有效性。我对比过WLS和快速解耦法的蒙特卡洛结果两者估计均值几乎重合差异主要在单次迭代耗时上——快速解耦法的平均计算量比WLS低60%以上在更大规模网络上优势更明显。但还有一个工程上容易忽略的前提无偏性成立的前提是量测误差均值为零。如果某条线路上CT存在角差或变比误差等于给量测叠加了非零均值误差那么无论WLS的统计性质多好估计结果必然有偏。所以“无偏估计”不只是算法问题更是量测系统质量问题。4. 实际调试中的问题与解决实录4.1 迭代发散最让人头疼的问题跑状态估计最常见的挫败感就是迭代不收敛或者收敛到离谱的点。我归纳了几个高频原因。第一雅可比矩阵符号错误。这是新手最高频的坑。线路有功潮流对θ_i和θ_j的偏导符号相反搞反了迭代方向就错。排查方法很简单把解析雅可比矩阵和数值差分做一次对照比较两者结果。我在调试阶段一定会做这个检查别嫌多花五分钟。第二参考节点相角处理不对。状态向量里如果还留着θ_1增益矩阵一定奇异因为所有量测函数对θ_1的共同平移不变。解决方式就是我前面写的把θ_1从状态向量里剔除恢复时再拼回去。或者用一个非常大的权重把θ_1钉死在0但不如直接剔除干净。第三量测权重差异过大导致病态。电压幅值量测方差可能只有1e-4注入功率方差是1e-2两者相差两个数量级增益矩阵条件数会变得很大。这时候可以加正则项或者检查量测权重设置是否合理。注意不要为了数值稳定把权重乱调那会破坏估计的统计特性。第四初值距离真值太远。状态估计整体上不像潮流计算那样对初值敏感但在重负荷工况下从平启动开始个别情况也可能振荡。解决方式是用潮流计算结果当初值或者先跑几步确认残差在下降再继续。4.2 坏数据的检测与处理工程现场的SCADA数据里坏数据是常态。所谓坏数据是严重偏离真值的量测通常来自通信误码、连接器松动或设备故障。WLS对坏数据特别敏感因为平方误差会把坏数据的异常偏差放大。经典检测方法是残差分析。先定义量测残差r z - h(x̂)在WLS框架下残差的协方差矩阵可以算出来Ω R - H G^{-1} H^T标准化残差就是r_i^N |r_i| / sqrt(Ω_ii)当某个量测的标准化残差超过阈值一般取3.0对应约99.7%置信区间就怀疑它是坏数据。实际调试中我的流程是先跑一次状态估计计算所有量测的标准化残差找出最大的那个如果超过阈值剔除它重新估计重复这个过程直到所有标准化残差都低于阈值。有个工程细节值得注意多个坏数据可能互相掩盖。比如两条相邻支路同时出现坏数据残差可能被压低单个量测的标准化残差都不超阈值但估计结果已经失真。这种情况靠单次残差检测往往查不出来需要更复杂的多坏数据辨识算法。不过在课程设计和入门阶段逐步剔除再估计已经够用。4.3 MATLAB编程层面的几个典型坑最后聊几个MATLAB环境下的实操细节都是我在调试中踩过或者帮别人排查过的。中文注释乱码问题。MATLAB 2023及更新版本默认用UTF-8读取脚本而旧版脚本可能是GBK编码保存的打开就会出现乱码。解决方式是在MATLAB预设项里把字符编码改成UTF-8或者用编辑器把旧脚本重新另存为UTF-8。这个和状态估计本身无关但确实很影响调试心情。矩阵维度问题。如果用state [theta(2:end); V]这种排列计算雅可比矩阵时列索引必须和state一一对应。我见过不少代码在算H时用了完整状态索引最后G的维度不匹配。建议在写量测计算函数时把状态排列方式写成注释放在文件头部并且用变量名明确区分完整状态和估计状态。矩阵求解方式。MATLAB里解G * dx rhs直接写dx G \ rhs别用inv(G) * rhs。左除在数值稳定性和速度上都优于显式求逆尤其当G是稀疏矩阵时左除走的是稀疏直接求解器内存占用和计算时长差异巨大。同理快速解耦法里B和B的LU分解只做一次别在循环里反复分解。尽量把量测组织成结构体不要用多个分散的向量。比如meas.p_inj_node、meas.q_inj_node、meas.p_branch_from这些字段计算h(x)和H时可以按类型分组处理可读性和可维护性都会好很多。我第一版代码就是所有量测揉在一个大向量里每次加一种量测类型就要改一遍索引非常痛苦。5. 从仿真到工程状态估计代码还能怎么扩展5.1 大规模系统的计算加速技巧在MATLAB里用14节点小算例验证正确性只是第一步真到实际电网几百上千个节点原来的代码可能会慢得离谱。我接手过省级电网模型后有几个加速心得。首先是稀疏化。MATLAB很多内置函数天然支持稀疏矩阵但如果你构建雅可比矩阵或增益矩阵时用的是全矩阵内存和时间都会爆炸。用sparse函数显式声明稀疏结构G的因子分解速度能提升一个数量级。另一个技巧是选用合适的稀疏求解器——MATLAB左除会自动选择但你可以利用矩阵的对称正定性用ichol做预条件共轭梯度法大规模情况下有时更快。其次是量测预过滤。实际电网里很多量测类型相似、位置相近可以在进状态估计前做等价合并减少量测方程数目。这个做法在工程上很常见对计算效率提升明显。第三是并行化。蒙特卡洛无偏性验证这种跑上千次的场景用parfor循环直接利用多核并行几乎零成本提速。我之前把无偏性验证从串行改成parfor1000次仿真时间从30分钟缩短到5分钟。5.2 无偏性之外鲁棒性和动态估计聊完无偏性想再延伸一层。WLS在量测误差严格高斯、没有坏数据时表现完美但真实电网里非高斯噪声和坏数据是常态。所以工业级状态估计器往往在WLS基础上加入不良数据辨识和鲁棒化处理比如用指数权函数迭代降权或者直接上最小绝对值估计。另一个进阶方向是动态状态估计。传统静态估计用一帧量测做一次计算而动态估计利用系统动态模型比如发电机惯性、负荷变化规律把时间维度引进来用卡尔曼滤波类方法递推估计状态。这样既能提前滤波也能给出比静态估计更平滑的结果。我在后续项目里实践过实现复杂度明显更高但在处理PMU量测时优势很大。把WLS和快速解耦法吃透之后这些都是很自然的下一步。最后说点个人体会。这个课题做下来我最深的感触是公式推导和代码实现之间的距离往往不在算法本身而在工程细节。参考节点的处理、雅可比矩阵的符号、权重矩阵的组织方式每一个细节都可能卡住你半天。把WLS在MATLAB里完整跑通是一段很值得的修炼它会逼着你把每个公式都转成具体的矩阵操作。快速解耦法则让你直观理解“用简化模型换计算速度”这种工程思想。至于节点电压无偏估计记住一句话量测通道不干净再漂亮的算法也给不了你无偏的结果。希望这篇整理能帮你少踩几个坑。