企业官网建设流程全解析

简介：电力系统状态估计是电网调度与运行分析的重要基础。该Matlab实现聚焦两大核心环节：基于最小二乘法的基本状态估计，以及针对采样数据中异常值的不良数据辨识与剔除，适合电力系统及其自动化方向的学生、研究人员与工程师快速理解状态估计的完整流程。资源共3个文件，以SEF.m脚本为主程序，配合2个txt数据文件（分别对应含误差与无误差的两组30节点输入数据），压缩包仅6KB，代码量精简，便于快速阅读算法逻辑并通过替换数据验证效果。目前已有2028人学习下载。借助源码可直接运行并观察最小二乘估计结果，对比含误差与无误差数据下的辨识差异，理解不良数据被检测后系统重新估计的处理方式；输出结果统一写入oStateEstimation，便于后续分析与拓展研究。

1. 先泼一盆冷水：最小二乘拟合在状态估计里不是万能药

学生阶段用 matlab 做曲线拟合时，最小二乘给我留下的印象是“无脑、好用”，数据丢进去、拟合完画个图就算交差。直到拿到这个电力系统状态估计 demo，第一次跑出来的结果差点以为自己看错了：量测里只混进一个不良数据，整个估计结果就被拉偏了一大截；而配合残差检验之后，它又能在同一套迭代框架里把这个坏点揪出来重新求解。这个项目干的活就是——对电网的 SCADA 量测做加权最小二乘状态估计，同时完成不良数据辨识。适合三类人：写电力系统方向毕设、研究 EMS 调度算法、想快速验证量测质量的开发者。跟着跑一遍之后你会发现，状态估计的核心难点根本不在“最小二乘”四个字，而在可观测性、权重设置和坏数据怎么剔。

2. 数学建模：WLS 目标函数、法方程与迭代收敛判据

2.1 状态向量与量测方程怎么搭

在打开代码之前，先把数学模型说清楚。N 个节点的电力系统，状态向量通常取：

x = [θ₂, θ₃, …, θ_N, V₁, V₂, …, V_N]ᵀ

维度是 2N-1。平衡节点的相角 θ₁ 不放进状态量，因为全网相角同时平移一个角度，量测方程完全不变；如果不挂一个参考节点，H 矩阵会缺秩，增益矩阵 G 奇异，迭代直接崩。电压幅值全部保留——虽然有的实现把平衡节点 V₁ 也当作已知量，但从量测冗余的角度讲，把 V₁ 保留更自然，反正 SCADA 里的母线电压量测是现成的，多一个约束多一分稳定。

量测方程统一写成 z = h(x) + ε，其中 ε 是量测噪声。h(x) 由四类基本量测构成：

• 支路有功 Pij = Vi² gij - Vi·Vj·(gij·cosθij + bij·sinθij)
• 支路无功 Qij = -Vi² (bij + bsh/2) - Vi·Vj·(gij·sinθij - bij·cosθij)
• 节点注入功率 Pi = Σ Pij，Qi = Σ Qij
• 节点电压幅值 Vi，极坐标下它本身就是状态量

这里 gij、bij 是支路串联导纳的实部和虚部，θij = θi - θj。对地电纳 bsh/2 在这类教学 demo 里通常很小，可以忽略，也可以放进量测方程的自导纳项里——后面讲代码时我会再说这个取舍。

状态估计和潮流计算有个本质区别：潮流问题是解一个确定的非线性方程组，而状态估计是在互有矛盾的冗余量测里求最优解。所以量测数量 m 必须大于状态量个数 n，才有“估计”可言。14 节点标准算例一般配 40 个左右量测，冗余度在 1.5 倍以上，这是能稳定迭代的前提。

2.2 WLS 目标函数与法方程的推导

加权最小二乘的出发点很直观：每个量测的“可信程度”不一样，误差方差大的量测，不该和方差小的平起平坐。目标函数写成：

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

