企业官网建设流程全解析

刚接触状态估计那会儿，我有个很深的困惑：既然电力系统的潮流方程能算出母线电压，为什么调度中心还要单独搞一套状态估计？后来在工程现场看多了就明白了，潮流计算输入的是确定的发电出力和负荷功率，算出来的是一个“理想断面”；而调度终端拿到的SCADA数据、PMU数据里面全是测量误差、坏数据和时标不一致，你要是直接把量测值拿去用，很可能被一条错误数据带到沟里。这个项目正好把两个关键问题串在了一起——一边是用加权最小二乘（WLS）结合PMU（相量测量单元）对电力系统的电压幅值和相角做状态估计，另一边是用经典的Newton-Raphson潮流计算结果作为参照基准，在Matlab里把两者逐母线比较，看同步相量测量到底能给状态估计带来多大的精度提升。

• 第1部分 H2 ->
• 第2部分 H2 ->

如果你正在学电力系统状态估计，或者想弄明白PMU量测到底怎么进WLS估计器，这篇博文的思路可以直接拿过去套。

1. 从潮流计算到状态估计：PMU量测带来的问题差异

1.1 状态估计到底在解决什么问题

调度中心常说的“状态”，指的是系统中每条母线的电压幅值 V 和相角 θ。潮流计算也能给出一组 V 和 θ，但潮流计算这件事有个前提：所有节点的注入功率必须准确已知，然后你去解一组确定性的非线性方程组。现场不是这样的，变压器变比会有偏差，互感器存在比差角差，遥测通道偶尔还会冒出个坏数。换句话说，你拿到的量测天生就带着噪声，你没法把它当作精确输入来用。

状态估计的定位恰好在这里：它把量测当成随机变量，在统计意义下寻找一个最有可能产生当前量测的系统状态。你一上来不需要“精确的注入功率”，你只需要一大堆带误差的量测，然后按它们的可靠程度加权，最后拟合出一组电压幅值和相角，让模型输出和量测之间的加权残差最小。

刚入门的时候容易把状态估计和潮流计算搞混，后来我用一句话区分它们：潮流计算是在输入精确的前提下求方程组的解；状态估计是在输入带误差的前提下求最可信的拟合解。这个区别直接决定了两个问题要用完全不同的数学工具。

1.2 PMU量测带来的最大变化

传统SCADA量测主要有三类：母线注入功率、线路潮流、电压幅值。唯独电压相角很难测，原因很直白——整个电网没有一个统一的时标基准，各变电站的遥测值之间无法比较相角关系。在WLS状态估计里，这就迫使你只能通过功率量测反推相角，量测冗余度不够的时候，估计精度和收敛性都会受到影响。

PMU（相量测量单元）改变了这个局面。它依靠卫星授时，对电压和电流波形做同步采样，可以直接输出带时标的电压相量，也就是同时给出幅值和相角。这意味着在WLS的测量方程里，原来最难观测的 θ 现在可以直接测到，雅可比矩阵里对应行的形式极其简单，权重也可以给得比较大。

这个项目里我们就是把这些优点全部用上：把PMU提供的电压幅值、电压相角量测混入WLS估计器，再和Newton-Raphson潮流解做对比。PMU测得的相角精度通常在0.1°级别，比通过功率量测反推出来的相角稳定得多，所以理论上WLS+PMU的估计结果应该非常接近潮流真值。

1.3 本文的整体路线

整个项目的实现路线我按下面四步走：

1. 在Matlab里搭一个测试系统，我先用IEEE 14节点系统做试验台，当然你也可以换成IEEE 30节点或者更简单的三节点系统。
2. 用Newton-Raphson潮流计算得到一组“参考状态”，把它当作系统在该运行方式下的“真值”。
3. 在参考状态的基础上生成带噪声的量测数据，模拟PMU的电压幅值、相角量测以及SCADA的功率量测。
4. 用WLS迭代估计系统状态，逐母线比较估计值和Newton-Raphson参考值之间的误差。

