企业官网建设流程全解析

电网状态估计，说白了就是调度中心拿到SCADA送上来一堆带有随机误差的遥测数据后，用数学方法把电网真实运行状态“反推”出来的过程。这个项目聚焦在两种最经典的解法上：加权最小二乘法（WLS）和快速解耦法，测试环境是MATLAB。别觉得这是教科书上的老题，真要在MATLAB里从零把它写成一个能跑、能收敛、能处理坏数据、还能对比两种方法优劣的程序，中间有不少文档里不会写的细节。这篇文章把我实际搭建和调试这套程序的经验完整复盘一遍，适合正在做电力系统课程设计、毕业设计，或者刚开始接触能量管理系统（EMS）状态估计模块的同行参考。

1. 状态估计到底在解一道什么题

1.1 量测冗余与随机误差

先把问题的本质说清楚。电网稳态运行可以用一组节点电压幅值V和相角θ来描述，全网有N个节点，待求的状态变量就是2N-1个（参考节点相角固定为0，少一个）。调度中心的SCADA系统会送来很多量测：节点有功注入、无功注入、线路有功潮流、线路无功潮流、母线电压幅值等。量测数量m通常会远大于2N-1，这正是状态估计能“估计”的前提——冗余。

为什么不能直接拿量测值当真相用？因为互感器、变送器、A/D采样、通信传输各个环节都会叠加误差，量测噪声一般近似为零均值高斯白噪声，方差记为σ²。用一个生活类比：你称体重，同一个秤称三次得到三个数，真正的体重只有一个，直接取平均值比随便取一次更靠谱。状态估计干的就是这件事，但它不是简单求平均，而是要让成千上万个量测在潮流方程约束下达到“整体最合理”。

量测方程的标准写法是z = h(x) + e。z是m维量测向量，h(x)是由节点电压算出来的潮流值，e是量测误差。如果x是真实的电压和相角，理论上z和h(x)之间的差就只包含噪声。状态估计的任务，就是找到那个让噪声“整体最小”的x̂。

1.2 为什么偏偏选了WLS和快速解耦法

WLS是这个问题的标准答案。它在高斯噪声假设下等价于最大似然估计，统计性质好，残差平方和最小，适合离线计算、事后分析和作为校验基准。但它的缺点是每次迭代都要重新计算雅可比矩阵H，再对信息矩阵G做一次因子分解，节点规模一大，单次迭代计算量就上去了。

快速解耦法则是在WLS迭代结构上“挖掉”不敏感耦合项得到的近似解法。高压电网里有功P主要跟相角θ耦合，无功Q主要跟电压幅值V耦合，把增益矩阵近似成两组与迭代过程无关的常数矩阵后，因子分解只需要做一次，每轮迭代只是两次回代，计算速度比WLS快好几倍，因此在线EMS里用得更多。

一个项目里同时写这两种方法，正好覆盖“精算”和“快算”两个定位。WLS的结果可以当标准答案，快速解耦法用来对比计算效率，顺便还能验证在特定测试系统下两种方法收敛性和精度的真实差异。这种对比做完了，你对状态估计的理解会扎实很多。

2. 加权最小二乘法的MATLAB实现细节

2.1 从目标函数到牛顿迭代公式

加权最小二乘法的目标函数是：

J(x) = (z - h(x))ᵀ W (z - h(x))

其中W = R⁻¹，R是对角量测协方差矩阵。权重W不是随手拍的，R对角元是各量测的方差σ_i²，量测精度越高，方差越小，权重越大。物理意思就是：相信高精度表计，不要让低精度表计把结果带偏。

目标函数对x求导并令其为零，整理后得到：

G Δx = Hᵀ W (z - h(x))，其中 G = Hᵀ W H

H是m×n的雅可比矩阵（n = 2N-1），G是n×n方阵，叫信息矩阵或增益矩阵。用牛顿法迭代求解，每次解这个线性方程组即可。

这里有个容易被忽略的点：量测误差的方差σ²是怎么来的？实际工程里依据的是量测装置精度等级。比如电压互感器0.2级，满量程100V，最大误差0.2V，标准差大约取0.067V；功率变送器类似。在标幺制下，电压量测的σ可以取0.004左右，功率量测的σ取0.01左右，这样权重W_i = 1/σ_i²的量级大约在10⁴，数值上不会因为权重差异太大导致病态。