其中 W 是对角权重矩阵，对角元取 1/σᵢ²。σᵢ 是第 i 个量测的标准差，权重取方差的倒数，意味着精度高的量测主导求解过程。这就是“加权”两个字的来源——普通最小二乘不区分量测精度，而电力系统量测的精度差异非常大，一台 0.2 级互感器和一台 2 级互感器给出的数据，权重差出两个数量级很正常。

对 J(x) 求梯度并令其为零，得到法方程：

Hᵀ·W·H·Δx = Hᵀ·W·(z - h(x))

其中 H 是 m×(2N-1) 的雅可比矩阵，G = Hᵀ·W·H 叫增益矩阵或信息矩阵。14 节点小系统的 G 也就是 27×27，matlab 里直接用左除求解，不需要求逆，速度极快；日后换成几百节点的配网模型，再用稀疏分解，差别也不过是毫秒级。

提示：工程上切忌直接写 inv(G)·HᵀW。因为 G 的条件数稍差时，inv 会把数值误差成倍放大，左除 G \ b 内部走的是 LU 分解，数值稳定性好一个数量级。

2.3 主迭代循环的 matlab 骨架

% 状态估计主迭代循环 % x0: 平启动初值，电压幅值1.0，相角0 % z: 量测向量，m×1 % W: 权重矩阵，对角元为 1/sigma^2 % bus: 节点数组 % branch: 支路数组，列=[首端,末端,R,X,对地电纳/2] x = x0; for iter = 1:max_iter h = compute_meas_function(x, bus, branch); % 量测函数预测值 h(x) H = compute_jacobian(x, bus, branch, Y); % 雅可比矩阵 G = H' * W * H; % 增益矩阵 dx = G \ (H' * W * (z - h)); % 法方程求解增量 x = x + dx; % 更新状态量 if max(abs(dx)) < tol break; end end

代码逻辑：每一轮先用当前状态计算量测预测值 h(x)，再求雅可比 H；法方程是线性方程组，matlab 左除一次解出 Δx；收敛判据取状态增量最大绝对值小于 tol，这是状态估计里最主流的判据，比看残差更直接。

参数建议：tol 取 1e-6 就够用，14 节点算例一般 5 到 8 轮收敛；max_iter 取 10 到 15，防止病态算例死循环。初值 x0 用平启动，全部电压幅值 1.0、相角 0。这个初值就算离真解远，只要量测配置正常，也能在几轮内拉回来。

2.4 量测权重怎么定

权重设置是新手最容易迷惑的环节，我也踩过。SCADA 量测精度的经验值：电压幅值量测方差取 (0.005)² 左右，功率量测方差取 (0.01)² 左右——标幺值体系下功率基准通常取 100MVA，量纲对得上。W 直接取 diag(1./sigma_2)，注意不要搞反，我见过有人把 σ 本身当权重填进去，结果高精度量测权重反而最小，估计出的电压幅值偏到天际。

另一个实用的校准技巧：先拿标准量测集跑一遍，观察各量测的标准化残差。如果某条支路功率量测的标准化残差总是大于 3，大概率不是量测坏，而是那个支路的参数（R、X 或变比）写错了。这时候改参数文件，比在权重里做文章更有效。

3. 项目结构与数据流：从量测文件到估计结果的调用链

3.1 目录文件清单与职责

这类 matlab 状态估计 demo 的文件组织方式很固定，核心文件一般不超过 8 个。我拆出来的这份结构如下：

其中 meas_data.m 值得多说一句：做状态估计验证时，最标准的套路是先跑一次潮流计算得到真解，在真解上叠加高斯噪声生成 SCADA 量测。这样我们手里有真值，才能量化评估估计误差。很多论文里的“仿真量测”，本质上都是这个思路。

3.2 量测数据矩阵的格式约定

量测矩阵怎么组织，直接决定后面代码的可读性。这类项目比较推荐五列格式：