下面每一章我都会把原理和代码对应起来讲，避免出现“公式看不懂、代码跑不通”两头空的状况。

2. WLS状态估计的数学基础：从测量方程到增益矩阵

2.1 测量模型和加权目标函数

电力系统状态估计的测量模型可以统一写成：

z = h(x) + e

这里 z 是m维量测向量，包括PMU电压幅值、PMU电压相角、注入有功无功、线路潮流等；x 是n维状态向量，通常取所有非参考母线的相角和所有母线的电压幅值；h(x) 是量测函数，由潮流方程决定；e 是量测误差向量，工程上一般假设它服从零均值高斯分布。

既然误差是随机的，我们没有办法逐点消除它，只能追求一个整体最优。WLS的做法是让加权残差平方和最小，目标函数写成：

J(x) = [z - h(x)]ᵀ R⁻¹ [z - h(x)]

其中 R 是量测误差协方差矩阵，典型形式是一个对角矩阵，对角元素是各量测方差σᵢ²。为什么用方差倒数做权重？原因很直观：方差大的量测可信度低，在目标函数里占的权重就小；方差小的量测精度高，理应获得更大的话语权。这个逻辑和日常加权平均完全一致，只不过这里是矩阵形式，可以把不同量测类型的相关性也考虑进去。

2.2 从一阶最优条件到迭代修正

要求 J(x) 最小，令其一阶导数为零：

∂J/∂x = -2 Hᵀ R⁻¹ [z - h(x)] = 0

其中 H = ∂h/∂x 是量测雅可比矩阵。h(x) 是x的非线性函数，没法直接解，需要用Newton类迭代来处理。把 h(x) 在当前状态 xₖ 附近做一阶泰勒展开，代入一阶最优条件后可以得到：

G Δx = Hᵀ R⁻¹ [z - h(xₖ)]

其中 G = Hᵀ R⁻¹ H 称为增益矩阵。每次迭代解这个线性方程组，得到状态修正量 Δx：

Δx = G⁻¹ Hᵀ R⁻¹ [z - h(xₖ)]

然后把 xₖ₊₁ = xₖ + Δx，继续算新的H和h，直到 Δx 的范数小于收敛阈值。

这个形式你可能看着眼熟，它就是Gauss-Newton法在加权最小二乘问题里的标准形态。增益矩阵G的维数是状态数×状态数，在每个迭代步只需要做一次矩阵分解，不需要像潮流那样反复求高阶导，所以计算效率是可以接受的。

2.3 PMU量测如何进入雅可比矩阵

我分别列出几类常见量测对应的H矩阵行向量，这样你在Matlab里搭矩阵时会很清晰。先假设状态向量排序为 x = [θ₁, θ₂, …, θ_N, V₁, V₂, …, V_N]。

第一类，母线i的电压幅值量测，测量方程是：

Vᵢᵐ = Vᵢ

对应雅可比行是：这一行除了第 N+i 个元素为1，其余全为0。

第二类，母线i的电压相角量测（PMU直接给出），测量方程是：

θᵢᵐ = θᵢ

对应雅可比行是：这一行除了第 i 个元素为1，其余全为0。

第三类，母线i注入有功量测，测量函数是：

Pᵢ = Vᵢ ∑ⱼ Vⱼ (Gᵢⱼ cosθᵢⱼ + Bᵢⱼ sinθᵢⱼ)

其中 θᵢⱼ = θᵢ - θⱼ。这一行就比较复杂，需要对 θ 和 V 分别求偏导，典型项包括：

∂Pᵢ/∂θᵢ = ∑ⱼ≠ᵢ VᵢVⱼ(-Gᵢⱼ sinθᵢⱼ + Bᵢⱼ cosθᵢⱼ) ∂Pᵢ/∂θⱼ = VᵢVⱼ(Gᵢⱼ sinθᵢⱼ - Bᵢⱼ cosθᵢⱼ) ∂Pᵢ/∂Vᵢ = 2VᵢGᵢᵢ + ∑ⱼ≠ᵢ Vⱼ(Gᵢⱼ cosθᵢⱼ + Bᵢⱼ sinθᵢⱼ) ∂Pᵢ/∂Vⱼ = Vᵢ(Gᵢⱼ cosθᵢⱼ + Bᵢⱼ sinθᵢⱼ)

