做电力系统分析的人大概都跟这个问题打过照面调度大厅里看到一屏的电压、有功、无功你敢直接拿它们做决策吗不敢。互感器有误差采集终端有量化误差通道上还可能丢帧哪怕噪声不大叠加起来也足以让一些薄弱断面被误判。这就是为什么调度自动化系统里必须有一层状态估计——它不是简单“测”一个值而是把全网冗余的量测数据拿来做统计拟合给出运行状态的最可信估计也就是各节点的电压幅值和相角。加权最小二乘法WLS是干这活几十年的经典主力。同步相量测量单元PMU这几年在工程现场越来越多它最大的本事是依靠统一时钟直接测出电压相量的幅值和绝对角度精度和更新速率都远超传统SCADA量测。量测手段变了估计算法的效果自然会变。于是大家关心一个很实际的问题把PMU量测喂给WLS估计器电压幅值和角度的估计精度到底能提升多少跟使用Newton-Raphson潮流计算出来的“标准状态”放一起比又能差多少本文就把这件事从头到尾做一遍用Matlab实现WLS状态估计并接入PMU量测同时用Newton-Raphson潮流求出系统状态作为基准真值最后比较两者的电压幅值和相角误差。文章适合正在做电力系统分析课程设计、在研究状态估计算法、或者想快速搭一套完整仿真流程来验证PMU效果的工程师和同学。整套流程在仿真里闭环不需要依赖任何商业软件纯手写代码就能跑出结果。1. 项目背景与整体设计思路1.1 状态估计从带噪声的量测里“拟合”出可信状态先解释为什么状态估计非做不可。一个n节点电力系统运行状态可以用节点电压幅值V_i共n个和相角θ_i共n-1个参考节点角度固定为0来描述。调度中心拿到的原始量测是SCADA采集的母线电压幅值、支路潮流、节点注入功率等这些量测有三个先天毛病一是数量冗余同一运行状态可以被很多量测反复描述二是含有误差测量误差、变送器误差、通信传输误差都混在其中三是分布不均有些母线量测特别全有些则几乎没有。直接把这些带噪声的量测当作“真值”来用轻则分析结果不可靠重则把调度员的判断带偏。状态估计的核心思路是把量测和状态之间的物理关系写成方程组然后利用冗余量测按照最小二乘准则求解出最符合所有量测的电压幅值和相角。通俗点说现场一堆“靠谱程度不同”的传感器在七嘴八舌描述同一个系统状态估计器做的就是加权投票误差小的传感器话权大误差大的话权小最后给出一个大家都能接受的状态共识。WLS把这个问题变成最优化找一个状态向量x使目标函数J(x)最小。J的含义是全部量测残差平方的加权和权重W取每个量测误差方差的倒数。这个过程不需要系统动态模型只需要量测函数和权重矩阵所以它特别工程化也特别稳定。1.2 PMU直接把相角量测“端上桌”的新队友传统SCADA时代有个尴尬现实电压幅值可以测有功无功潮流也可以测但相角几乎没法直接测。调度系统里的相角基本是靠状态估计间接“推算”出来的。PMU的出现改变了这一点它利用统一时钟同步采样在同一个时间基准下测量电压和电流的相量也就是同时给出幅值和绝对相角。相角量测不再是“推算值”而是一个和电压幅值同级别的直接量测。PMU给状态估计带来的好处最直观的可以总结为两条。一是角度信息直接可用量测方程中电压幅值和角度与状态量的关系变得线性或近线性收敛难度大大下降。二是精度高典型PMU的电压幅值测量精度在0.1%量级相角测量精度能达到0.01到0.05度比传统SCADA高出约一个数量级。在WLS的加权框架下误差方差小意味着权重高PMU量测会在估计中占据更大话语权从而把整体估计精度拉起来。当然PMU不是没有代价。它的数据是高速流式的需要和SCADA慢速数据做时间对齐它测的是相量工程上还要考虑互感器幅值和相位误差、卫星信号丢失等实际问题。仿真中这些问题可以先简化为理想的加噪模型但心里得有数不能把仿真结论直接照搬到现场。1.3 为什么拿WLS状态估计和NR潮流结果硬碰硬很多第一次接触这个题目的人会疑惑状态估计和潮流计算明明是两个不同流程为什么要把它们放在一起比我在仿真里的做法是这样的先用Newton-Raphson潮流在某一组确定的负荷和发电出力下解出全网电压幅值和相角这个解满足潮流方程且没有任何测量噪声我在仿真中把它当作“系统真实状态”。然后基于这个真实状态模拟出PMU和SCADA量测给量测加上符合实际精度的随机噪声。接下来用WLS从这些带噪声的量测中反推状态。最后把WLS的估计结果与NR潮流的真值比较统计电压幅值和相角误差。这个“潮流真值→加噪量测→状态估计→误差统计”的流程本质上是蒙特卡洛验证估计器性能的标准做法。它的合理性在于NR潮流是精确代数方程的解不依赖测量设备所以在仿真环境中可以做“裁判”。真实系统里永远拿不到真值只能用估计值去贴近量测但仿真环境给了我们一个上帝视角这正是用来检验算法精度的黄金机会。我再多提醒一句这种比较的结论只说明“在该运行点、该噪声模型下WLS的估计贴近潮流解的程度”不要把它理解成“WLS能替代潮流计算”。2. 核心算法原理拆解WLS和NR的数学本质2.1 WLS状态估计让加权残差平方和说话WLS的出发点很朴素。设系统的状态向量为x包含所有节点的电压幅值以及除参考节点外的相角总量是2n-1。量测向量z包含所有量测量测函数h(x)描述“状态x下这些量测应有的取值”。残差就是z - h(x)。WLS的目标是让下面的二次型最小J (z - h(x))ᵀ W (z - h(x))其中W是对角权重矩阵对角线元素是量测误差方差的倒数也就是1/σᵢ²。权重越大代表这个量测越可信估计结果会越照顾它的残差。h(x)对状态量是非线性的因为电压幅值、支路潮流这些量测表达式里都含sin、cos项所以WLS要用高斯-牛顿法迭代求解。每一轮迭代在当前状态点x处把h(x)线性化求量测函数的雅可比矩阵H然后解下面的正规方程Δx (Hᵀ W H)⁻¹ Hᵀ W (z - h(x)) x x Δx其中G HᵀWH称为增益矩阵。迭代到Δx足够小或目标函数不再下降为止。整个过程很像“用泰勒展开反复修正猜测值”每轮都沿着加权梯度的方向走一步直到残差平方和不再显著变化。2.2 量测函数与雅可比矩阵的构建要点实现WLS的关键是把每一类量测对应的h(x)和H写对。最常用的有下面几类。SCADA电压幅值量测h Vᵢ。对应H的行只有对状态量中Vᵢ的偏导为1其他为0。非常简单这也是状态估计里最直观的一行。SCADA节点注入功率以节点i注入有功为例h Vᵢ Σⱼ Vⱼ(Gᵢⱼcosθᵢⱼ Bᵢⱼsinθᵢⱼ)其中θᵢⱼ θᵢ - θⱼG和B取自导纳矩阵的实部和虚部。偏导公式有四组对θᵢ、θⱼ、Vᵢ、Vⱼ分别求导这部分是H矩阵里最容易出错的地方。我的习惯是先手写解析式再用小步长数值差分校验一遍确认无误再固化到代码里。SCADA支路潮流量测支路i-j的有功可以写成Pᵢⱼ Vᵢ²g - VᵢVⱼ(g cosθᵢⱼ b sinθᵢⱼ)其中g和b是支路导纳的实部和虚部。无功量测的表达式类似只是cos和sin的位置以及前面的系数不同。偏导同样按四个变量展开公式虽然繁琐但很有规律注意符号别搞反就行。PMU电压相量量测这是线性量测h [Vᵢ; θᵢ]。对应H的行是稀疏的常数对Vᵢ的偏导为1对θᵢ的偏导为1其余为0。这也是为什么PMU一加入增益矩阵的条件数会明显改善估计收敛更快。如果PMU还带电流相量量测电流和电压相量在复数域里也是线性关系整行量测依然是关于状态量的线性函数只是把实部虚部分别拆成两个量测方程而已。混合量测下H矩阵由线性块和非线性块拼装而成。第一次实现时如果H矩阵的拼接顺序和z的排列顺序不一致迭代很容易直接发散。建议在代码里给每个量测定义一个结构体包含类型、首末端节点、权重、量测值然后用统一接口计算h和H避免手滑。2.3 Newton-Raphson潮流解功率平衡方程NR潮流解决的问题是已知各节点的注入有功和无功求全网电压幅值和相角。它要解的方程是节点功率平衡方程Pᵢ Vᵢ Σⱼ Vⱼ(Gᵢⱼcosθᵢⱼ Bᵢⱼsinθᵢⱼ) Qᵢ Vᵢ Σⱼ Vⱼ(Gᵢⱼsinθᵢⱼ - Bᵢⱼcosθᵢⱼ)把失配量F定义为给定注入功率与当前状态计算出的注入功率之差即F [ΔP; ΔQ]其中ΔP针对PV和PQ节点的有功失配ΔQ只针对PQ节点的无功失配。然后求F对状态量的雅可比矩阵J解方程J Δx -F更新状态后重复迭代直到失配量小于阈值。J的结构与WLS的H完全不同J是方阵因为方程数等于未知数而WLS的H是长方形因为量测数远大于状态数。NR潮流的收敛特点是二阶收敛初值靠近解时收敛非常快通常几轮迭代就能达到1e-10精度。但初值离解太远时可能振荡甚至发散。所以在做基准真值生成时我会把平启动值设成所有节点电压幅值1.0、相角0再配合合理的节点归类通常不会有问题。这里有一个容易漏的细节如果发电机无功出界PV节点要转成PQ节点重新计算否则迭代可能在越限区间反复震荡。2.4 两种方法的核心区别对比维度WLS状态估计NR潮流计算要解决的问题从含噪声冗余量测中拟合运行状态从给定功率注入中精确求解运行状态输入PMU/SCADA量测 权重矩阵负荷、发电机出力、PV节点电压设定值方程维度量测数远大于状态数超定方程数等于未知数方阵核心数学工具加权最小二乘 高斯-牛顿迭代牛顿法解非线性代数方程组输出电压幅值和相角的最优估计电压幅值和相角的理论精确解对噪声的态度容忍噪声通过冗余抑制误差不考虑噪声功率注入即真值工程角色EMS/调度自动化核心功能运行方式计算、离线分析核心功能这个表能帮初学者快速分清两个概念。二者不是竞争关系而是互补关系潮流计算告诉我们系统在给定运行方式下应该是什么状态状态估计告诉我们系统在当前量测下最可能是什么状态。做仿真验证时前者常被拿来当后者的“标准答案”用。3. Matlab实现全流程从潮流量测到误差对比3.1 测试系统建模与数据准备我这次用的是某14节点测试系统节点数适中支路20条左右既有PQ节点又有PV节点还带一个平衡节点足够把WLS和NR的细节暴露出来。数据准备阶段主要做两件事构建导纳矩阵Ybus以及明确节点类型索引。构建Ybus时按支路逐条处理每条支路对非对角元素减去串联导纳y对两端节点自导纳加上y如果是带对地导纳的支路还要在对角元素上加上对应的半电容导纳。代码非常简单核心循环就十几行。但我习惯在构建好之后把对角线元素和互动关系打印出来肉眼检查一遍防止原始数据里单位不一致这类低级错误。节点类型索引要分成三组平衡节点ref、PV节点数组、PQ节点数组。function Ybus build_ybus(branch, nb) Ybus zeros(nb, nb); for k 1:size(branch, 1) f branch(k,1); t branch(k,2); r branch(k,3); x branch(k,4); bs branch(k,5); y 1/(r 1j*x); Ybus(f,f) Ybus(f,f) y 1j*bs; Ybus(t,t) Ybus(t,t) y 1j*bs; Ybus(f,t) Ybus(f,t) - y; Ybus(t,f) Ybus(t,f) - y; end end这个函数不复杂但它是后面所有计算的基石。如果Ybus错了NR潮流和WLS估计会同时出错而且错误很难从结果里一眼看出来所以一定要在前期卡好质量。3.2 NR潮流生成基准真值接下来实现NR潮流。核心代码是失配量计算和雅可比矩阵求解。收敛判据我用失配量的最大绝对值一般设到1e-10。这个精度对后面的“真值”非常重要——如果潮流只收敛到1e-6那么它本身的误差就会淹没后面要观察的WLS估计误差。潮流迭代核心如下function [V, theta, iter] nr_solve(Ybus, Sbus, V0, ref, pv, pq, tol, max_iter) nb length(V0); order [pv; pq]; % 待求解相角的节点 np length(order); % 相角未知量个数 V V0; theta zeros(nb, 1); for iter 1:max_iter S V .* conj(Ybus * V); dP real(Sbus(order)) - real(S(order)); dQ imag(Sbus(pq)) - imag(S(pq)); F [dP; dQ]; if norm(F, inf) tol break; end J build_jacobian(V, theta, Ybus, ref, pv, pq, order); dx -J \ F; % 注意负号F定义为给定值-计算值 theta(order) theta(order) dx(1:np); V(pq) V(pq) .* (1 dx(np1:end)); end end雅可比函数build_jacobian把J分成四个大块dP/dθ、dP/dV只对PQ节点电压、dQ/dθ、dQ/dV只对PQ节点电压。这里有一个容易忽略的细节雅可比矩阵里的电压幅值修正如果采用ΔV/V形式和采用ΔV形式的偏导公式会差一个V倍。我在代码里统一用ΔV形式偏导公式配套省得调试时怀疑人生。另一个细节是dx -J \ F的负号失配量定义为“给定值减去计算值”时牛顿修正量方向必须反转。这个负号我丢过一次结果迭代两三步直接发散。3.3 量测模拟给真值加噪声有了潮流真值V_true和theta_true接下来就可以模拟量测了。这一步的核心是知道PMU和SCADA各自的精度用高斯白噪声模拟实测误差。SCADA量测我做四类母线电压幅值、支路有功潮流、支路无功潮流、节点注入有功。通常还会加入注入无功但不一定所有母线都有注入量测要按可观测性要求配置。PMU量测做两类电压幅值相量和电压相角。仿真中可以指定若干关键枢纽母线加装PMU。PMU误差方差比SCADA小一个量级具体来说SCADA电压幅值标准差约0.005到0.01pu对应0.5%到1%的误差。SCADA潮流标准差约0.01到0.03pu对应1%到3%的误差。PMU电压幅值标准差约0.0005到0.001pu对应0.05%到0.1%。PMU相角标准差约0.0002到0.001rad对应约0.01到0.05度。量测向量在代码里可以这样生成每一种量测都叠加正态分布随机数z zeros(num_meas, 1); sigma zeros(num_meas, 1); % 假设第k个量测是某母线的PMU电压幅值 z(k) V_true(bus_k); sigma(k) 0.0005; % 假设第k1个量测是同一母线的PMU相角 z(k1) theta_true(bus_k); sigma(k1) 0.0002; % 加噪声 z(k) z(k) sigma(k) * randn; z(k1) z(k1) sigma(k1) * randn;需要注意randn每次运行结果不同所以要做多次蒙特卡洛试验然后取统计误差别拿单次结果下结论。权重矩阵W就是对角线元素放1./sigma.^2。我会先把所有sigma存好再一次性构建W这样后面改动PMU配置时只需要改sigma的分配逻辑。还有一个坑参考节点相角为0时PMU角度量测必须在同一坐标系下真值也基于该参考否则两者之间会出现一个伪偏差这一点我后面在问题排查部分还会展开。3.4 WLS状态估计核心迭代WLS的状态向量排列我采用这样的方式把所有非参考节点的相角放在前面所有节点的电压幅值放在后面总共2n-1个状态量。这一步没有PV/PQ的区别所有电压幅值都是状态这正是状态估计和潮流建模思路的重要差别。迭代代码大致如下nb length(V0); x [theta0(2:end); V0]; % n-1个相角 n个电压幅值 for iter 1:max_iter hx meas_func(x, branch, Ybus, meas_type); H jacob_meas(x, branch, Ybus, meas_type); r z - hx; G H * W * H; dx G \ (H * W * r); x x dx; if norm(dx, inf) 1e-6 break; end end % 还原电压幅值和相角 theta zeros(nb, 1); theta(2:end) x(1:nb-1); V x(nb:end);meas_func需要根据每个量测的类型去计算对应的h(x)。我习惯用switch-casecase 1是电压幅值量测case 2是PMU相角量测case 3是支路有功case 4是支路无功case 5是节点注入有功case 6是节点注入无功。每个case里写对应的h和H的行向量。这里要特别说下收敛判据。我习惯用Δx的无穷范数小于1e-6同时检查目标函数J是否在下降。如果J反而增加说明H矩阵或者权重有问题要停下来查量测类型和索引。还有一种情况是迭代来回震荡那基本就是初值距离解太远或者PMU角度量测的参考没对齐。WLS在状态估计中不需要像潮流那样做PV/PQ类型切换因为所有节点电压幅值都被当状态量处理发电机的无功出力和PV节点的电压设定值在状态估计的模型里通常是作为等式约束或伪量测处理不直接出现在状态向量中。3.5 结果对比与误差指标跑完NR和WLS以后就是比较环节。我做的误差指标有三个电压幅值最大绝对值误差、相角最大绝对值误差、所有节点电压幅值和相角的均方根误差RMSE。RMSE的计算公式是sqrt(mean((V_est - V_true).^2))角度同理。如果需要体现PMU的增益可以把同一套潮流真值分别跑“只有SCADA量测”的WLS和“SCADAPMU量测”的WLS对比两条误差曲线。这种对照在论文里常见实操里也最能看出PMU的价值。画图的时候我习惯把每个节点的电压幅值估计误差画成柱状图相角误差画成另一个子图。柱状图能一眼看出误差集中在哪些母线尤其能暴露“离参考节点远的角度误差大”这个典型现象。RMSE数值则在图注里给出。误差统计代码很简单err_V abs(V_est - V_true); err_theta abs(theta_est - theta_true); rmse_V sqrt(mean(err_V.^2)); rmse_theta sqrt(mean(err_theta.^2)); figure; subplot(2,1,1); bar(err_V); ylabel(电压幅值误差(pu)); subplot(2,1,2); bar(err_theta * 180/pi); ylabel(相角误差(度));实验结论通常是这样加PMU以后角度误差显著下降因为传统SCADA完全无法直接测量相角角度主要靠支路潮流间接推算误差会被放大PMU直接把角度“读”出来误差自然小一个量级。电压幅值因为SCADA本来就能直接量测改善幅度相对有限但当PMU数量足够多时对关键母线的幅值约束也会明显增强。4. 实操中常见问题与排查增益奇异、坏数据、参考角4.1 增益矩阵奇异或接近奇异这是WLS里最常见的“啪一下死了”的故障。G奇异意味着状态量中至少有一个方向上没有任何量测信息换句话说系统不可观测。常见原因量测配置不全某些母线既没有本地量测也没有相连支路量测或者参考节点角度方向被重复约束。排查办法是在迭代前算一下rank(G)和cond(G)如果秩缺了就打印H矩阵看看哪一列全零。如果遇到接近奇异但不是完全奇异的情况可以给HWH加一个很小的对角扰动比如1e-6倍单位阵这在实际工程中就是Levenberg-Marquardt修正。如果加了PMU仍然奇异多半是PMU安装位置太集中比如全装在同一区域另一侧母线仍然没有量测覆盖。这说明PMU配置不是简单地“越多越好”还需要考虑覆盖范围。4.2 初值设置不合理导致不收敛WLS的迭代本质和NR潮流一样是局部收敛算法初值离解太远就发散。常见的错误是把初值设置成全0电压幅值为0会导致复功率公式里出现除以零问题或者雅可比矩阵奇异。我在实际中统一用平启动所有节点电压幅值1.0相角0。这个初值在绝大多数运行点下都能收敛。如果遇到重负荷系统平启动也可能失败那就分两步走先用直流潮流或一个粗略的高电压估计做预热再切到WLS精确迭代。另外提一个容易被忽略的点如果角度量测的参考和状态量的参考不一致初值等于白设。我在仿真中统一把参考节点角度固定在0PMU角度量测真值也基于这个参考坐标系生成两边天然对齐。4.3 坏数据会把估计结果带偏多少WLS的一个短板是对坏数据敏感。因为目标函数是残差平方和一个大误差量测的残差平方会以平方的速度放大影响直接把状态拉偏离。仿真里我故意在某个支路有功量测上加了10倍标准差的大误差很快看到该支路端节点的电压估计值明显偏离真值。这提醒我们PMU虽然精度高但如果某个PMU通道坏了且没有被识别它的“高权重”反而会放大坏数据的影响——因为权重高残差会被压得极低污染更隐蔽。排查坏数据的实用手段是残差分析计算残差r z - h(x̂)然后用每个量测的标准差标准化得到标准化残差。标准化残差超过3.0的量测大概率是坏数据。识别出来后把这路量测权重置零或剔除再重新做WLS。这个“估计-检验-剔除-再估计”的循环是工程上最常见的坏数据处理流程。实际项目里这一步比算法本身更花时间。4.4 参考角不一致带来的虚假误差我在不少入门复现里看到过这样的情况WLS估计出的角度和NR潮流计算出的角度差值图上所有节点都偏了一个常数角。这通常不是算法错了而是两者的参考坐标系不一致。NR潮流把某节点角度固定为0WLS如果把另一个节点固定为0或PMU角度量测采用了绝对坐标系两套结果就会整体平移一个角度。解决方法是在比较前做一次角度平移计算参考节点处的角度差把它从所有估计角度里减掉然后再画误差图。这一步虽然简单但处理不好会让人浪费很多时间查代码。我通常是这样处理theta_est theta_est - theta_est(ref) theta_true(ref);顺便把参考节点的绝对误差变成0剩下的节点误差才是真正的估计质量体现。为了方便以后快速定位问题我把调试中反复遇到的现象和原因整理成一张速查表现象可能原因处理办法G奇异或rank缺量测不足、覆盖不均检查可观测性补量测或调整PMU位置迭代发散初值差、参考角不一致用平启动统一参考坐标系估计值偏离真值坏数据污染标准化残差检验剔除后重估角度整体偏移参考坐标系不一致比较前做角度平移J值来回跳动H矩阵拼接顺序与z不对应用结构体统一管理量测定义这张表的场景我基本都踩过尤其是“J值来回跳动”那一条当时查了好几天最后发现是H矩阵有一列索引错位导致梯度方向完全不对。5. 扩展方向与实操体会5.1 往深走的方向这套Demo做完之后还能往几个方向走。第一个方向是全PMU线性估计如果系统所有母线都装了PMU量测方程完全线性WLS就不需要迭代一次矩阵求逆就出结果实时性大幅提升。第二个方向是动态状态估计把PMU的时序数据喂给卡尔曼滤波一类的动态估计器可以跟踪状态变化趋势而不仅仅是稳态值。第三个方向是加权策略的自适应WLS权重固定工程中可以根据残差的统计特性在线调整权重抑制坏数据对估计的影响。第四个方向是把PMU电流相量量测也加进去让PMU覆盖的区域形成更强的量测冗余H矩阵的稀疏结构和数值稳定性又会有新变化。每个方向都会带来新的公式和新的坑但核心的WLS框架不变所以把基础版跑通非常值得。5.2 我的实操体会我个人的体会是状态估计这个课题最大的门槛不在数学公式而在“把量测、权重、雅可比矩阵、索引这些环节串成一个不出错的程序”。我在复现过程中踩过不少坑最深的一个是H矩阵排列顺序和量测向量不一致迭代三四轮后J值来回跳。后来我把量测定义成结构体数组每个量测自己带上类型和索引彻底解决了这个隐患。如果读者想从零复现我强烈建议按这个顺序推进先抛弃PMU用纯SCADA量测把WLS跑通再把NR潮流结果作为真值确认估计误差在可接受范围最后加入PMU量测做对比。这样每一步出错都能定位到具体模块不至于一团乱麻。这个顺序本身就是我在反复调试里总结出来的经验照着走可以少走不少弯路。