% 量测矩阵格式: [类型, 首端节点, 末端节点, 量测值, 方差] % 类型编码: % 1 = 支路有功 Pij % 2 = 支路无功 Qij % 3 = 节点注入有功 Pi % 4 = 节点注入无功 Qi % 5 = 节点电压幅值 Vi（末端节点列填0） meas = [ 1, 1, 2, 0.8883, 1.0e-4; % 支路1-2有功 2, 1, 2, 0.4267, 1.0e-4; % 支路1-2无功 3, 2, 0, 0.4561, 2.5e-4; % 节点2注入有功 5, 1, 0, 1.0082, 2.5e-5; % 节点1电压幅值 ];

这种格式的好处是主程序遍历量测时，只要读 type 字段就能分派到不同的量测函数分支，不用为每个量测类型单独建数组。末端节点列对于节点注入和电压量测是无效字段，约定填 0，防止出现未初始化的随机值。

参数注意三点：量测值必须全部是标幺值；方差是标幺值的平方；如果后续接实际 SCADA 工程数据，先把有名值除以对应基准值转成标幺值再放进矩阵，否则数值量级越界，迭代会异常缓慢甚至发散。

3.3 主流程调用链

run_main.m 的调用链很短，但每一步都不可省：

% 主流程：数据准备 -> 导纳矩阵 -> 状态估计 -> 不良数据辨识 bus = case_bus(); % 读母线数据 branch = case_branch(); % 读支路参数 z_meas = meas_data(); % 读量测数据 Y = build_ybus(branch); % 节点导纳矩阵 x0 = init_flat_start(bus); % 平启动初值 x_hat = wls_se(z_meas, Y, x0); % 加权最小二乘估计 bad = bad_data(x_hat, z_meas); % 不良数据辨识 if ~isempty(bad) % 剔除坏量测后必须重新估计 x_hat = wls_se(z_meas_clean, Y, x0); end

这段流程里最容易被忽略的是“剔除不良数据之后要重新估计”。第一次估计用的坏量测已经污染了状态量，如果只把坏数据打个标记输出而不重新求解，后面所有分析、比对都建立在污染后的状态上。把剔除后的重新求解连带做掉，输出的才是可用的干净结果。

3.4 为什么不直接调 matpower

既然 matpower 自带 run_se 函数，为什么还要拆这个 demo？matpower 封装得确实完整，但同时也像个黑匣子：你很难看清雅可比矩阵怎么组装，也没法在每一轮迭代里打断看残差变化。而这个 demo 把迭代循环、雅可比组装、残差检验全部摊在脚本里。对想做二次开发的开发者来说，裸实现的改写成本比反编译 matpower 低得多——换一个估计准则、加一层抗差权重，动手改也就几十行代码。

4. 核心模块拆解：导纳矩阵、雅可比与不良数据辨识的代码实现

4.1 节点导纳矩阵构建

所有状态估计计算都建立在 Y 矩阵上。下面这段代码用稀疏矩阵构建 Y 阵：

function Y = build_ybus(branch) % branch: 每行一条支路，列 = [首端, 末端, R, X, 对地电纳/2] nb = max(max(branch(:,1:2))); Y = sparse(nb, nb); for k = 1:size(branch,1) f = branch(k,1); t = branch(k,2); R = branch(k,3); X = branch(k,4); y = 1 / (R + 1j*X); % 串联导纳 Y(f,f) = Y(f,f) + y; % 首端自导纳 Y(t,t) = Y(t,t) + y; % 末端自导纳 Y(f,t) = Y(f,t) - y; % 互导纳 Y(t,f) = Y(t,f) - y; end end

关键在三点：一是用 sparse 预分配，14 节点看不出差别，换几百节点的配网算例，稀疏和稠密求解速度差出几十倍；二是每条支路要同时更新自导纳和互导纳，漏掉任何一个，H 矩阵算出的结果都对不上真值；三是 R 和 X 直接传入，用复数导纳 y = 1/(R+jX)，规避了分别算电导电纳时容易出现的正负号错误。