这里 Gᵢⱼ、Bᵢⱼ 是节点导纳矩阵Y的实部和虚部。功率量测的雅可比公式和Newton-Raphson潮流里的子矩阵非常像，但注意在WLS里所有量测行都要放进同一个H里，行数等于量测数量，列数等于状态数量。

第四类，线路 ij 的有功潮流量测，测量函数要考虑线路对地电容和变压器变比，公式比注入功率略长，但求导思路一样。

放PMU相角量测的最大好处，从H矩阵就能看出来：它的行是高度稀疏的单位向量，不像功率量测那样每个元素都是三角函数堆出来的。也就是说，PMU不仅提供了直接可观性，还让增益矩阵G的稀疏结构和数值条件都变得更友好。

3. 试验台搭建：用Newton-Raphson潮流给状态估计器提供参考基准

3.1 为什么要先算一次潮流

状态估计做得好不好，得有一个参照物。PMU和SCADA量测都是带噪声的，你不能拿带噪声的量测直接去验证自己的结果。工程上最常用的做法是：先在Matlab里搭建一个标准测试系统，用Newton-Raphson法把该运行方式下的潮流算出来，把潮流解当作“真值”，然后在真值基础上叠加模拟噪声生成量测，这样你就可以精确地比较WLS估计值和真值之间的差异。

这个思路本质上是一种蒙特卡洛仿真思路，比拿现场数据验证要干净得多。现场数据你永远不知道真正的运行状态是什么，仿真算例却可以做到“心里有数”。这也是为什么论文和项目里普遍采用IEEE标准测试系统的原因。

3.2 组装测试系统数据和量测数据

我在Matlab里用的是IEEE 14节点系统的标准参数，数据组织采用Matpower风格，用矩阵存母线数据、支路数据和发电机数据。母线数据矩阵 bus 每一行对应一条母线，列分别表示母线编号、类型、有功负荷、无功负荷、电导、电纳、电压幅值初值、电压相角初值等。支路数据矩阵 branch 每一行对应一条线路，列分别表示首端母线、末端母线、电阻R、电抗X、对地电纳B、变比等。

如果你手头没有现成的IEEE 14节点数据，直接在Matlab命令行输入下面这行也行，Matpower工具箱自带这个案例：

mpc = loadcase('case14');

不过要注意，这需要安装Matpower工具箱。要是你只是想跑一个最小可用版本，我建议先搭一个三节点系统，手写导纳矩阵也就几十行代码，调起问题来方便得多。本文后面的代码逻辑不依赖具体节点数，你换成IEEE 14或者IEEE 30都只需要换数据矩阵，核心算法不用动。

量测数据的生成分两步。第一步，用Newton-Raphson潮流算出参考状态 x_true；第二步，在参考状态上叠加高斯噪声。比如要模拟一条电压幅值量测，就写成：

zV(i) = V_true(i) + sigmaV * randn;

相角量测类似，但注意要把角度换算成弧度。功率量测的噪声标准差通常按满量程的百分数取，我用的是0.02 pu，电压幅值用0.005 pu，相角用0.01 rad。

3.3 Newton-Raphson潮流的Matlab实现要点

Newton-Raphson潮流本身也是一个迭代求解过程。它的核心思路是把PQ节点的有功无功失配量和PV节点的有功失配量，通过雅可比矩阵映射到电压幅值和相角的修正量，反复迭代直到失配量小于阈值。

核心循环大致是：

for iter = 1:maxIter % 计算失配量 [dP, dQ] = mismatch(bus, branch, V, theta); % 组装雅可比矩阵J J = jacobian_nr(bus, branch, V, theta); % 求解修正方程 dx = J \ [dP; dQ]; % 更新状态 theta = theta + dx(1:nPQ); V = V + dx(nPQ+1:end); % 判断收敛 if norm([dP; dQ], inf) < tol break; end end

