基于MATLAB的电力系统状态估计：WLS与不良数据辨识实战

简介：这份资源面向电力系统分析与电气工程专业学习者，聚焦状态估计这一核心课题，提供基于最小二乘法的完整MATLAB实现方案。最小二乘是电力网系统辨识中最基础的状态估计方法，而实际量测数据往往夹杂不良数据，因此资源进一步引入不良数据辨识环节：先检测异常数据点位置，再将其剔除后重新进行状态估计，形成闭环处理流程。包内共3个文件，包含2个txt数据文件与1个m脚本文件，压缩包约6KB，其中两份导入数据分别对应有误差与无误差场景，便于对比验证算法效果，估计结果统一输出至指定变量中。目前已有2029人学习下载，说明该案例在电力系统状态估计教学中具有较高的参考价值。读者可借此掌握最小二乘状态估计的建模思路、不良数据检测与剔除的实现细节，并直接运行脚本观察不同数据条件下的估计结果差异，适合课程设计、毕业设计及科研入门阶段参考使用。

1. 电力系统状态估计：从量测噪声里把真实运行状态“抠”出来

电网调度中心屏幕上那组母线电压、相角数据，看着挺稳，其实背后全是量测噪声和坏数据。SCADA 传上来的遥测值，偶尔会蹦出一个明显偏离物理规律的数——可能是 CT 饱和，可能是通道干扰，也可能是通信误码。如果直接拿这些数去做潮流计算或者安全分析，结果轻则偏差大，重则直接发散。这个 MATLAB 项目干的事，就是电力系统状态估计：用最小二乘法从一堆带噪声的量测里，估计出最接近真实运行状态的电压幅值和相角，同时把那些“捣乱”的不良数据揪出来。适合电力系统自动化、调度自动化方向的开发者，或者正在做课程设计、想搞懂 WLS 和不良数据辨识到底怎么落地的人。整套代码在 MATLAB 里跑，不需要额外工具箱，拿到手改改量测数据就能用。

2. 核心架构拆解：量测方程、WLS 迭代与残差检测怎么串起来

2.1 状态估计的数学骨架：量测方程 z = h(x) + v

电力系统状态估计的本质，是找一个状态向量 x，让量测向量 z 和由 x 算出来的理论量测 h(x) 之间的差距最小。状态向量 x 通常取所有母线的电压幅值和相角（平衡节点相角作为参考，不参与估计）。量测 z 包括支路有功/无功潮流、母线注入有功/无功、母线电压幅值。量测方程 h(x) 就是潮流方程那套非线性关系。

项目里用到的核心迭代公式是加权最小二乘（WLS）的修正方程：

% 计算量测函数 h(x) 和雅可比矩阵 H [H, h] = compute_h_and_H(x, Ybus, measurement_config); % 计算残差向量 r = z - h; % 增益矩阵 G = H' * W * H G = H' * W * H; % 右端项 H' * W * r rhs = H' * W * r; % 求解修正量 delta_x delta_x = G \ rhs; % 更新状态量 x = x + delta_x;

这段代码是整个估计器的“心脏”。compute_h_and_H根据当前状态 x 算出理论量测值和雅可比矩阵；W是量测权重矩阵，通常取量测误差方差的倒数，对角阵。G是增益矩阵，维度是状态数乘状态数，规模不大（几十到几百节点），直接求逆或者用\解线性方程组都行。delta_x是状态修正量，迭代到max(abs(delta_x))小于收敛阈值（比如 1e-6）就停。

参数说明：Ybus是节点导纳矩阵，由网络参数生成；measurement_config定义了哪些量测存在、对应哪条支路或哪个节点；W的对角元素直接决定估计结果更“信”哪个量测。常见做法是电压量测权重给高一点，潮流量测权重按精度等级给。

2.2 雅可比矩阵的稀疏组装：别用循环硬拼

雅可比矩阵 H 的维度是“量测数 × 状态数”，电力系统里量测数通常是状态数的 1.5 到 2 倍。如果每个元素都用 for 循环去填，节点数一多就慢得离谱。项目里用的是稀疏组装方式，先按支路和节点分别计算局部雅可比块，再塞进全局稀疏矩阵。

% 初始化稀疏矩阵的三元组 rows = []; cols = []; vals = []; for k = 1:num_branches % 计算支路潮流对两端母线电压幅值和相角的偏导 [dP_dVa, dP_dVm, dQ_dVa, dQ_dVm] = branch_jacobian(branch_data(k), V, theta); % 把局部偏导填入全局行列索引 % ... 省略索引映射细节 end % 用 sparse 一次性组装 H = sparse(rows, cols, vals, num_measurements, num_states);

逻辑说明：branch_jacobian返回的是 2×4 或 4×4 的小块，对应有功/无功对电压幅值/相角的偏导。索引映射要把支路两端母线在状态向量里的位置找出来。最后用sparse组装比动态扩容快一个数量级。参数方面，V和theta是当前迭代的电压幅值和相角向量，branch_data包含支路阻抗、对地导纳等。

2.3 不良数据辨识：残差范数与假设检验

