电力系统加权最小二乘状态估计与不良数据辨识实战解析

简介电力系统状态估计是电网调度与运行分析的重要基础。该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 Σ PijQi Σ 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×27matlab 里直接用左除求解不需要求逆速度极快日后换成几百节点的配网模型再用稀疏分解差别也不过是毫秒级。提示工程上切忌直接写 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 个。我拆出来的这份结构如下文件职责run_main.m主脚本按顺序调用数据加载、导纳计算、状态估计、不良数据辨识case_bus.m母线原始数据节点编号、类型、电压基准case_branch.m支路参数首端/末端节点、R、X、对地电纳meas_data.m量测数据生成在真值上叠加高斯噪声按需注入坏量测build_ybus.m由支路参数构建节点导纳矩阵wls_se.m加权最小二乘求解核心迭代体bad_data.m不良数据辨识计算标准化残差并循环剔除其中 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 函数为什么还要拆这个 demomatpower 封装得确实完整但同时也像个黑匣子你很难看清雅可比矩阵怎么组装也没法在每一轮迭代里打断看残差变化。而这个 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/(RjX)规避了分别算电导电纳时容易出现的正负号错误。对地电纳在这个教学版本里先不纳入——带上的话每个 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 的条件数直接爆到 1e14dx 解出来是 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); endOmega 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·dxalpha 取 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 到 1e0dx 降到 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 的范数判断收敛趋势最后把标准化残差排名前十的量测列出来扫一眼。这三步做完状态估计代码基本不会翻车。希望帮到你。本文还有配套的精品资源点击获取