对地电纳在这个教学版本里先不纳入——带上的话，每个 Q 量测方程得多加一项 Vi²·bsh/2 的偏导，提升的精度却很有限，因为对地充电功率在 14 节点规模下通常只有无功量测噪声的零头。你的数据里如果有明显的对地电纳，正确做法是在量测函数里显式补上这一项，而不是把它塞进 Y 矩阵里指望自动生效。

4.2 支路量测方程的雅可比偏导

雅可比矩阵是迭代的核心，以支路有功 Pij 对状态量的偏导为例：

∂Pij/∂θi = Vi·Vj·(gij·sinθij - bij·cosθij) ∂Pij/∂θj = -∂Pij/∂θi ∂Pij/∂Vi = 2·Vi·gij - Vj·(gij·cosθij + bij·sinθij) ∂Pij/∂Vj = -Vi·(gij·cosθij + bij·sinθij)

注意 gij、bij 是支路串联导纳的实虚部，不是 Y 矩阵对角元素的实虚部。组装时要用 Y(f,t) 那个互导纳元素来取——这也是新手最容易写错的地方。

function H = compute_jacobian(x, bus, branch, Y) % x: 状态向量 [相角(除平衡节点); 电压幅值(所有节点)] nb = size(bus,1); theta = [0; x(1:nb-1)]; % 平衡节点相角补0 V = x(nb:end); H = sparse(n_meas, 2*nb-1); for k = 1:size(branch,1) f = branch(k,1); t = branch(k,2); g = real(Y(f,t)); % 串联导纳实部 b = imag(Y(f,t)); % 串联导纳虚部 df = V(f)*V(t)*(g*sin(theta(f)-theta(t)) - b*cos(theta(f)-theta(t))); dt = -df; dVf = V(t)*(g*cos(theta(f)-theta(t)) + b*sin(theta(f)-theta(t))); dVt = V(f)*(g*cos(theta(f)-theta(t)) + b*sin(theta(f)-theta(t))); % 填入H矩阵对应行列 H(row_k, col_theta(f)) = df; H(row_k, col_theta(t)) = dt; H(row_k, col_V(f)) = dVf; H(row_k, col_V(t)) = dVt; end end

代码逻辑：对每条支路，把四个偏导分别填到 H 矩阵的对应列。col_theta 和 col_V 是两个索引映射表——相角列索引是节点编号减 1，因为平衡节点相角不进入状态量；电压列索引是 nb 加上节点编号再减 1。H 矩阵的列跨度必须和状态向量定义严格一致，映射错一位，法方程就张冠李戴。这是状态估计代码里最容易翻车的地方，调试时可以打印 H 某一行的非零元，人工核对支路两端节点号。

4.3 增益矩阵与迭代求解的数值细节

% 法方程求解 G = H' * W * H; % 检查G的条件数，条件数过大会导致求解精度崩掉 if condest(G) > 1e12 warning('增益矩阵病态，请检查量测配置是否充足'); end dx = G \ (H' * W * (z - h)); x = x + dx;

这里加了一个 condest 检查，这是我的血泪经验总结出来的：某次我少配了一条支路的功率量测，节点 5 成为事实上的孤岛，H 缺列，G 的条件数直接爆到 1e14，dx 解出来是 1e4 量级，状态量一轮迭代就飞到离谱位置。加了检查之后，问题在迭代早期就能暴露，不用等结果跑飞再回头翻数据。

1e12 这个阈值是有依据的：IEEE 双精度浮点约 16 位有效数字，线性方程求解会损失大约 log10(cond) 位有效数字，条件数到 1e10 就该警惕，1e12 基本是“还能算但结果不可信”的红线。小算例的 G 正常在 1e3 到 1e5 之间。