这里最值得注意的是雅可比矩阵的组装，它的子矩阵形式和WLS里的H矩阵很像，但维度和符号约定略有不同。很多人在Matlab里写牛顿法，公式明明抄对了，结果一直不收敛，往往是平衡节点和PV节点的方程排布搞错了。我的建议是在开始写代码前，先把节点类型标清楚：平衡节点不参与失配量计算，PV节点没有无功失配方程。

把Newton-Raphson潮流结果存下来，作为后面比较的基准。在理想状况下，WLS使用无限精度的量测和正确的权重，估计结果会和潮流解完全重合；加入测量噪声后，两者会有一个合理的偏差，这个偏差的大小正是我们评估PMU量测价值的关键指标。

4. Matlab实现WLS+PMU状态估计的完整过程

4.1 初始化：母线数据、支路数据、量测数据

状态估计代码的第一步是读入系统数据并组织成程序所需的结构。我在代码里用结构体存数据，这样可读性更好。基本流程是：

nb = size(bus, 1); % 母线数 nl = size(branch, 1); % 支路数 Y = ybus(bus, branch); % 节点导纳矩阵 G = real(Y); B = imag(Y);

量测数据我统一存在一个向量z里，同时用一个量测类型标签数组type记录每条量测的种类：1代表电压幅值，2代表电压相角，3代表节点注入有功，4代表节点注入无功。这样在计算H矩阵和h(x)时可以循环遍历每一条量测，根据类型调用不同的计算分支，代码结构非常清晰，也方便你后续加入线路潮流量测。

对于这个小项目，我建议把所有量测组织成一个结构体或者cell数组，不要散落成多个变量。后面调试的时候，你只需要看这一个结构体就能搞清楚整套量测配置。

4.2 构建量测雅可比矩阵H和权重矩阵R

权重矩阵R我用对角阵Rdiag表示，这样内存开销小，求逆也方便，直接用逐元素除法代替完整矩阵求逆。Rdiag里第k个元素就是第k条量测的方差σₖ²。构建方式如下：

Rdiag = ones(m, 1); for k = 1:m switch type(k) case 1 % 电压幅值量测 Rdiag(k) = sigmaV^2; case 2 % 电压相角量测 Rdiag(k) = sigmaTheta^2; case 3 % 注入有功量测 Rdiag(k) = sigmaP^2; case 4 % 注入无功量测 Rdiag(k) = sigmaQ^2; end end invR = 1 ./ Rdiag;

H矩阵的构建是WLS实现的核心，也是出错率最高的地方。我的做法是先初始化一个全零的m行2N列矩阵，然后逐条量测填充对应行。以三相角量测为例：

% 第k条量测是母线i的电压相角 H(k, i) = 1; % 如果母线i是平衡节点，并且我们把平衡节点相角固定在0，则这一行需要删除

注意这里有个特别容易踩的坑：如果状态向量里保留平衡节点的相角θ₁，但量测模型里没有相角参考，增益矩阵G就会奇异。处理办法有几种，最简单的是直接把平衡节点相角从状态向量中剔除，也就是把θ₁从x里拿掉，在每次迭代收敛后再补回0。另一种办法是给G矩阵中对应平衡节点相角的位置加一个很大的权重惩罚项。我推荐第一种，逻辑最干净。

4.3 WLS迭代主循环

在主循环里，每轮迭代做四件事：计算当前状态下的量测函数h(x)，计算雅可比矩阵H，组装增益矩阵G，求解状态修正量Δx。代码主线如下：

x = zeros(2*nb - 1, 1); % 不含平衡节点相角 x(1:nb-1) = 0; % 非平衡节点相角初值 x(nb:end) = 1; % 所有母线电压幅值初值 maxIter = 20; tol = 1e-8; for iter = 1:maxIter h = calc_h(x, bus, branch, type, nb); H = calc_H(x, bus, branch, type, nb, G, B); dz = z - h; Gmat = H' * (invR .* H); % 注意这是逐行加权再求和 dx = Gmat \ (H' * (invR .* dz)); x = x + dx; if norm(dx, inf) < tol fprintf('WLS收敛于第%d次迭代\n', iter); break; end end

