做电力系统状态估计的人早晚会遇到一个经典问题PMU相量测量单元不是已经把电压幅值和相角都测出来了吗为什么还要用加权最小二乘WLS再估一遍干脆直接读 PMU 的数据当结果不就行了这个问题我最早也琢磨了很久。后来我把 WLS 状态估计和 Newton-RaphsonNR潮流放到同一个仿真闭环里用 PMU 输出的带噪声量测去重构电压状态再和 NR 给出的潮流状态做对比才真正明白二者之间那种“互补又不互相替代”的关系。这篇内容不是单纯的理论堆砌而是把整个项目的思路、原理、Matlab 实现和踩坑记录放在一起。你会看到PMU 的相量量测是怎么进入 WLS 的测量模型的为什么 NR 潮流结果在仿真里可以充当“真值”代码里每一步在干什么以及我在调试过程中踩过的几个典型坑。无论你是刚做课程设计还是正在写状态估计相关的算法验证程序这个闭环架构都可以直接拿过去改。1. 项目核心到底在解决什么问题1.1 状态估计和潮流计算的分工很多人刚接触这个项目时会把“状态估计”和“潮流计算”混成一回事。两者目标函数看着都像“求电压幅值和相角”但本质完全不同。潮流计算是给定了网络拓扑、发电机出力、负荷功率然后求解一个确定的非线性方程组。它的前提是“我知道系统的运行方式”输入是注入功率输出是一组满足潮流方程的状态解。NR 潮流本身就是用来解这个方程组的。状态估计则是反过来的问题调度中心并不知道系统里每一个负荷和发电机到底是多大功率它只知道 SCADA 和 PMU 传回来的带噪声的量测数据。这些量测有冗余、有误差、可能有坏数据。状态估计要做的是从一堆不完全可靠的数据里按统计最优的方式推算出系统当前最可能的状态。放到这个项目里就有了一个清晰的仿真链条先用 NR 潮流计算出“真实状态”再在真实状态下叠加传感器噪声生成 PMU 量测然后用 WLS 状态估计从这些带噪声的量测里还原状态最后把还原结果与潮流真值对比。这样你就能直接看出 WLS 到底能把误差压到多少也能量化不同量测配置对估计精度的影响。1.2 PMU 改变了什么传统 SCADA 量测主要是电压幅值、有功、无功和电流。它最大的痛点是没有同步相角量测。你想知道某个节点的相角只能靠状态估计间接算出来。PMU 出现后靠卫星授时信号让不同变电站的采样严格同步从而直接给出电压相量的幅值和角度。这意味着相角从一个“隐含状态”变成了“直接量测”。但要注意PMU 的直接量测不等于绝对准确。它一样有噪声比如电压互感器和电流互感器的传变误差、A/D 转换误差、相量算法本身的误差等。而且万一某个 PMU 出现丢帧或通道漂移坏数据带来的影响比传统量测更隐蔽因为它的量测值看起来“太准了”。所以 PMU 进入 WLS 不是“替代估计”而是“给估计提供更高品质的量测”。在 WLS 框架里PMU 量测对应的是非常简单的量测方程电压幅值 h(x)Vi相角 h(x)θi。对应的雅可比矩阵元素基本就是 0 或 1计算量小收敛很快。这也是我比较喜欢把 PMU 量测单独建模的原因调试的时候你能一眼看出问题出在哪个量测通道。1.3 为什么用 NR 潮流结果作为基准现实中你永远不知道电力系统某个瞬时的“真实状态”是多少。线路参数有误差负荷在波动传感器有偏差。但在仿真里你可以人为设置一个确定性的运行条件用 NR 潮流算出满足潮流方程的唯一状态解然后把它当成“真值”。这就是这个项目能成立的前提。运行潮流时我习惯把一台发电机设成平衡节点给定电压幅值角度设为零另一台发电机设成 PV 节点给定有功出力和电压幅值其余负荷节点都设成 PQ 节点。然后 NR 迭代算出的 V 和 θ 就是后面的基准状态。之后生成的 PMU 量测等于在这个基准状态上叠加高斯噪声。最终 WLS 还原出的状态和基准状态之间的误差就是状态估计器对测量噪声的抑制效果。有一点值得强调NR 潮流得到的“真值”只是仿真意义上的真值它本身也依赖线路参数模型。如果你把线路参数设错了整个闭环里的“真值”都是错的但这一步通常被很多人忽略。2. 算法原理拆解WLS 和 NR 的定位差异2.1 WLS 状态估计的数学模型WLS 的核心是下面这个测量模型状态向量 x 里放所有需要估计的节点电压幅值和相角z 是全部量测组成的向量h(x) 是量测函数表示“如果系统状态是 x理论上量测应该读到什么”e 是量测误差向量通常假设为零均值正态分布协方差矩阵为 R。目标函数写成J(x) [z - h(x)]ᵀ W [z - h(x)]其中 W 是权重矩阵取 R 的逆矩阵。工程直观理解是某个量测的噪声方差越小说明它越可靠权重就越大噪声大的量测权重越小状态估计不会为了迁就一个烂数据而牺牲整体精度。令 J(x) 对 x 的梯度为零经过高斯牛顿迭代得到每一步的修正量Δx (Hᵀ W H)⁻¹ Hᵀ W [z - h(x)]其中 H 是量测雅可比矩阵也就是 h(x) 对 x 的偏导数矩阵。G Hᵀ W H 通常称为信息矩阵。G 是否满秩直接决定了这个系统在当前量测配置下是否可观测。实际写代码时我从来不显式计算 (Hᵀ W H)⁻¹而是用左除dx (H * W * H) \ (H * W * r);MATLAB 的左除运算在数值稳定性上比手动求逆好很多。这一点对状态估计尤其重要因为 H 矩阵如果包含功率量测往往条件数不小显式求逆容易把误差放大。2.2 PMU 量测在 h(x) 里的写法如果把状态向量排成 x [θ1, V1, θ2, V2, θ3, V3, ...]那么 PMU 量测的写法会非常直观第 i 个节点电压幅值量测 z_Vi对应方程 h(x) Vi第 i 个节点相角量测 z_θi对应方程 h(x) θi。在雅可比矩阵 H 里电压幅值量测那一行只在 ∂h/∂Vi 位置取 1其他位置取 0相角量测那一行只在 ∂h/∂θi 位置取 1其他位置取 0。如果还加入传统有功、无功注入量测对应行就要展开成交流潮流的功率方程里面会有 Gij、Bij、cos、sin 这些项。我在实际项目里最常见的做法是PMU 量测负责电压幅值和相角传统量测负责注入功率和支路功率。这样 H 矩阵中结构简单和结构复杂两类行混在一起信息矩阵 G 的填充度更饱满可观测性也更容易满足。有一点需要特别小心PMU 的相角量测初始基准不同。NR 潮流里你通常把平衡节点相角设为 0而 PMU 输出的相角是相对卫星授时产生的全局参考相角二者之间往往差一个常数。比较估计值和真值之前必须统一参考基准否则角度误差会出现一个整体偏移看起来像系统误差其实只是参考没对齐。2.3 NR 潮流在项目里到底是什么角色NR 潮流不是状态估计它是这个项目里的“真值生成器”。NR 求解的是潮流方程组P_i V_i Σ V_j (G_ij cosθ_ij B_ij sinθ_ij) Q_i V_i Σ V_j (G_ij sinθ_ij - B_ij cosθ_ij)给定部分节点的 P、Q 和电压幅值约束迭代求解所有节点的 V 和 θ。它的雅可比矩阵是潮流方程对状态变量的偏导和状态估计里的量测雅可比 H 是完全不同的两个矩阵。一个最容易混淆的地方是NR 潮流雅可比矩阵的维数由节点类型决定PQ 节点占多数平衡节点和 PV 节点数量决定方程数目而状态估计的 H 矩阵维度由量测数量决定量测越多H 行数越多G 矩阵越容易满秩。初学者如果直接把 NR 潮流程序改个目标函数就当状态估计用很容易在可观测性和权重处理上栽跟头。我在这个项目里宁可多花半小时把“生成真值”和“状态估计”两条链路分开写也不愿意把它们耦合在一起。分开之后后续想换系统参数、加量测、加噪声场景都只需要改一个入口不用动主循环。3. Matlab 实操从建模型到出结果3.1 仿真系统与量测生成我这次用的是一个小型四节点环形网络。为什么不直接上一个大型标准测试系统因为四节点模型足够简单每个节点的注入功率、每条支路阻抗都能手写在纸面上调试雅可比矩阵时可以手工推演两行定位问题很快。先构造 Ybus。节点 1 设平衡节点节点 2 设 PV 发电机节点 3、4 设 PQ 负荷。示例 Ybus 矩阵可以是Ybus [ ... 5.0-20.0j -3.010.0j -2.010.0j 0; -3.010.0j 6.0-25.0j -3.010.0j 0; -2.010.0j -3.010.0j 7.0-30.0j -2.010.0j; 0 0 -2.010.0j 3.0-15.0j ];所有量都换算到标幺值这样电压幅值基本在 1.0 p.u. 附近数值上比较稳定。跑一次 NR 潮流之后得到基准状态 V_true 和 theta_true。然后生成 PMU 量测最简单的方式就是对真值叠加独立高斯噪声sigma_V 0.001; % 电压幅值量测标准差单位 p.u. sigma_a 0.001; % 相角量测标准差单位 rad zV V_true sigma_V * randn(nBus, 1); za theta_true sigma_a * randn(nBus, 1);这里的 sigma_V 和 sigma_a 要与 PMU 的实际相量测量误差水平大致对应。如果你希望对比不同 PMU 精度对估计结果的影响就把 sigma 按多组值扫描分别跑 WLS 并统计误差。这也是后面做蒙特卡洛分析的基础。我更建议把量测生成放在一个独立函数里返回 z 和 R 矩阵而不是在主脚本里硬编码。这样后续增加传统功率量测、改变 PMU 布点都只需要改函数内部逻辑。3.2 WLS 迭代流程与代码骨架WLS 的主循环其实非常短。关键在于 measurement_fun 函数里要正确返回量测预测值 h 和雅可比矩阵 H。我习惯把这段迭代写成x [theta0; V0]; % 初值平启动 for iter 1:30 [h, H] measurement_fun(x, Ybus, pmu_list, trad_list); r z - h; % 左除代替求逆 dx (H * (R \ H)) \ (H * (R \ r)); x x dx; if max(abs(dx)) 1e-6 break; end end theta_est x(1:nBus); V_est x(nBus1:2*nBus);注意 R 矩阵是量测误差协方差矩阵W R⁻¹。代码里写 R \ H 和 R \ r本质上是左除比 inv(R) * H 更稳。measurement_fun 的实现如果只考虑 PMU 电压幅值量测和相角量测会非常简单。对电压幅值量测那一行if strcmp(type, V) h(idx) V(i); H(idx, nBus i) 1; end对相角量测那一行if strcmp(type, theta) h(idx) theta(i); H(idx, i) 1; end如果你把传统注入功率量测也放进来情况会复杂很多。以有功注入为例量测函数是h_Pi V_i * sum( V_j * (G_ij * cos(theta_i - theta_j) B_ij * sin(theta_i - theta_j)) )偏导项要分别对 θ_i、θ_j、V_i、V_j 展开。写这种表达式时最容易出错的是符号G 是电导B 是电纳P 对应 G 在前、B 在后Q 对应 B 在前、G 在后取负号的位置也要格外小心。我在调试雅可比矩阵时有个笨但有效的办法用数值差分验证解析雅可比。构造一个小扰动对 h(xδ) 和 h(x-δ) 做差分再和解析 H 比较。只要两者最大误差在 1e-6 量级就说明偏导公式和索引位置基本没问题。这个习惯帮我避免了好几次“看着程序能跑、结果误差却始终降不下去”的诡异故障。3.3 结果对比误差指标怎么算估计完成之后自然要对 V_est、theta_est 和 V_true、theta_true 做对比。最直接的指标是绝对误差和均方根误差V_err abs(V_est - V_true); theta_err abs(theta_est - theta_true); V_rmse sqrt(mean((V_est - V_true).^2)); theta_rmse sqrt(mean((theta_est - theta_true).^2));如果只跑一次噪声的随机性会掩盖规律。所以我一般会叠加 500 次蒙特卡洛仿真每次重新生成 PMU 量测再跑 WLS最后统计误差的标准差和均值。单次结果只能用来验证代码有没有 bug统计结果才能说明算法性能。下面是一个示意性的结果表真实数值取决于具体系统参数和噪声种子节点NR 真值 VWLS 估计 V幅值绝对误差11.06001.05980.000221.04001.04010.000130.98210.98170.000440.96040.96110.0007相角误差同理。如果你发现误差明显大于量测噪声标准差先检查是不是参考相角没对齐如果误差比噪声标准差还小很多也要警惕是不是权重设置过于集中导致个别量测被过度信任。我在做误差对比时还习惯把“无 PMU 纯 SCADA 量测”和“有 PMU 混合量测”两种情况放在一起比较。这样能直观看出 PMU 对角度估计精度的提升效果而这个效果恰恰是论文和报告中很关注的结论点。4. 我踩过的坑与排查思路4.1 初值乱设导致 WLS 迭代发散WLS 本质上是牛顿类迭代对初值敏感。我刚开始试过把电压幅值初值设成全 0结果雅可比矩阵里很多偏导项直接变成 0迭代第一步就出现 NaN。后来改成平启动所有节点电压幅值取 1.0角度取 0问题立刻消失。需要注意平启动不是万能钥匙。重负荷系统、串联补偿线路等情况下平启动也可能不收敛。这时我会先用 NR 潮流算一遍把潮流解作为 WLS 初值或者直接进行分段电压启动。实际工程里状态估计前通常会有一份基于历史数据的“上一次状态”用它做初值收敛速度比平启动快得多。4.2 权重矩阵 R 设不好导致估计结果漂移权重矩阵的设置直接决定 WLS 倾向于信任哪些量测。如果误差标准差用错了单位会把量测权重放大到不合理程度。举个例子相角量测噪声如果写成 0.001 度而量测值本身用弧度那么这个噪声对应的弧度其实是约 1.745e-5 rad。相比之下如果你把另一个通道的标准差写成 0.001 rad两个量测在权重上的差距会非常悬殊最终估计结果会严重偏向那个“被误认为更准”的相角量测。所以我给自己定了一条规矩所有量测、所有误差标准差统一在同一个单位体系下计算。电压幅量测用 p.u.相角量测用 rad电阻电抗用 p.u.这样 R 矩阵的每个元素物理意义才一致。R 矩阵的构造方式通常是R diag([sigma_V1^2; sigma_V2^2; sigma_a1^2; sigma_a2^2; ...]);如果你有不同精度等级的 PMU 或传统量测就把对应的方差填到对应位置W 自然会给高精度量测更大话语权。4.3 可观测性检查与信息矩阵奇异量测太少时G Hᵀ W H 会奇异MATLAB 左除会给出警告估计值变得异常大或直接 NaN。可观测性是状态估计的前提不是“跑不出来再调”而是应该在写代码前就检查量测配置。我的检查方法很简单G H * (R \ H); fprintf(rank(G) %d, size(G) %d\n, rank(G), length(G));如果 rank 小于状态个数说明系统不可观测。这时候要么增加 PMU 布点要么补充传统功率量测和伪量测。伪量测在配电网里特别常用比如把某个节点的负荷预测当成一个有较大方差的量测也能帮助满足可观测性。在一个小系统里我遇到最典型的情况是4 号节点既没有 PMU也没有注入功率量测完全靠线路另一端的信息推测结果 G 矩阵的秩差 1。补上一个 4 号节点的电压幅值量测或注入量测后秩立即恢复估计结果也恢复正常。这个现象说明这类项目改可观测性比改算法本身更有效。4.4 常见问题速查表我把调试过程中遇到的典型问题整理成一张速查表方便你对照排查。现象可能原因排查办法迭代发散或 NaN初值不合理改用平启动或潮流解做初值角度误差整体偏移参考相角基准不一致将两套角度统一平移到同一参考节点G 矩阵奇异量测不足增加量测或伪量测检查 rank(G)估计误差远大于噪声权重单位或方差设置错误统一单位检查 R 矩阵对角线结果稳定但误差不下降可观测性弱或 PMU 布点覆盖不足增加关键节点 PMU做蒙特卡洛评估雅可比矩阵数值错误偏导公式符号或索引写错用数值差分校验解析雅可比还有一条容易忽视的规则不要一开始就迷信“量测越多越好”。量测冗余可以提高精度但如果你把两个强相关量测同时放进 R 矩阵而又没有正确处理它们的协方差反而可能导致信息矩阵条件数变差。我通常的做法是先看主对角元再看量测之间的相关性相关性高时只保留一个或做等值变换。5. 后续可扩展的方向从对比走向混合估计5.1 把 PMU 和传统 SCADA 量测一起融合进 WLS这个项目做到最后最自然的扩展就是做“混合状态估计”。传统量测覆盖全但缺少同步相角PMU 量测精度高但通常装设数量有限。把两部分量测放进同一个 WLS 模型H 矩阵里既有结构化简单的相量量测行也有传统功率量测行信息矩阵 G 的秩更容易满足估计精度通常也比单独用任一类型量测更好。实现时不要复制两份量测函数而是维护一个统一的量测清单。每条量测记录包含类型、所在节点或支路、数值、方差。measurement_fun 根据清单逐条生成 h 和 H。这样加一个量测只是往清单里添一行不需要动主循环。我实际比较过的结论是在大多数中等负荷水平下只要 PMU 覆盖到关键边界节点混合估计的角度最大误差可以从传统量测的 0.01 rad 量级降到 0.001 rad 量级。但这也取决于系统结构和量测冗余度不能一概而论。5.2 加入坏数据检测与辨识PMU 数据在真实系统里偶尔会出现异常帧或者通道漂移。WLS 对坏数据非常敏感一个突然跳变的相角量测可能把附近节点的估计结果全部拉偏。常用做法是先做目标函数 J(x) 的卡方检验再用标准化残差逐点判断r_i (z_i - h_i) / sqrt(R_ii)超过阈值就认为该量测可能是坏数据剔除后重新估计。实际编码时要注意残差检验是基于“无坏数据”假设的如果坏数据太多标准化残差检验也会失真。所以现场经验是“先整体检验再逐个剔除”而不是一次性把所有超阈值量测全部删除。5.3 动态状态估计与误差偏置修正如果继续往下走可以把静态 WLS 扩展成动态状态估计用卡尔曼滤波或扩展卡尔曼滤波处理时序 PMU 数据。这时候 NR 潮流真值仍然可以用于生成仿真轨迹比较对象也从“单一时刻估计值”变成“轨迹跟踪误差”。动态模型会比静态模型复杂很多但收益也明显能利用时间相关性在量测短暂丢失时维持状态输出。如果你未来要做 PMU 布点优化或广域监测系统这条路径会更有实践价值。最后说一个我自己现在还在用的习惯每次跑完 WLS不要只盯着最终误差把迭代过程中信息矩阵 G 的条件数和最大残差打出来看一眼。条件数异常升高时即使估计结果没发散也说明当前量测配置里存在病态因素。早发现早处理比最后拿着误差曲线猜原因省事得多。这个习惯一开始只是为排查方便后来却成了我评估量测配置好坏的快速指标也算这个项目留下的最大后遗症。