4.4 不良数据辨识：最大标准化残差法

这是整个项目最有价值的部分。估计完成之后，残差 r = z - h(x̂) 并不能直接比较，因为不同量测的残差方差差别很大，必须标准化：

function bad_idx = bad_data_identify(x_hat, z, H, sigma) % 计算标准化残差并标记不良数据 r = z - compute_meas_function(x_hat); % 残差向量 R = diag(sigma.^2); % 量测噪声协方差 Omega = R - H * ((H' * (R \ H)) \ H'); % 残差协方差矩阵 r_norm = abs(r) ./ sqrt(diag(Omega)); % 标准化残差 thresh = 3.0; % 3σ准则 bad_idx = find(r_norm > thresh); end

Omega = R - H·G⁻¹·Hᵀ 是残差协方差矩阵的理论表达式，它的对角元开根号就是每个残差的标准差。标准化残差服从近似标准正态分布，超过 3.0 的量测，判定为不良数据的置信度很高——正常量测残差超过 3σ 的概率不到 0.3%，而坏数据会把残差推到 5σ 以上。

实现上注意两点：一是 H'·(R\H) 用来拼 G，再用左除求逆，避免显式 inv；二是这里的 H 是最后一次迭代的雅可比，不是初始雅可比，直接用初始 H 算出来的 Omega 会偏乐观。剔出 bad_idx 后，标准流程是一次只剔一个，下面第 5 章会详细讲为什么要这么迂回。

5. 避坑指南：可观测性、迭代发散与残差污染的四个场景

5.1 G 矩阵奇异，dx 直接算出 NaN

现象：wls_se 里报 warning，矩阵接近奇异，dx 全是 NaN，估计结果崩掉。

原因：量测配置不足导致可观测性问题。我有一次想当然地删掉了某条支路的所有 P、Q 量测，以为少一条线不影响全局。结果那个支路末端的节点没有任何量测关联，H 矩阵对应位置缺列，G 奇异。

解决：补回量测，或者给该节点加一对伪量测——注入功率 0、方差取大一点（比如 0.1），把矩阵撑满。从那以后，我每次换算例的第一件事就是打印 rank(full(G))，确认满秩再往下跑。这一步只要 1 秒钟，能省下半天排查时间。

5.2 状态量迭代震荡，几十轮不收敛还越跑越歪

现象：状态量在真解附近来回跳动，残差一会儿好一会儿差，max_iter 跑满还是不收敛。

原因：两类情况最常见。一是初值离真解太远，尤其当线路重载、相角差较大时，支路方程的非线性更强，平启动会震荡。二是存在“电压控制孤岛”式的量测薄弱区域，某一片区域的量测冗余度刚好卡在临界，每轮迭代都在小幅修正之间摇摆。

解决：给更新加阻尼，x = x + alpha·dx，alpha 取 0.5 到 0.8；或者用潮流结果做初值，把真解附近作为起点。检查时把前 5 轮的 dx 范数打印出来，正常应该单调递减，如果出现 1e-1 → 1e-3 → 1e-1 的锯齿状，就是初值或阻尼的问题。

5.3 结果能收敛但明显偏真值，残差却不大

现象：估计出的电压幅值和潮流真值差 0.02 以上，功率估计也偏，但残差检验没报警。

原因：八成是权重矩阵填反了。把 σ 当成权重直接填进去，或者把方差 σ² 和标准差 σ 搞混，电压量测在 W 里的相对权重就会差出 4 个数量级，整个解被功率量测主导。这属于典型的“看着没问题，结果错得离谱”，比显式的报错难查得多。

解决：核对 W 对角元，1/σ² 是权重不是 σ。顺手检查量测方差的单位——很多实际工程数据里功率方差写的是有名值 MW，忘了除以基准值转标幺值，数量级直接差出 1e4。这两处排掉后，估计会立刻回到正常水平。

5.4 存在多个不良数据时，剔除一个后另一个“消失”了

现象：第一次剔除后重新估计，之前超阈值的另一个量测残差反而降到了阈值以下，漏检。

原因：这就是残差污染——一个大的坏数据会把相邻量测的残差“吸”过去，盖住了旁边正常量测的真实表现。如果一次把所有超阈值的量测全删了，很可能误删正常数据；但反向极端——只删一个就被“掩盖”的另一个则要迭代多次才能暴露。

解决：坚持逐条剔除策略。每轮只删标准化残差最大的那一个量测，然后重新估计、重算残差，循环到没有超阈值为止。工程经验是，在冗余度低于 1.3 的区域，可以结合量测的物理值域做初筛，比如功率量测超过线路热稳定极限，直接判死，不需要跟残差较劲。

5.5 同一份数据，别人四轮收敛你跑二十轮

现象：结果最终能收敛，但迭代次数明显偏多，每轮 dx 都只降一点点。

原因：收敛阈值设得太严。有人觉得 tol 越小越精确，直接设 1e-12。但双精度浮点下，状态量量级在 1e-1 到 1e0，dx 降到 1e-10 就进入机器精度抖动区，再小只是空转。另外每轮都重新计算 H 矩阵，虽然 14 节点无所谓，但换大算例就白白多花时间。

解决：tol 取 1e-6 足够工程使用。想提速的话，在迭代收敛附近固定雅可比矩阵，用修改牛顿法，后续每轮只做一次前代回代，速度能再快一倍。严谨起见，修改牛顿法会牺牲少量收敛速度，适合那些每轮 H 计算代价高于前代回代的场景。

6. 二次开发：从 14 节点算例到自定义量测数据的改造路径

6.1 换算例的最小改动集

自定义算例只需要动两个数据文件。第一步，把 case_bus.m 和 case_branch.m 里的母线数、支路数、R/X/B 换成新网络；第二步，按新拓扑重建量测矩阵，这里有个硬性冗余度要求——每个非平衡节点至少有一条注入量测或电压量测，每条支路最好有一对 P、Q 量测；第三步，跑一遍 rank(full(G)) 确认满秩。这三步走完，主流程一行不用改。

6.2 从 WLS 扩展到抗差估计

想抵抗少量坏数据干扰，最常见的改造是把固定权重换成迭代重加权，IRLS。在每轮迭代后用残差更新权重，Huber 函数是工程上默认的起点：

% IRLS：根据上一轮残差更新权重，抵抗少量坏数据 c = 1.345; % Huber 阈值常数 for i = 1:m ri_norm = abs(r(i)) / sigma(i); if ri_norm <= c w_i = 1; else w_i = c / ri_norm; % 大残差自动降权 end W(i,i) = w_i / sigma(i)^2; end

逻辑不复杂：残差小的量测保持原权重，残差超过阈值的量测权重线性衰减，这样坏数据对目标函数的拉动被压住，状态量不会跑偏。这个改造大概 20 行代码，却能显著提升对粗差数据的鲁棒性。效果验证时，可以故意注入两个相邻支路的坏数据，对比 WLS 和 IRLS 谁更稳。

6.3 估计结果怎么量化验证

改造完成后，用三个指标评价质量：最大绝对误差、平均绝对误差、以及不良数据辨识准确率。最好线性画出量测值和估计值的散点图，偏离 45 度线的点就是问题量测。这套验证流程比单纯看残差是否收敛更有说服力，论文和工程报告里都会用到。

从那以后，我每次拿到新的量测案例，都会强制走一遍三件事：先看 rank(G) 是不是满秩，再打印前两轮 dx 的范数判断收敛趋势，最后把标准化残差排名前十的量测列出来扫一眼。这三步做完，状态估计代码基本不会翻车。希望帮到你。

本文还有配套的精品资源，点击获取