这里面有一个Matlab性能小技巧：不要直接构造完整的R⁻¹矩阵再和H做矩阵乘法，而是用 invR .* H 这种按行缩放的方式，在大规模系统下能省不少内存和时间。

h(x)的计算函数里，需要对每条量测判断类型，再代入当前x计算对应量测函数值。V_i和V_j从x的对应位置取，θ_i和θ_j同理，注意关节角度差θᵢⱼ的符号，符号反了结果完全不对。

4.4 Matlab运行效果与常见坑位

我在Matlab R2023b上跑这个程序，IEEE 14节点系统，量测数量配置为14个电压幅值、5个PMU相角、14个注入有功、14个注入无功，迭代基本在4到6次内收敛，单次迭代时间不到0.1秒，整个估计过程感觉不到任何延迟。

遇到的第一个坑是量测单位不统一。导纳矩阵是标幺值，功率量测是标幺值，这些没问题，但相角我一开始习惯用度数存，结果H矩阵里三角函数的输入和量测残差全部错位。后来统一规定：内部计算全部用弧度，只在最后输出对比表时把角度转成度数。

第二个坑是平衡节点相角的约束问题。我之前试过在每次迭代后手动把x(1)强制置0，但这样做的效果等同于是硬性改了状态更新方向，导致某些迭代步修正量异常。后来改成从状态向量里直接剔除平衡节点相角，反而简单可靠。

第三个坑在前面第2部分也提过，H矩阵的行和量测z的行顺序必须严格一致。我在调试早期犯过把z排序和type排序搞反的低级错误，结果估计出来的电压幅值整体偏移了一个母线位置，数值看起来“很合理”但实际错得离谱。建议在循环前加一行assert检查：

assert(length(z) == size(H,1), '量测向量和H矩阵行数不匹配');

5. 对比结果：电压幅值与相角的逐母线差异分析

5.1 无噪声理想量测下的估计精度

我先跑了一个“理想化”实验：量测噪声全设为0，也就是z完全等于h(x_true)。这种情况下WLS估计值应该和Newton-Raphson潮流解几乎完全一致，差距只来自数值计算的舍入误差。实测下来，各母线电压幅值最大误差在1e-9 pu量级，相角最大误差在1e-8 rad量级，基本就是浮点数精度极限。

这个实验的意义在于验证代码本身没有漏洞。如果你的程序在无噪声情况下都不能复现潮流解，那肯定不是估计方法的问题，而是H矩阵、h函数或者权重矩阵哪里写错了。所以遇到结果不对，我第一个建议就是先跑这个理想量测测试。

5.2 加入高斯量测噪声后的鲁棒性表现

接下来是更有实际意义的实验：在量测上叠加高斯噪声，标准差按前面说的方案设置，电压幅值0.005 pu，相角0.01 rad，功率0.02 pu。加入噪声后，WLS估计结果不再和潮流解重合，而是呈现统计意义上的接近。

我在IEEE 14节点系统上取一个典型运行断面，得到母线电压幅值和相角的对比结果如下表：

从结果看，幅值估计误差集中在正负0.001 pu以内，相角估计误差集中在正负0.06°以内。这个精度对调度运行来说已经相当不错了。你可能会问，既然量测噪声标准差有0.005 pu，为什么估计误差比噪声小这么多？这就是WLS“加权平均”的功劳，多条量测相互冗余，把随机噪声给平均掉了。

5.3 误差来源与量测冗余的关系

误差能压到多低，取决于两个因素：单条量测本身的精度，以及量测冗余度。用通俗的话讲，一条量测是“一个人的意见”，十条量测就是“十个证人证词”，十个证词综合下来比单个证词可靠得多。