WLS 估计完之后，不能直接信结果，得看残差。如果某个量测的残差特别大，它很可能就是不良数据。项目里用的是标准化残差法：先算残差协方差矩阵，再算每个量测的标准化残差，超过阈值（通常 3σ）就标记为可疑。

% 计算残差灵敏度矩阵 S = eye(num_measurements) - H * (G \ (H' * W)); % 残差协方差矩阵 Omega = S * inv(W) * S'; % 标准化残差 r_norm = abs(r) ./ sqrt(diag(Omega)); % 找出超过阈值的量测 bad_idx = find(r_norm > 3.0);

这里S是残差灵敏度矩阵，反映量测误差如何传递到残差。Omega的对角线是各量测残差的方差。r_norm就是标准化残差，大于 3 的认为是可疑不良数据。实际用的时候，通常会迭代：先估计，再检测，把最大的不良数据剔除或降权，重新估计，直到没有超限的。

注意：阈值 3.0 不是绝对的，量测冗余度低的时候要放宽到 3.5 甚至 4.0，否则容易误杀正常量测。

3. 跑起来：从数据准备到结果输出的完整操作链

3.1 环境与文件清单

项目在 MATLAB R2018b 及以上版本测试通过，不需要任何工具箱。把压缩包解压后，目录结构大致如下：

常见做法是先把data/里的算例换成自己的电网数据，再跑main.m。如果只是验证算法，直接用自带的 IEEE 9 节点或 14 节点算例就行。

3.2 主流程逐段拆解

main.m里的流程很直白，按顺序调用几个核心函数：

% 1. 加载算例数据 case_data = load_case('data/ieee14.m'); % 2. 生成节点导纳矩阵 Ybus = makeYbus(case_data); % 3. 生成量测数据（含噪声和不良数据） [z, measurement_config] = generate_measurements(case_data, Ybus, noise_level, bad_data_ratio); % 4. 初始化状态量（平启动） x0 = init_state(case_data); % 5. WLS 状态估计 [x_est, converged, iter] = wls_state_estimation(z, measurement_config, Ybus, x0); % 6. 不良数据辨识 [bad_idx, r_norm] = bad_data_detection(z, x_est, measurement_config, Ybus); % 7. 剔除不良数据后重新估计 if ~isempty(bad_idx) z_clean = z; z_clean(bad_idx) = []; measurement_config_clean = remove_measurements(measurement_config, bad_idx); [x_est, ~, ~] = wls_state_estimation(z_clean, measurement_config_clean, Ybus, x_est); end % 8. 输出结果 print_results(x_est, case_data); plot_results(x_est, case_data, r_norm);

逻辑说明：第 3 步的generate_measurements是关键，它模拟了真实量测的噪声和不良数据。noise_level控制高斯噪声标准差，bad_data_ratio控制不良数据比例。第 5 步的wls_state_estimation就是第 2 章讲的迭代过程。第 6 步检测出不良数据后，第 7 步把它们剔掉再估一次，结果会干净很多。

参数怎么改：noise_level默认 0.01，想测试算法鲁棒性可以调到 0.05；bad_data_ratio默认 0.02，调大能看检测算法的极限。收敛阈值在wls_state_estimation里，默认 1e-6，一般不用动。

3.3 结果怎么看：电压幅值对比与残差图

跑完之后会弹两个图。第一个是估计电压幅值和真实值的对比，横轴是母线编号，纵轴是标幺值。如果估计曲线和真实曲线基本重合，说明估计器工作正常。第二个是标准化残差分布图，横轴是量测编号，纵轴是标准化残差，超过红色阈值线的点就是被标记的不良数据。

常见做法是先把bad_data_ratio设为 0，看纯噪声下的估计精度；再逐步调大，观察不良数据检测的召回率和误报率。如果发现某个正常量测老是被误判，检查一下它的权重是不是设得太低了。

4. 避坑与排查：那些让估计结果“翻车”的细节

4.1 现象：迭代不收敛，delta_x 越算越大

原因：最常见的是量测配置和网络拓扑不匹配。比如某条支路的量测存在，但支路本身在导纳矩阵里被设成了断开；或者参考节点选错了，导致相角不可观测。另一个原因是权重矩阵 W 里某个量测权重设成了 0 或者负数，增益矩阵直接奇异。

解决：先检查measurement_config里每条量测对应的支路和节点是否都在case_data里存在。再把 W 的对角元素打印出来，确保全是正数。如果还不行，把收敛阈值放宽到 1e-4 试试，有时候是数值精度问题。

4.2 现象：不良数据检测把所有量测都标红了

原因：残差协方差矩阵Omega计算时用了inv(W)，如果 W 接近奇异（比如某些量测权重极小），Omega的对角线会异常大或异常小，导致标准化残差失真。另一个可能是量测冗余度太低，残差灵敏度矩阵 S 接近单位阵，残差直接反映量测误差，没有“平滑”效果。

解决：把 W 的条件数打出来，如果大于 1e12，说明权重设置有问题，把最小权重调大一点。量测冗余度低的话，阈值从 3.0 放宽到 4.0，或者改用最大标准化残差法，每次只剔除最大的那个，迭代几次。

4.3 现象：估计结果和真实值偏差很大，但残差看起来正常

原因：这种情况通常是“杠杆量测”在作怪。某些量测对状态估计的影响特别大，但本身误差也大，WLS 会被它带偏，而残差因为灵敏度矩阵的“掩盖”效应反而不大。常见于注入功率量测配置不合理的时候。

解决：检查注入功率量测的权重是不是给得太高了。常见做法是把注入功率量测权重设得比支路潮流低一个数量级。另外可以用最大标准化残差法替代整体检测，对杠杆量测更敏感。

4.4 现象：MATLAB 报错“矩阵维度不一致”

原因：z向量和measurement_config的长度不匹配，通常是剔除不良数据后忘了同步更新配置。或者H矩阵组装时行列索引超出了预分配的大小。

解决：在每次修改z之后，立刻检查length(z) == length(measurement_config)。组装 H 之前，先算一下num_measurements和num_states，确保sparse的维度参数对得上。

提示：调试的时候把dbstop if error打开，报错时直接停在出错行，看工作区变量比猜快得多。

5. 二次开发：换数据、改权重、加检测算法

5.1 换成自己的电网数据

项目自带的算例是标准 IEEE 数据，格式是 MATPOWER 风格的。要换成自己的数据，最省事的做法是照着data/ieee14.m的格式改。需要准备三样东西：节点数据（母线编号、电压等级、负荷）、支路数据（首末节点、阻抗、对地导纳）、量测配置（量测类型、位置、权重）。

% 自定义算例模板 case_data.baseMVA = 100; case_data.bus = [ 1 3 0 0 0 0; % 节点1，平衡节点 2 1 50 20 0 0; % 节点2，PQ节点，负荷50MW+20MVar % ... 更多节点 ]; case_data.branch = [ 1 2 0.01 0.05 0.02 0; % 支路1-2，阻抗0.01+j0.05，对地导纳0.02 % ... 更多支路 ];

改完之后，把main.m里的load_case参数换成你的文件名。量测配置可以先用全量测（每条支路两端都有潮流，每个节点都有注入和电压），跑通了再逐步删减，看最少需要多少量测才能保证可观测性。

5.2 改权重矩阵：让估计器更“信”谁

权重矩阵 W 直接决定估计结果。默认是量测误差方差的倒数，但实际工程里，不同量测的精度不一样。电压量测通常比潮流量测准，所以权重给高一点。注入功率量测因为聚合了多个负荷，误差大，权重给低一点。

% 按量测类型设置权重 W = speye(num_measurements); for i = 1:num_measurements switch measurement_config(i).type case 'V' % 电压幅值 W(i,i) = 1 / (0.005^2); case 'P' % 支路有功 W(i,i) = 1 / (0.01^2); case 'Q' % 支路无功 W(i,i) = 1 / (0.01^2); case 'Pinj' % 注入有功 W(i,i) = 1 / (0.03^2); end end

这段代码把电压量测的标准差设为 0.005，支路潮流 0.01，注入功率 0.03。标准差越小，权重越大，估计器越“信”这个量测。改的时候注意别把某个权重设得比别的量测高好几个数量级，否则增益矩阵条件数爆炸，迭代容易出问题。

5.3 加一个自己的不良数据检测算法

项目里用的是标准化残差法，你可以换成假设检验法、量测突变检测法，或者基于机器学习的异常检测。接口很简单：输入是z、x_est、measurement_config、Ybus，输出是bad_idx和r_norm。只要保持这个签名，直接替换bad_data_detection.m就行。

function [bad_idx, score] = my_detector(z, x_est, measurement_config, Ybus) % 计算理论量测值 [~, h] = compute_h_and_H(x_est, Ybus, measurement_config); % 计算残差 r = z - h; % 用马氏距离替代标准化残差 % ... 计算残差协方差矩阵的逆 score = sqrt(diag(r' * inv(Omega) * r)); % 卡方检验阈值 threshold = chi2inv(0.99, 1); bad_idx = find(score > threshold); end

逻辑说明：马氏距离考虑了残差之间的相关性，比单独看每个残差的标准化值更严格。chi2inv(0.99, 1)是自由度为 1 的卡方分布 99% 分位数，约等于 6.63。用这个阈值，误报率更低，但召回率也会降。实际用的时候可以先用标准化残差粗筛，再用马氏距离精筛。

注意：换检测算法之后，记得把main.m里对应的函数名也改掉，别只改文件内容。

5.4 一个我踩过的坑

刚开始做这个项目的时候，我觉得权重矩阵随便设设就行，反正 WLS 会自己收敛。结果有一次用了一个实际电网的数据，迭代了 50 次还不收敛，delta_x在 0.1 附近来回跳。查了一下午才发现，有个注入功率量测的权重设成了 1e6，而电压量测权重只有 1e2，增益矩阵条件数直接爆了。从那以后我每次设完权重，都强制走一遍cond(G)检查，条件数超过 1e10 就回去调权重。这个习惯帮我省了很多玄学调试时间。希望帮到你。

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