2.2 组装H矩阵：列顺序对齐状态量

写代码时最繁琐的就是H矩阵。先给每类量测定义类型号，比如：1代表节点电压幅值量测，2代表节点有功注入量测，3代表节点无功注入量测，4代表支路有功潮流量测，5代表支路无功潮流量测。然后逐行填充雅可比矩阵。

节点电压幅值量测对应的行很简单，只有本节点电压V一个非零导数，值为1。支路潮流量测那一行就复杂了，会同时涉及到首末端节点的V和θ。比如线路i-j从i端看过去的有功潮流Pij = Vi²g - ViVj(gcosθij + bsinθij)，它对θi、θj、Vi、Vj分别求导，就会占四个列位置。

组装时的核心教训是：状态量x的内部顺序要固定且全程序统一。比如规定前N-1个元素是除参考节点外的相角，后N个元素是全部节点电压幅值。H矩阵的列就按这个顺序排列，参考节点相角不是状态量，求导后那一列直接删掉，不然G矩阵必然奇异。

```matlab %% WLS核心迭代循环（伪代码框架） x = [zeros(N-1,1); ones(N,1)]; % 相角初值0，电压初值1（平启动） for iter = 1:15 [hx, H] = calc_meas(x, meas); % 计算量测函数和雅可比矩阵 r = z - hx; % 残差向量 G = H' * W * H; % 增益矩阵 dx = G \ (H' * W * r); % 求解增量，直接用左除 x = x + dx; if max(abs(dx)) < 1e-5 break; end end

还有一个非常影响调试效率的细节：雅可比矩阵要建成稀疏矩阵。MATLAB里用sparse函数构建，或者先预分配满阵到循环里再转稀疏。14节点系统满阵也就27×27，无所谓；但随便扩到几百节点，满阵G运算速度会明显变慢，尤其多次因子分解时差别巨大。

2.3 量测函数的计算细节

量测函数h(x)的核心是潮流公式。节点注入有功、节点注入无功、支路两端有功无功潮流，一共就这几类。在MATLAB里写几个子函数分别算，比全塞在主循环里可读性好得多。

节点注入量的公式是：

Pi = Vi * Σ Vj (Gij cosθij + Bij sinθij)

Qi = Vi * Σ Vj (Gij sinθij - Bij cosθij)

支路潮流的公式要看方向。从i端到j端的有功潮流是：

Pij = Vi² g - Vi Vj (g cosθij + b sinθij)

无功潮流：

Qij = -Vi² (b + yc/2) - Vi Vj (g sinθij - b cosθij)

注意其中g、b是支路导纳的实部和虚部，yc/2是对地导纳的一半。变压器的变比k也要折算进去，否则算出来的Pij会跟标准测试系统的潮流结果对不上。这个“对不上”排查起来很痛苦，因为解可能收敛，但是估计值整体偏移。

我建议在状态估计之前，先做一次潮流计算，把量测函数子程序用潮流的电压和相角去算一组“伪量测”，跟潮流结果对比。如果误差在1e-8以下，说明量测函数没写错；如果对不上，抓紧修公式，不要等到状态估计跑偏了再去怀疑量测函数。

3. 快速解耦法为什么能省下一大截计算

3.1 P-Q解耦的工程基础

高压输电网里有个非常实用的物理特性：线路电抗远大于电阻，R/X一般小于0.1。这种参数结构下，潮流方程里节点有功主要与相角差相关，无功主要与电压幅值差相关。也就是说P对V的偏导很小，Q对θ的偏导也很小。

WLS在算H矩阵时把这四块耦合关系全部算进去了，信息矩阵G经常要重算。快速解耦法直接预先把弱耦合块清零，只保留P-θ和Q-V两条主通道，得到两组独立方程：

ΔP/V = B' Δθ

ΔQ/V = B'' ΔV

写得更通俗一点：有功不平衡量除以电压幅值，用常数矩阵B'求相角修正量；无功不平衡量除以电压幅值，用常数矩阵B''求电压修正量。两个方程交替迭代，这就是“快速解耦”的全部秘密。

3.2 B'、B''两个常数矩阵的构造

B'和B''的构造是本方法最需要严谨对待的地方。它们不是直接从H矩阵里抽出来的，而是从系统节点导纳矩阵Y的虚部里取出来的。

具体做法是：先构建节点导纳矩阵Ybus，Ybus的虚部取负号，得到B0 = -imag(Ybus)。然后删掉参考节点对应的行和列，得到B'。B''的做法类似，但为了更贴合“忽略充电电容”的简化假设，在构建Ybus时可以把所有支路对地导纳yc/2先去掉，再取虚部删参考节点行和列。

为什么B''要忽略充电电容？因为充电电容主要影响无功和电压，但高压线路充电功率在轻载时对电压幅值影响不小，如果B''太粗略，Q-V迭代会变慢甚至震荡。工程上有多种变体，最常用的是B''中包含充电电容影响，B'不包含。你可以在代码里做成开关选项，对比哪种收敛更快。

%% 构造B'和B'' % 支路数据: [首端 末端 r x 对地导纳1/2 变比k] nb = size(branch,1); Ybus = zeros(N,N); for k = 1:nb i = branch(k,1); j = branch(k,2); y = 1 / (branch(k,3) + 1j*branch(k,4)); Ybus(i,i) = Ybus(i,i) + y; Ybus(j,j) = Ybus(j,j) + y; Ybus(i,j) = Ybus(i,j) - y; Ybus(j,i) = Ybus(j,i) - y; Ybus(i,i) = Ybus(i,i) + 1j*branch(k,5); Ybus(j,j) = Ybus(j,j) + 1j*branch(k,5); end Bmat = -imag(Ybus); id = setdiff(1:N, ref); % 去掉参考节点 Bp = Bmat(id, id); % B'' 忽略充电电容版本 Ybus2 = Ybus; for k = 1:nb % 把对地导纳去掉 i = branch(k,1); j = branch(k,2); Ybus2(i,i) = Ybus2(i,i) - 1j*branch(k,5); Ybus2(j,j) = Ybus2(j,j) - 1j*branch(k,5); end Bpp = -imag(Ybus2); Bpp = Bpp(id, id);

注意Bp和Bpp都是常数矩阵，整个迭代过程只需要构建一次，也只需要做一次因子分解。MATLAB里直接用左除也行，但更高效的做法是先用chol做一次Cholesky分解，之后每轮迭代只回代两次。测试时14节点系统看不出来太大差距，节点数到几百时这个优化很值。

3.3 快速解耦法的迭代主循环

迭代逻辑是：先算当前状态下的有功不平衡量ΔP，解B'方程得到Δθ，更新相角；再用更新后的相角算无功不平衡量ΔQ，解B''方程得到ΔV，更新电压幅值。如此交替往复。

%% 快速解耦法主循环 theta = zeros(N-1,1); % 参考节点以外相角 V = ones(N,1); % 全部节点电压幅值 [Lp, ~] = chol(Bp); % 预分解，Bp是对称正定 [Lpp, ~] = chol(Bpp); for iter = 1:15 x = [theta; V]; % 计算节点注入有功Pcal和支路潮流 Pcal = calc_node_inject(x, meas, 'P'); dP = (meas.P - Pcal) ./ V_active; % 去掉参考节点行 dtheta = Bp \ dP; % 实际可改用 Lp theta = theta + dtheta; x = [theta; V]; Qcal = calc_node_inject(x, meas, 'Q'); dQ = (meas.Q - Qcal) ./ V_active; dV = Bpp \ dQ; V = V + dV; if max([abs(dtheta); abs(dV)]) < 1e-5 break; end end

有个细节：如果Bp和Bpp不是正定矩阵（原因通常是某些节点孤立或量测配置不足），chol会报错。这时就用普通左除，虽然慢一点，至少能出结果。工程实现中更稳妥的办法是调用MATLAB的lu分解或直接用符号分解工具去处理稀疏矩阵。

3.4 两种方法在测试系统上的对比结果

我在标准14节点测试系统上跑了同样的量测数据，统计结果大致如下：

快速解耦法虽然迭代次数多，但每次迭代代价极低，总耗时反而少一半以上。这里要注意：如果线路R/X过高（比如配电网），P-Q解耦的假设就不成立了，快速解耦法可能迭代十几次还不收敛。我做测试时用的输电网数据没问题，但换成低电压网络就栽过跟头。

4. 从平启动到收敛：完整跑通一次状态估计

4.1 测试数据准备与量测生成

环境用MATLAB即可，版本没有太多要求，R2016以后都行。数据规模建议先拿标准14节点测试系统练手，节点数不大不小，支路潮流和电压分布有变化，能清楚看出两种方法的差异。

数据准备分三步。第一步，把系统母线数据和支路数据整理成统一格式，自行构建导纳矩阵Ybus。第二步，用潮流程序算出一组节点电压和相角，作为“真值”。第三步，在每个量测的真值上叠加高斯噪声，生成模拟遥测数据。电压幅值量测的σ取0.004标幺，功率量测的σ取0.01标幺，这样冗余量测下的估计结果应该能恢复到真值附近。

量测配置要有冗余。我给每条支路两端各配一组有功无功潮流量测，在大部分非参考节点配注入功率量测，在所有节点配电压幅值量测。算下来量测数m要多于状态量n，最好做到1.5到2倍以上，否则估计结果会被少量量测主导，坏数据辨识也难做。

4.2 主程序框架和执行流程

整个程序建议拆成四个文件：数据读取与导纳矩阵构建、量测函数与雅可比计算、WLS循环、快速解耦循环。主脚本负责调用和展示结果。这样结构清晰，后期加坏数据辨识也容易。

%% 主脚本示例 clear; clc; run('data14.m'); % 载入母线、支路、量测数据 [Ybus] = build_ybus(bus, branch); N = size(Ybus,1); ref = 1; % 参考节点编号 meas = build_meas_cases(N, Ybus, sigma); % 生成模拟量测 % 运行WLS [x_wls, hist_wls] = run_wls(Ybus, meas, N, ref); % 运行快速解耦法 [x_fd, hist_fd] = run_fd(Ybus, meas, N, ref); % 画收敛曲线 figure; plot(hist_wls.dx_max, '.-', 'LineWidth', 1.2); hold on; plot(hist_fd.dx_max, 's-', 'LineWidth', 1.2); legend('WLS', '快速解耦法'); xlabel('迭代次数'); ylabel('max|Δx| (p.u.)'); set(gca, 'YScale', 'log');

迭代记录dx_max是我强烈建议保留的调试信息。把每次迭代的最大状态量变化存下来，最后画成对数坐标曲线，一眼就能看出收敛速度。WLS是典型的三次左右陡降，快速解耦法是阶梯式下降，每一轮交替消去相角和电压的误差。

4.3 结果怎么看：收敛曲线与误差评估

判敛标准建议用max|Δx|小于1e-5标幺。为什么不建议用目标函数J的变化量判敛？因为J在接近极小值的时候变化很平缓，可能迭代没收敛J就不怎么变了，判断过早结束，误差比较大。用状态量变化量最直接。

正常收敛时，WLS一般在第2到第3次迭代达到1e-5以下，快速解耦法大概需要4到7次。如果WLS超过10次才收敛或者根本不收敛，不用怀疑公式，先回头查数据：参考节点有没有删干净、量测残差量级是否异常、导纳矩阵是否对称。大概率就这三个地方出错。

误差评估可以把估计值跟潮流真值比一下平均绝对误差。正常情况下电压幅值误差在0.001到0.005标幺，相角误差在0.0001到0.001弧度。如果你看到某个节点的功率估计误差明显偏大，先检查该节点的量测配置和权重设置，不要急着调收敛阈值。

5. 坏数据识别：别让一个坏遥测带偏整张网

5.1 残差检测的原理

状态估计做完之后，残差向量r = z - h(x̂)包含丰富的信息。在量测噪声正常、模型正确的前提下，残差服从已知的统计分布。假如某个量测因为变送器故障或者通信误码混入了严重偏差，状态估计会试图“照顾”这个坏数据，导致周围节点的估计结果都被拉偏，残差结构出现明显异常。

检测手段最常用的是标准化残差。把每个残差除以它的标准差，得到：

r_i^std = |r_i| / σ_i

如果某个量测是坏数据，它的标准化残差会远大于正常水平。阈值一般取3到5，对应误检概率大约0.27%到0.00006%。实际调试中我习惯先用5筛选，再用3复查。阈值取得太小容易把好数据误删，取得太大则坏数据漏检。

严格意义上的标准化残差还要考虑量测函数和增益矩阵的影响，公式稍微复杂一点。但对于课程设计和工程初版，直接用残差除以σ_i已经能解决绝大多数坏数据问题，等发现问题再升级算法也来得及。

5.2 坏数据检测与剔除的迭代流程

坏数据处理要迭代进行，不能一次把所有超限的量测全删掉。因为删除量测会改变估计结果和残差分布，原来超限的量测可能恢复正常，原来正常的量测反而可能超限。稳妥的做法是每次删一个标准化残差最大的量测，删完重新跑状态估计，再重复检测，直到没有超限量测为止。

%% 坏数据辨识循环（伪代码） while true [x, ~] = run_wls(Ybus, meas, N, ref); % 重跑WLS [hx, ~] = calc_meas(x, meas); r = meas.z - hx; r_std = abs(r) ./ meas.sigma; % 标准化残差 [max_r, idx] = max(r_std); if max_r < 3 break; end fprintf('剔除量测 %d，标准化残差 %.2f\n', idx, max_r); meas.z(idx) = []; meas.sigma(idx) = []; end

剔除坏数据之后，最好对比一下剔除前后的估计结果。如果整网电压、相角变化很小，说明坏数据的影响被冗余量测压制住了，这是高冗余带来的好处；如果变化很大，说明量测冗余不足，坏数据已经污染了估计结果，此时要考虑补充伪量测或者增加量测点。

实际做测试时我给某个支路功率量测叠加了20倍σ的偏差，WLS初版结果里该支路两端节点的电压估计都被拉偏了0.02标幺左右，标准化残差检测出来后，剔除一次就恢复正常。这个实验强烈建议你也做一遍，对理解“坏数据为什么可怕”非常直观。

6. 调试过程中我踩过的坑与经验

6.1 不收敛，先别怀疑公式

我在第一次跑WLS时，程序怎么都不收敛，迭代五六次后误差反而增大。排查了半天，最后发现是平启动时相角用了角度制，量测函数里却用了弧度制，等于自己给自己制造初值误差。这个错误太隐蔽了，因为程序不报错，数值也不至于NaN，就是收敛得很慢。

另一个常见坑是参考节点相角忘记固定。状态向量里包含参考节点相角，H矩阵里也没删那一列，G矩阵必然奇异，左除出来一堆Inf或者巨大的数。排查方法很简单：在迭代前打印一次量测残差r = z - h(x0)，如果量级在1e-1以上，说明量测函数或单位有错；如果r很小，说明问题在迭代环节。

6.2 矩阵奇异与可观测性问题

可观测性的数学表现是G矩阵奇异或病态。实际调试中可以做个特征值分析，把G的特征值排序看看有几个接近0的小特征值。N=14的系统，如果存在多个远小于1的特征值，基本可以断定存在不可观测岛。

补救方法有三条：一是增加该区域的注入功率量测或支路潮流量测；二是加伪量测，比如把母线电压初值当作一个低权重的量测塞进去；三是删掉不可观测岛，把问题缩小到可观测部分。课程设计里最常用的是加伪量测，权重给很小，比如W = 1e-2，这样既不干扰正常量测，又能保证G可逆。

6.3 快速解耦法在重负荷和低压网下的失效经验

当支路R/X比值过高（比如配电网），快速解耦法的P-Q解耦假设就不成立了，迭代次数急剧增加，甚至发散。我在测试低电压网络时踩过这个坑，现象是前两次迭代还挺正常，第三次开始dTheta和dV来回震荡，怎么都压不下去。

解决思路有两个。一是把B'矩阵从“只取电抗倒数”改成“用实际导纳矩阵虚部”，提高近似精度；二是干脆回到WLS，别硬撑着用快速解耦法。对于输电网标准测试系统，快速解耦法表现稳定，我建议在项目报告里注明：本方法的适用边界是高压输电网，R/X较低的场景需要谨慎。

6.4 调试过程的记录习惯

最后分享一个实用的调试习惯：每次迭代都把残差范数和最大状态量变化记录下来，最后画成图。如果观察到残差先降后升，考虑步长过大；如果残差单调下降但很慢，考虑量测权重设置不合理；如果残差完全没变化，检查量测数据有没有加载成功。状态估计的调试就是这么一层一层查出来的，记录充分的话，问题定位会快很多。

这个项目做完之后，我最大的感受是：状态估计的难点不在公式推导，而在工程实现的一致性。单位、列顺序、参考节点、权重设置，任何一个细节没对齐，结果都能“优雅地错”。如果你正打算在MATLAB里复现这套程序，建议从14节点系统开始，先跑WLS验证正确性，再对比快速解耦法，最后加上坏数据辨识，循序渐进，每一步都确认无误后再往下走。