在只接入PMU电压幅值和相角量测、没有功率量测的情况下，因为每个状态量至少有一条直接量测覆盖，WLS照样能估计，但冗余度低，误差偏大。当我把注入功率量测也加进去，多条量测相互交叉约束，估计误差就明显缩小。所以工程上很少只靠PMU做状态估计，通常是把PMU和SCADA量测混合使用，让它们互相校验。

另一个值得注意的现象是：相角估计误差整体比幅值误差小。原因有两条，一是PMU的相角量测精度高（0.01 rad对应0.57°），加上我们对相角量测设置了相对较小的方差；二是相角在H矩阵中的行就是简单的单位向量，不经过三角函数放大，数值稳定性更好。这其实就是PMU在状态估计里最直观的优势体现。

6. 实测中的避坑经验：权重设置、可观测性与收敛判断

6.1 权重矩阵R的标定远比你想的重要

很多第一次写WLS的人，把R简单设成单位阵，结果发现估计结果和最小二乘没有区别，而且坏数据的影响非常大。实际上R的对角元素直接决定每条量测在目标函数里的“话语权”，如果你把PMU相角量测的权重设得太小，它就会被一堆低精度的功率量测淹没，PMU带来的精度优势就体现不出来。

我的经验是，权重设置要按照仪表精度等级来。工程上PMU相角量测误差一般在0.1°到0.5°之间，折算成弧度大约0.002到0.009 rad；电压幅值量测精度在0.1%到0.5%之间，对应0.001到0.005 pu。SCADA功率量测误差则相对大一些，按满量程2%左右估计。把方差设成这些值的平方，基本就是合理的起点。

如果你设置的权重之间差距超过1e6倍，增益矩阵G的条件数会变得很大，数值上容易出现振荡甚至不收敛。这种情况在小规模系统上还能勉强算出来，遇到几百条母线的实际系统，基本必炸。所以我的建议是：一开始用比较保守的方差，跑通以后再慢慢调精。

6.2 相角参考基准的统一

PMU测量的是绝对相角，它依赖GPS时标，全网有一个统一的参考坐标系。但Newton-Raphson潮流解里的相角是相对于平衡节点的相对相角，平衡节点通常取θ=0。这两个参考基准如果不统一，直接比较会整体偏移一个常数角度，整个对比结果全部作废。

我在代码里的处理方法是：先选择1号母线作为平衡节点，然后对所有PMU相角量测都做一次相对化处理，即：

θᵢᵐ_rel = θᵢᵐ - θ₁ᵐ

这样就把PMU的绝对相角换算成相对于平衡节点的相角，再进入WLS估计器。这样做既保留了PMU高精度的相角差信息，又避开了绝对相角参考不一致的问题。

6.3 收敛判据和迭代次数

WLS迭代的收敛判据，我采用的是状态修正量的无穷范数小于阈值，即：

‖Δx‖∞ < 1e-8

有同学喜欢用目标函数J(x)的变化量来判断收敛，这在理论上也可以，但J的变化未必能及时反映状态修正情况，特别是在权重悬殊的时候，J可能变化很小了但Δx还在大踏步移动。我的建议是两种判据同时看，状态修正量为主，J的变化作为辅助参考。

迭代次数方面，IEEE 14节点系统从平启动（所有相角为0，所有电压幅值为1）开始，通常在5次左右收敛。如果你的程序超过15次还在迭代，大概率是H矩阵算错了。我建议在程序里打印每次迭代的Δx范数，观察它是不是单调下降。如果出现来回震荡不下降的情况，优先检查H矩阵的符号，尤其是θᵢⱼ的正负号，这个我在调试中至少翻了三次车。

我在实际工程项目里还有一个经验：状态估计真正花时间的往往不是WLS本身的迭代，而是量测数据清洗和权重标定。程序跑通只是第一步，把量测结构设计合理、让坏数据能被检测出来，那才是工程落地时最磨人的部分。如果你在这个项目基础上继续扩展，下一步可以试着加入坏数据检测与辨识、量测配置优化，或者把PMU数据流和SCADA数据的时延对齐问题一起考虑进来，这些都是电力系统状态估计在实际调度系统里绕不开的主题。