1. 项目概述与核心思路拆解1.1 这个项目到底在解决什么问题电力系统状态估计是EMS能量管理系统最核心的功能模块之一。调度中心每天面对海量遥测数据这些数据来自SCADA系统受CT/PT变比误差、通道噪声、测量装置老化等因素影响量测值并不精确。如果直接用这些有误差的数据去做在线分析、安全评估结果很可能失真。状态估计的作用就是利用一组带误差的冗余量测按照统计学准则估算出系统最可能处于的运行状态这个状态通常用所有节点的电压幅值和相角来表示。这个项目的标题很直白它要做的事情是搭一个可复现的IEEE标准测试系统用加权最小二乘法WLS配合PMU相量量测估算系统中各节点的电压幅值和相角同时用Newton-Raphson潮流计算求解同一系统的状态把PMU/WLS的估计结果和NR潮流的解放在一起对比检验状态估计算法的准确性。项目用Matlab实现适合电力系统方向的研究生、刚接触EMS算法的工程师以及研究PMU在电网中应用的同学参考。做完这个项目你不仅仅会运行一段代码你会把量测方程、雅可比矩阵、权重矩阵、可观测性这几块硬骨头一次啃透。之所以拿NR潮流做基准不是因为潮流比状态估计更高级而是因为潮流求解的问题本身就是确定的给定各节点注入功率方程组数量等于未知量数量解是确定的。状态估计处理的是量测冗余问题方程数多于未知数解是统计意义上的最优估计。如果把量测噪声控制得很小WLS的估计结果理论上会非常逼近NR潮流的解。所以用NR潮流解作为标准答案去验证状态估计在含噪声量测条件下的表现逻辑上说得通工程上也是常用的验证手段。1.2 PMU与传统量测的代差体现在哪里PMU即相量测量单元和传统RTU远程终端单元最本质的区别在于PMU能直接测量电压和电流的相量也就是幅值和相角。深入研究PMU在状态估计中的作用首先要理解传统量测的局限。SCADA系统里大量使用的是功率量测有功P、无功Q和电压幅值量测相角信息无法直接测到只能靠状态估计算法从功率方程里间接推算出来。而且传统RTU的扫描周期是秒级各个量测点之间没有统一的时间基准同步性差遇到系统动态过程时量测数据反映的其实不是同一个断面。PMU依靠GPS/北斗卫星授时采样率通常做到几十帧每秒所有PMU上报的数据都有统一时标天然是同步断面的。更关键的是PMU能直接给出电压相角这正好补齐了传统量测的短板。在这个项目里把PMU量测纳入WLS状态估计后相角不再完全依赖间接推算信息矩阵的形态会更好估计精度也会提升。不过需要注意PMU在工程上还没有做到全网全覆盖更现实的方案是传统量测为主、PMU量测为辅的混合量测架构本项目的量测配置也遵循这个思路。1.3 为什么WLS依然是状态估计的主流选择尽管学术界提出了不少新算法比如抗差估计M估计、最大指数绝对值估计等但WLS凭借模型简单、计算效率高、在正态噪声假设下具有最优线性无偏性质至今仍然是工程界应用最广的状态估计内核。WLS的核心思想并不复杂把所有量测方程写成一个超定方程组然后寻找一组状态变量使量测残差的加权平方和最小。权重取什么值很关键通常取量测噪声方差的倒数也就是说精度越高的量测在目标函数里的话语权越大。具体到这个项目WLS要处理三种量测支路功率线路首末端P和Q、节点注入功率P_i和Q_i、电压幅值量测再加上PMU提供的节点电压相角量测。结果就是在原有量测方程组里增加一组线性方程相角量测方程简单雅的可比矩阵那一行也简单但它在算法收敛性和解算精度上带来的提升实测下来非常明显。整个迭代求解过程本质上是Gauss-Newton法每轮迭代利用线性化的量测方程构造法方程解出状态修正量后更新直到修正量的最大绝对值小于收敛阈值。2. 状态估计与潮流计算的数学对照2.1 WLS状态估计的数学框架状态估计的数学模型普遍采用如下量测方程描述z h(x) e其中z是量测向量x是系统状态向量通常取所有节点的电压幅值和相角参考节点相角除外h(x)是量测函数表示状态到量测的非线性映射e是量测噪声向量一般假设为零均值高斯分布协方差矩阵为R。量测噪声的方差由仪表精度、变送器误差、通信量化误差等因素共同决定。目标函数定义为J(x) [z - h(x)]^T R^{-1} [z - h(x)]对J(x)求梯度并令其为零得到迭代公式Δx (H^T R^{-1} H)^{-1} H^T R^{-1} [z - h(x)]其中H是量测雅可比矩阵第i行第j列元素表示第i个量测函数对第j个状态变量的偏导数。这个式子看着吓人实际上每一步迭代都是在解一个线性最小二乘问题G H^T R^{-1} H称为信息矩阵或增益矩阵工程上也直接叫它G矩阵。当量测冗余度足够、系统可观测时G矩阵是对称非奇异的可以稳定求逆。说白了WLS做的就是在量测残差的加权平方和最小的意义下把状态量反推出来。如果量测全是PMU提供的电压相量那么h(x)就是线性的一步就能得到解。但实际系统中功率量测占据大头所以还是得走迭代。2.2 Newton-Raphson潮流计算的求解框架潮流计算的问题描述要简洁得多给定平衡节点电压幅值与相角、PV节点有功注入和电压幅值、PQ节点有功和无功注入求所有节点的电压幅值和相角。核心方程组是节点功率平衡方程P_i V_i Σ_{j} V_j (G_ij cosθ_ij B_ij sinθ_ij) Q_i V_i Σ_{j} V_j (G_ij sinθ_ij - B_ij cosθ_ij)其中G_ij和B_ij是节点导纳矩阵元素的实部和虚部θ_ij是节点i和j的相角差。NR法把这个非线性方程组在初值附近做泰勒展开忽略二阶以上项得到线性修正方程组[ΔP; ΔQ] J_NR [Δθ; ΔV]J_NR是潮流雅可比矩阵由P、Q对θ、V的偏导数构成。每轮迭代计算当前状态下的注入功率与给定注入做差得到失配量然后解线性方程组求修正量更新状态后继续迭代直到失配量足够小。NR潮流的收敛特性好正常条件下平启动也能在几次迭代内收敛所以在工程软件里基本是默认选择。这个项目里NR潮流还有一个特殊用途就是产生标准答案先用潮流解出系统真实状态再在这个真实状态上叠加高斯噪声模拟量测值这样我们就知道了量测对应的真值评估WLS估计误差才有依据。2.3 两个方法到底有什么本质区别很多初学者会把潮流计算和状态估计搞混因为两个算法都在解电压幅值和相角公式长得像迭代框架也像。但它们的本质区别非常明显。从输入数据看潮流计算用的是发电和负荷的计划值或预测值是一种模型驱动状态估计用的是实际量测是数据驱动。从方程数量看潮流是方阵求解方程数和未知数数量一致状态估计是超定求解量测方程数量远大于状态变量数量。从解的性质看潮流解是模型输入下的确定性响应状态估计解是统计意义下的最优估计它能同时给出量测残差供坏数据检测使用。这个对照反映在代码层面就是NR潮流失配量计算直接对应功率平衡方程没有权重概念WLS的右端项则是加权残差H^T W (z - h(x))包含权重矩阵。两者的雅可比矩阵也不一样潮流的雅可比是功率对电压幅值和相角的偏导状态估计的雅可比是各类型量测函数对状态的偏导维度通常比潮流大得多。在项目实现时我会把两者放在同一个脚本框架里量测雅可比和潮流雅可比分开构建避免搞混。3. Matlab实现流程与关键环节解析3.1 测试系统选型与数据准备做状态估计不能上来就在自己拼的随机网络上跑否则哪里算错都说不清。靠谱的做法是直接用IEEE标准节点系统推荐IEEE 14节点系统。它的规模适中既不会因为网络太大导致可视化困难又能体现出状态估计对全网状态的感知能力。IEEE 14节点系统包含5个发电节点、11条负荷支路、3台变压器和1条并联电容器支路量测配置可以做得比较丰富。在Matlab中我建议直接借助Matpower工具箱加载数据它可以省去手动构造线路参数和节点注入数据的麻烦mpc loadcase(case14); Ybus makeYbus(mpc); baseMVA mpc.baseMVA; bus mpc.bus; branch mpc.branch; gen mpc.gen;这里loadcase会返回一个结构体里面包含节点数据矩阵、支路数据矩阵、发电机数据矩阵。节点数据矩阵的各列含义需要弄清楚其中第8列是电压幅值标幺值第9列是电压相角角度制后续构造量测时会用到。makeYbus则根据节点和支路数据生成全系统的节点导纳矩阵这是所有后续计算的基础。有一点要强调Matpower里的负荷和发电机数据都是标幺值量测构造时要保持同样的标幺体系。如果直接拿有名值去算导纳矩阵的量级会对不上迭代计算容易出现数值发散。做项目时先在脚本开头统一单位的处理后续所有量测和状态量都用标幺值或弧度制只在输出结果时换算成有名值或角度制。3.2 量测数据生成与PMU配置方案WLS状态估计需要量测数据作为输入。在真实系统中量测来自RTU和PMU在仿真项目里需要人为构造量测。我的做法是分三步。第一步用NR潮流计算系统的精确状态得到所有节点电压幅值V_true和相角θ_true弧度制。这个状态就是后续评估的基准。第二步根据要量测的物理量和拓扑关系计算理想的量测值比如支路潮流P_ij通过功率方程由V_true和θ_true算出。第三步在理想量测值上叠加高斯噪声模拟真实量测。Matlab里生成带噪声量测的代码如下% 设定量测误差标准差标幺值 sigma_P 0.01; % 有功量测误差标准差约为1% sigma_Q 0.01; % 无功量测误差标准差 sigma_V 0.005; % 电压幅值量测误差标准差0.5% sigma_theta 0.001; % PMU相角量测误差标准差约0.057度 % 假设已通过潮流计算出P_true, Q_true, V_true, theta_true z_P P_true sigma_P * randn(size(P_true)); z_Q Q_true sigma_Q * randn(size(Q_true)); z_V V_true sigma_V * randn(size(V_true)); z_theta theta_true sigma_theta * randn(size(theta_true_pmu));PMU配置方案我推荐这样设计不是所有节点都装PMU那不符合工程实际。典型做法是在部分关键节点配置PMU比如发电机出口节点和电压中枢节点。在每个配置了PMU的节点能量测到的量包括该节点的电压幅值和相角以及与该节点相连的支路电流幅值相角。为了项目简单可以只把PMU的电压相角量测加入量测集也就是在传统量测基础上增加一组相角量测方程。量测冗余度这个概念需要重点理解。系统状态变量的数量是2n - 1n个节点每个节点有幅值和相角但参考节点相角固定量测总数必须显著大于状态数才能保证可观测性并具备抗差能力。一般工程实践中量测冗余度在2到3以上比较好。IEEE 14节点系统状态变量共27个我建议构造至少50到60个量测既包含支路潮流、节点注入、电压幅值再叠加PMU相角量测这样信息矩阵条件数才好看。3.3 WLS状态估计迭代实现核心代码WLS的核心迭代循环并不复杂但有几个细节必须处理到位否则程序很容易跑飞或不收敛。下面给出一个可以直接改用的框架function [V_est, theta_est, iter] wls_estimator(Ybus, z, W, idx, loadflow_solved, tol, max_iter) % Ybus: 节点导纳矩阵 % z: 量测向量按固定顺序组织支路有功、支路无功、注入有功、注入无功、电压幅值、PMU相角 % W: 权重对角矩阵 % idx: 结构体记录各类量测在z中的索引位置 % 初始化平启动 n size(Ybus, 1); V ones(n, 1); theta zeros(n, 1); x [theta(2:end); V]; % 状态向量参考节点相角固定为0 for iter 1:max_iter % 构建量测雅可比矩阵 H 和量测函数值 h(x) [H, h] build_measurement_jacobian(V, theta, Ybus, idx); % 增益矩阵 G H * W * H; % 右端项 rhs H * W * (z - h); % 求解修正量 dx G \ rhs; % 更新状态 dtheta [0; dx(1:n-1)]; dV dx(n:end); theta theta dtheta; V V dV; if max(abs(dx)) tol break; end end theta_est theta; V_est V; endbuild_measurement_jacobian函数是整个程序里最容易写错的地方。它要对每种量测类型分别求偏导。以支路有功量测为例P_ij对θ_i、θ_j、V_i、V_j都有偏导公式不复杂但容易漏项。一个实用的检查办法是用有限差分法在相同状态下数值求雅可比和解析表达式做对比最大误差在1e-4以下基本可以确认偏导公式写对了。信息矩阵G H^T W H的维度是(2n-1)×(2n-1)在IEEE 14节点系统上规模很小直接用Matlab的左除运算符\即可求解。如果系统规模很大则需要利用G矩阵的稀疏性进行因子分解避免对稠密矩阵求逆。IEEE 14节点系统上直接稠密求解没有性能压力。还有一个细节PMU量测的相角是绝对相角而状态估计的参考节点相角固定为0。如果PMU配置在参考节点之外的节点它的相角量测方程是θ_i加上一个系统参考偏移在算法实现上要么把参考节点同PMU对齐要么在量测方程中把参考相角作为一个附加状态变量来处理。本项目里建议直接选择参考节点也配置PMU这样处理最简洁避免了相角参考不一致的问题。3.4 Newton-Raphson潮流求解实现要点NR潮流在这个项目里的定位是真值生成器同时也作为状态估计的对照基准。用Matpower的runpf可以直接得到潮流结果但如果想深入理解算法还是建议自己写一遍核心迭代。自写NR潮流有几个关键点。第一节点分类要清晰ref find(bus(:, 2) 3); % 平衡节点 pv find(bus(:, 2) 2); % PV节点 pq find(bus(:, 2) 1); % PQ节点第二雅可比矩阵的构建。NR潮流的雅可比矩阵分成4个子块dP/dθ、dP/dV、dQ/dθ、dQ/dV。注意PV节点的无功方程不参与迭代因为PV节点的无功注入是待定量平衡节点的有功和无功方程都不参与迭代。最终雅可比矩阵的维度是(2×PQ数量 PV数量)的方阵。第三迭代修正量的更新方式。电压幅值采用相乘修正V_new V ΔV相角采用相加修正θ_new θ Δθ。收敛判据推荐使用功率失配量的最大绝对值通常设为1e-8。自己写一遍NR潮流会让你对功率平衡方程的理解从名词变成肌肉记忆后续构建量测雅可比时很多偏导数公式其实是互通的会省不少力气。3.5 结果对比与误差指标体系完成WLS状态估计和NR潮流求解后就来到了验收环节。需要对比的物理量主要是两类电压幅值V的估计值和电压相角θ的估计值。通常的控制台输出方式是把所有节点的数据按三列并排打出来真值、WLS估计值、NR潮流计算值。我习惯用如下方式组织结果err_V abs(V_est - V_true); err_theta abs(theta_est - theta_true); mean_err_V mean(err_V); max_err_V max(err_V); rmse_V sqrt(mean((V_est - V_true).^2));从实际运行结果来看当量测误差标准差设置为典型值有功1%、无功1%、电压幅值0.5%、PMU相角0.001弧度时WLS估计的电压幅值平均误差通常能控制在0.1%以内相角平均误差大约在0.01到0.05度之间。如果PMU数量增加相角估计精度会有更明显的提升。图形对比建议画两个图。第一个图画出各节点的电压幅值真值、WLS估计值和NR潮流值的三条曲线能直观看到WLS估计与真值之间的偏差大小。第二个图画出相角误差的柱状图从图上能清晰看出哪些节点相角估计误差偏大通常是没有PMU覆盖的负荷节点误差更大这是符合直觉的。4. 调试心得、收敛问题与工程避坑指南4.1 初值选择与收敛性判断WLS状态估计的初值一般取平启动也就是所有节点电压幅值取1.0相角取0。对绝大多数正常量测配置来说平启动就能收敛。但有两种情况容易出问题一是量测冗余度不够G矩阵接近奇异迭代会出现振荡不收敛二是量测中存在不良数据导致残差过大迭代方向脱离可行域。我在调试时遇到过一种典型现象程序在头几次迭代里修正量很大看起来要发散但其实是因为初值状态距离真实解太远线性化模型误差较大。处理方法很简单不要急着判定发散观察一下修正量的曲线趋势。如果修正量虽然在波动但整体在减小就让它继续迭代如果修正量持续增大或者出现NaN再考虑调整初值或检查雅可比矩阵。经验做法是先用NR潮流算一遍把潮流解作为WLS的初值。这样做的好处是初值距离真值足够近WLS迭代次数会明显减少而且能有效避免因为初值太差导致的数值问题。代价是潮流本身要先解出来但在这个项目里潮流本来就要算所以白赚一个稳定初值。4.2 权重矩阵W的设置到底有什么讲究权重矩阵W是量测噪声协方差矩阵R的逆矩阵。理论上讲W的取值应该反映量测的真实精度。但工程上量测精度并不总是准确已知权重的设置往往成了影响估计质量的关键因素。从实践来看有几条经验可以直接用。各类量测误差标准差的典型参考值大致如下量测类型典型误差标准差权重参考值支路有功P_ij0.011%10000支路无功Q_ij0.011%10000节点注入有功P_i0.011%10000节点注入无功Q_i0.011%10000电压幅值V_i0.0050.5%40000PMU相角θ_i0.0010.057°1000000注意PMU相角的权重比传统量测高出两个数量级这在数学上是合理的因为PMU相角的测量精度确实远高于由功率量测间接推算的相角。但一个隐患是如果PMU量测本身存在系统误差比如相角基准偏差过高的权重反而会把估计结果往错误方向带。所以工程上对PMU数据一般会先做合理性校验再进入状态估计器。权重矩阵的病态问题也值得留意。当权重数值跨越5到6个数量级时增益矩阵G的条件数会变得很大求解线性方程组时数值误差可能掩盖真实解。遇到这种情况可以对权重做适度缩放比如统一除以最大权重把权重范围压到10^(-3)到1之间数值稳定性会好很多。4.3 量测冗余度对估计精度的影响量测冗余度的定义是量测数量与状态数量之比。这个指标直接影响估计精度和抗差能力。我做过一组测试把量测从最少的刚好可观测冗余度接近1逐步增加到冗余度超过3对比WLS的估计误差。结果非常有意思冗余度从1提升到2估计误差下降非常显著差不多能下降50%再从2提升到3误差下降幅度变缓大概再降20%左右。这背后的道理也不难理解。冗余量测提供了更多关于系统状态的约束信息最小二乘估计本质上是多个量测信息的折中量测越多随机噪声被平滑掉的效果越好。但冗余度超过3之后新增量测带来的边际收益越来越小而量测系统的建设和维护成本却直线上升。工程上设计量测配置时要找到性价比平衡点不是量测越多越好。在IEEE 14节点系统上做实验时一个很直观的做法是先只布置最小可观测量测比如各节点电压幅值平衡节点和若干支路功率跑一遍WLS记录误差然后逐步添加支路功率量测、注入量测、PMU相角量测每次添加后重新运行记录误差变化。你会看到一个清晰曲线这也是理解量测冗余度最好的实验。4.4 不良数据检测的处理思路实际系统的量测数据难免出现坏数据比如采集模块故障发出一个特别离谱的值或者通信传输偶发丢包导致数据错乱。如果不做处理一个坏数据就足以让WLS估计结果整体偏移污染整个状态估计质量。这就是工程上为什么状态估计必须配不良数据检测模块。最常见的方法是残差分析法。WLS求解收敛后可以计算每个量测的归一化残差r_i (z_i - h_i(x_hat)) / sqrt(R_ii)在正态噪声假设下如果量测正常r_i应当服从标准正态分布大概率落在[-3, 3]区间内。如果某个量测的归一化残差显著超过3基本可以判定它是坏数据。工程实现的时候比较耗时的是剔除坏数据后的重新估计。一种策略是先把所有超限的量测全部清除然后重新跑WLS再检查残差另一种策略是每次只剔除残差最大的那个量测重跑后再检查迭代执行。后者的计算量更大但更稳健因为同时剔除多个量测可能误杀正常数据。这个项目里可以在量测集里刻意注入一两个坏数据比如把某个P_ij量测增大50%然后观察不良数据检测机制能否准确把它揪出来。这个实验做完你对状态估计的工程细节会有一个质的理解提升。5. 从仿真到工程现场的延伸思考5.1 混合量测状态估计的工程意义前面所有工作都是在一个仿真环境里完成的。但仿真项目的价值不只是交一份报告更重要的是理解工程实现背后的逻辑。当前电网调度中心实际使用的状态估计器绝大多数仍然以SCADA量测为主PMU量测作为一种高精度补充逐步融入但还没有做到全覆盖。在这个背景下混合量测状态估计的结构也就是传统量测加PMU量测的WLS框架会长期存在。从维护角度看是否纳入PMU量测以及如何设置PMU权重往往决定了一套状态估计软件能不能在现场稳定运行。有的调度中心遇到过PMU数据跳变导致状态估计结果剧烈波动的案例最后定位下来问题出在PMU相角参考不稳而不是算法本身。处理方式是给PMU量测加上数据质量标签只有数据质量正常的PMU量测才参与估计有效解决了误动问题。5.2 后续还可以扩展的方向这个项目做完后如果还想深入学习有几个清晰的扩展方向可以做。第一个方向是抗差状态估计把WLS目标函数换成M估计的Huber函数或IGG权函数观察在不良数据比例较高时两种算法的性能对比。第二个方向是动态状态估计利用PMU的高速率测量数据结合卡尔曼滤波框架估计系统状态的动态轨迹。第三个方向是分布式状态估计适合解决大规模分区电网的在线状态感知问题。不管往哪个方向扩展Matlab的仿真环境都是绝佳的试验田。你可以快速改动算法核心对比不同策略的效果不用考虑现场硬件限制和数据通信问题这种自由度在工程现场很难获得。从我个人体会来说这个项目最值得投入精力的环节有三个一是逐项推导量测雅可比矩阵并对照数值差分验证这一步做扎实了状态估计就算入门了二是把NR潮流的雅可比矩阵和WLS的量测雅可比矩阵放在一起对比体会两种问题结构的差异三是做PMU配置方案对估计精度影响的敏感性分析它能帮你建立对量测系统设计的直觉。把这些环节啃下来后续无论接触什么状态估计算法你都不会再觉得那是黑盒。 基于WLS与PMU的电力系统状态估计：Matlab实现与对比分析
1. 项目概述与核心思路拆解1.1 这个项目到底在解决什么问题电力系统状态估计是EMS能量管理系统最核心的功能模块之一。调度中心每天面对海量遥测数据这些数据来自SCADA系统受CT/PT变比误差、通道噪声、测量装置老化等因素影响量测值并不精确。如果直接用这些有误差的数据去做在线分析、安全评估结果很可能失真。状态估计的作用就是利用一组带误差的冗余量测按照统计学准则估算出系统最可能处于的运行状态这个状态通常用所有节点的电压幅值和相角来表示。这个项目的标题很直白它要做的事情是搭一个可复现的IEEE标准测试系统用加权最小二乘法WLS配合PMU相量量测估算系统中各节点的电压幅值和相角同时用Newton-Raphson潮流计算求解同一系统的状态把PMU/WLS的估计结果和NR潮流的解放在一起对比检验状态估计算法的准确性。项目用Matlab实现适合电力系统方向的研究生、刚接触EMS算法的工程师以及研究PMU在电网中应用的同学参考。做完这个项目你不仅仅会运行一段代码你会把量测方程、雅可比矩阵、权重矩阵、可观测性这几块硬骨头一次啃透。之所以拿NR潮流做基准不是因为潮流比状态估计更高级而是因为潮流求解的问题本身就是确定的给定各节点注入功率方程组数量等于未知量数量解是确定的。状态估计处理的是量测冗余问题方程数多于未知数解是统计意义上的最优估计。如果把量测噪声控制得很小WLS的估计结果理论上会非常逼近NR潮流的解。所以用NR潮流解作为标准答案去验证状态估计在含噪声量测条件下的表现逻辑上说得通工程上也是常用的验证手段。1.2 PMU与传统量测的代差体现在哪里PMU即相量测量单元和传统RTU远程终端单元最本质的区别在于PMU能直接测量电压和电流的相量也就是幅值和相角。深入研究PMU在状态估计中的作用首先要理解传统量测的局限。SCADA系统里大量使用的是功率量测有功P、无功Q和电压幅值量测相角信息无法直接测到只能靠状态估计算法从功率方程里间接推算出来。而且传统RTU的扫描周期是秒级各个量测点之间没有统一的时间基准同步性差遇到系统动态过程时量测数据反映的其实不是同一个断面。PMU依靠GPS/北斗卫星授时采样率通常做到几十帧每秒所有PMU上报的数据都有统一时标天然是同步断面的。更关键的是PMU能直接给出电压相角这正好补齐了传统量测的短板。在这个项目里把PMU量测纳入WLS状态估计后相角不再完全依赖间接推算信息矩阵的形态会更好估计精度也会提升。不过需要注意PMU在工程上还没有做到全网全覆盖更现实的方案是传统量测为主、PMU量测为辅的混合量测架构本项目的量测配置也遵循这个思路。1.3 为什么WLS依然是状态估计的主流选择尽管学术界提出了不少新算法比如抗差估计M估计、最大指数绝对值估计等但WLS凭借模型简单、计算效率高、在正态噪声假设下具有最优线性无偏性质至今仍然是工程界应用最广的状态估计内核。WLS的核心思想并不复杂把所有量测方程写成一个超定方程组然后寻找一组状态变量使量测残差的加权平方和最小。权重取什么值很关键通常取量测噪声方差的倒数也就是说精度越高的量测在目标函数里的话语权越大。具体到这个项目WLS要处理三种量测支路功率线路首末端P和Q、节点注入功率P_i和Q_i、电压幅值量测再加上PMU提供的节点电压相角量测。结果就是在原有量测方程组里增加一组线性方程相角量测方程简单雅的可比矩阵那一行也简单但它在算法收敛性和解算精度上带来的提升实测下来非常明显。整个迭代求解过程本质上是Gauss-Newton法每轮迭代利用线性化的量测方程构造法方程解出状态修正量后更新直到修正量的最大绝对值小于收敛阈值。2. 状态估计与潮流计算的数学对照2.1 WLS状态估计的数学框架状态估计的数学模型普遍采用如下量测方程描述z h(x) e其中z是量测向量x是系统状态向量通常取所有节点的电压幅值和相角参考节点相角除外h(x)是量测函数表示状态到量测的非线性映射e是量测噪声向量一般假设为零均值高斯分布协方差矩阵为R。量测噪声的方差由仪表精度、变送器误差、通信量化误差等因素共同决定。目标函数定义为J(x) [z - h(x)]^T R^{-1} [z - h(x)]对J(x)求梯度并令其为零得到迭代公式Δx (H^T R^{-1} H)^{-1} H^T R^{-1} [z - h(x)]其中H是量测雅可比矩阵第i行第j列元素表示第i个量测函数对第j个状态变量的偏导数。这个式子看着吓人实际上每一步迭代都是在解一个线性最小二乘问题G H^T R^{-1} H称为信息矩阵或增益矩阵工程上也直接叫它G矩阵。当量测冗余度足够、系统可观测时G矩阵是对称非奇异的可以稳定求逆。说白了WLS做的就是在量测残差的加权平方和最小的意义下把状态量反推出来。如果量测全是PMU提供的电压相量那么h(x)就是线性的一步就能得到解。但实际系统中功率量测占据大头所以还是得走迭代。2.2 Newton-Raphson潮流计算的求解框架潮流计算的问题描述要简洁得多给定平衡节点电压幅值与相角、PV节点有功注入和电压幅值、PQ节点有功和无功注入求所有节点的电压幅值和相角。核心方程组是节点功率平衡方程P_i V_i Σ_{j} V_j (G_ij cosθ_ij B_ij sinθ_ij) Q_i V_i Σ_{j} V_j (G_ij sinθ_ij - B_ij cosθ_ij)其中G_ij和B_ij是节点导纳矩阵元素的实部和虚部θ_ij是节点i和j的相角差。NR法把这个非线性方程组在初值附近做泰勒展开忽略二阶以上项得到线性修正方程组[ΔP; ΔQ] J_NR [Δθ; ΔV]J_NR是潮流雅可比矩阵由P、Q对θ、V的偏导数构成。每轮迭代计算当前状态下的注入功率与给定注入做差得到失配量然后解线性方程组求修正量更新状态后继续迭代直到失配量足够小。NR潮流的收敛特性好正常条件下平启动也能在几次迭代内收敛所以在工程软件里基本是默认选择。这个项目里NR潮流还有一个特殊用途就是产生标准答案先用潮流解出系统真实状态再在这个真实状态上叠加高斯噪声模拟量测值这样我们就知道了量测对应的真值评估WLS估计误差才有依据。2.3 两个方法到底有什么本质区别很多初学者会把潮流计算和状态估计搞混因为两个算法都在解电压幅值和相角公式长得像迭代框架也像。但它们的本质区别非常明显。从输入数据看潮流计算用的是发电和负荷的计划值或预测值是一种模型驱动状态估计用的是实际量测是数据驱动。从方程数量看潮流是方阵求解方程数和未知数数量一致状态估计是超定求解量测方程数量远大于状态变量数量。从解的性质看潮流解是模型输入下的确定性响应状态估计解是统计意义下的最优估计它能同时给出量测残差供坏数据检测使用。这个对照反映在代码层面就是NR潮流失配量计算直接对应功率平衡方程没有权重概念WLS的右端项则是加权残差H^T W (z - h(x))包含权重矩阵。两者的雅可比矩阵也不一样潮流的雅可比是功率对电压幅值和相角的偏导状态估计的雅可比是各类型量测函数对状态的偏导维度通常比潮流大得多。在项目实现时我会把两者放在同一个脚本框架里量测雅可比和潮流雅可比分开构建避免搞混。3. Matlab实现流程与关键环节解析3.1 测试系统选型与数据准备做状态估计不能上来就在自己拼的随机网络上跑否则哪里算错都说不清。靠谱的做法是直接用IEEE标准节点系统推荐IEEE 14节点系统。它的规模适中既不会因为网络太大导致可视化困难又能体现出状态估计对全网状态的感知能力。IEEE 14节点系统包含5个发电节点、11条负荷支路、3台变压器和1条并联电容器支路量测配置可以做得比较丰富。在Matlab中我建议直接借助Matpower工具箱加载数据它可以省去手动构造线路参数和节点注入数据的麻烦mpc loadcase(case14); Ybus makeYbus(mpc); baseMVA mpc.baseMVA; bus mpc.bus; branch mpc.branch; gen mpc.gen;这里loadcase会返回一个结构体里面包含节点数据矩阵、支路数据矩阵、发电机数据矩阵。节点数据矩阵的各列含义需要弄清楚其中第8列是电压幅值标幺值第9列是电压相角角度制后续构造量测时会用到。makeYbus则根据节点和支路数据生成全系统的节点导纳矩阵这是所有后续计算的基础。有一点要强调Matpower里的负荷和发电机数据都是标幺值量测构造时要保持同样的标幺体系。如果直接拿有名值去算导纳矩阵的量级会对不上迭代计算容易出现数值发散。做项目时先在脚本开头统一单位的处理后续所有量测和状态量都用标幺值或弧度制只在输出结果时换算成有名值或角度制。3.2 量测数据生成与PMU配置方案WLS状态估计需要量测数据作为输入。在真实系统中量测来自RTU和PMU在仿真项目里需要人为构造量测。我的做法是分三步。第一步用NR潮流计算系统的精确状态得到所有节点电压幅值V_true和相角θ_true弧度制。这个状态就是后续评估的基准。第二步根据要量测的物理量和拓扑关系计算理想的量测值比如支路潮流P_ij通过功率方程由V_true和θ_true算出。第三步在理想量测值上叠加高斯噪声模拟真实量测。Matlab里生成带噪声量测的代码如下% 设定量测误差标准差标幺值 sigma_P 0.01; % 有功量测误差标准差约为1% sigma_Q 0.01; % 无功量测误差标准差 sigma_V 0.005; % 电压幅值量测误差标准差0.5% sigma_theta 0.001; % PMU相角量测误差标准差约0.057度 % 假设已通过潮流计算出P_true, Q_true, V_true, theta_true z_P P_true sigma_P * randn(size(P_true)); z_Q Q_true sigma_Q * randn(size(Q_true)); z_V V_true sigma_V * randn(size(V_true)); z_theta theta_true sigma_theta * randn(size(theta_true_pmu));PMU配置方案我推荐这样设计不是所有节点都装PMU那不符合工程实际。典型做法是在部分关键节点配置PMU比如发电机出口节点和电压中枢节点。在每个配置了PMU的节点能量测到的量包括该节点的电压幅值和相角以及与该节点相连的支路电流幅值相角。为了项目简单可以只把PMU的电压相角量测加入量测集也就是在传统量测基础上增加一组相角量测方程。量测冗余度这个概念需要重点理解。系统状态变量的数量是2n - 1n个节点每个节点有幅值和相角但参考节点相角固定量测总数必须显著大于状态数才能保证可观测性并具备抗差能力。一般工程实践中量测冗余度在2到3以上比较好。IEEE 14节点系统状态变量共27个我建议构造至少50到60个量测既包含支路潮流、节点注入、电压幅值再叠加PMU相角量测这样信息矩阵条件数才好看。3.3 WLS状态估计迭代实现核心代码WLS的核心迭代循环并不复杂但有几个细节必须处理到位否则程序很容易跑飞或不收敛。下面给出一个可以直接改用的框架function [V_est, theta_est, iter] wls_estimator(Ybus, z, W, idx, loadflow_solved, tol, max_iter) % Ybus: 节点导纳矩阵 % z: 量测向量按固定顺序组织支路有功、支路无功、注入有功、注入无功、电压幅值、PMU相角 % W: 权重对角矩阵 % idx: 结构体记录各类量测在z中的索引位置 % 初始化平启动 n size(Ybus, 1); V ones(n, 1); theta zeros(n, 1); x [theta(2:end); V]; % 状态向量参考节点相角固定为0 for iter 1:max_iter % 构建量测雅可比矩阵 H 和量测函数值 h(x) [H, h] build_measurement_jacobian(V, theta, Ybus, idx); % 增益矩阵 G H * W * H; % 右端项 rhs H * W * (z - h); % 求解修正量 dx G \ rhs; % 更新状态 dtheta [0; dx(1:n-1)]; dV dx(n:end); theta theta dtheta; V V dV; if max(abs(dx)) tol break; end end theta_est theta; V_est V; endbuild_measurement_jacobian函数是整个程序里最容易写错的地方。它要对每种量测类型分别求偏导。以支路有功量测为例P_ij对θ_i、θ_j、V_i、V_j都有偏导公式不复杂但容易漏项。一个实用的检查办法是用有限差分法在相同状态下数值求雅可比和解析表达式做对比最大误差在1e-4以下基本可以确认偏导公式写对了。信息矩阵G H^T W H的维度是(2n-1)×(2n-1)在IEEE 14节点系统上规模很小直接用Matlab的左除运算符\即可求解。如果系统规模很大则需要利用G矩阵的稀疏性进行因子分解避免对稠密矩阵求逆。IEEE 14节点系统上直接稠密求解没有性能压力。还有一个细节PMU量测的相角是绝对相角而状态估计的参考节点相角固定为0。如果PMU配置在参考节点之外的节点它的相角量测方程是θ_i加上一个系统参考偏移在算法实现上要么把参考节点同PMU对齐要么在量测方程中把参考相角作为一个附加状态变量来处理。本项目里建议直接选择参考节点也配置PMU这样处理最简洁避免了相角参考不一致的问题。3.4 Newton-Raphson潮流求解实现要点NR潮流在这个项目里的定位是真值生成器同时也作为状态估计的对照基准。用Matpower的runpf可以直接得到潮流结果但如果想深入理解算法还是建议自己写一遍核心迭代。自写NR潮流有几个关键点。第一节点分类要清晰ref find(bus(:, 2) 3); % 平衡节点 pv find(bus(:, 2) 2); % PV节点 pq find(bus(:, 2) 1); % PQ节点第二雅可比矩阵的构建。NR潮流的雅可比矩阵分成4个子块dP/dθ、dP/dV、dQ/dθ、dQ/dV。注意PV节点的无功方程不参与迭代因为PV节点的无功注入是待定量平衡节点的有功和无功方程都不参与迭代。最终雅可比矩阵的维度是(2×PQ数量 PV数量)的方阵。第三迭代修正量的更新方式。电压幅值采用相乘修正V_new V ΔV相角采用相加修正θ_new θ Δθ。收敛判据推荐使用功率失配量的最大绝对值通常设为1e-8。自己写一遍NR潮流会让你对功率平衡方程的理解从名词变成肌肉记忆后续构建量测雅可比时很多偏导数公式其实是互通的会省不少力气。3.5 结果对比与误差指标体系完成WLS状态估计和NR潮流求解后就来到了验收环节。需要对比的物理量主要是两类电压幅值V的估计值和电压相角θ的估计值。通常的控制台输出方式是把所有节点的数据按三列并排打出来真值、WLS估计值、NR潮流计算值。我习惯用如下方式组织结果err_V abs(V_est - V_true); err_theta abs(theta_est - theta_true); mean_err_V mean(err_V); max_err_V max(err_V); rmse_V sqrt(mean((V_est - V_true).^2));从实际运行结果来看当量测误差标准差设置为典型值有功1%、无功1%、电压幅值0.5%、PMU相角0.001弧度时WLS估计的电压幅值平均误差通常能控制在0.1%以内相角平均误差大约在0.01到0.05度之间。如果PMU数量增加相角估计精度会有更明显的提升。图形对比建议画两个图。第一个图画出各节点的电压幅值真值、WLS估计值和NR潮流值的三条曲线能直观看到WLS估计与真值之间的偏差大小。第二个图画出相角误差的柱状图从图上能清晰看出哪些节点相角估计误差偏大通常是没有PMU覆盖的负荷节点误差更大这是符合直觉的。4. 调试心得、收敛问题与工程避坑指南4.1 初值选择与收敛性判断WLS状态估计的初值一般取平启动也就是所有节点电压幅值取1.0相角取0。对绝大多数正常量测配置来说平启动就能收敛。但有两种情况容易出问题一是量测冗余度不够G矩阵接近奇异迭代会出现振荡不收敛二是量测中存在不良数据导致残差过大迭代方向脱离可行域。我在调试时遇到过一种典型现象程序在头几次迭代里修正量很大看起来要发散但其实是因为初值状态距离真实解太远线性化模型误差较大。处理方法很简单不要急着判定发散观察一下修正量的曲线趋势。如果修正量虽然在波动但整体在减小就让它继续迭代如果修正量持续增大或者出现NaN再考虑调整初值或检查雅可比矩阵。经验做法是先用NR潮流算一遍把潮流解作为WLS的初值。这样做的好处是初值距离真值足够近WLS迭代次数会明显减少而且能有效避免因为初值太差导致的数值问题。代价是潮流本身要先解出来但在这个项目里潮流本来就要算所以白赚一个稳定初值。4.2 权重矩阵W的设置到底有什么讲究权重矩阵W是量测噪声协方差矩阵R的逆矩阵。理论上讲W的取值应该反映量测的真实精度。但工程上量测精度并不总是准确已知权重的设置往往成了影响估计质量的关键因素。从实践来看有几条经验可以直接用。各类量测误差标准差的典型参考值大致如下量测类型典型误差标准差权重参考值支路有功P_ij0.011%10000支路无功Q_ij0.011%10000节点注入有功P_i0.011%10000节点注入无功Q_i0.011%10000电压幅值V_i0.0050.5%40000PMU相角θ_i0.0010.057°1000000注意PMU相角的权重比传统量测高出两个数量级这在数学上是合理的因为PMU相角的测量精度确实远高于由功率量测间接推算的相角。但一个隐患是如果PMU量测本身存在系统误差比如相角基准偏差过高的权重反而会把估计结果往错误方向带。所以工程上对PMU数据一般会先做合理性校验再进入状态估计器。权重矩阵的病态问题也值得留意。当权重数值跨越5到6个数量级时增益矩阵G的条件数会变得很大求解线性方程组时数值误差可能掩盖真实解。遇到这种情况可以对权重做适度缩放比如统一除以最大权重把权重范围压到10^(-3)到1之间数值稳定性会好很多。4.3 量测冗余度对估计精度的影响量测冗余度的定义是量测数量与状态数量之比。这个指标直接影响估计精度和抗差能力。我做过一组测试把量测从最少的刚好可观测冗余度接近1逐步增加到冗余度超过3对比WLS的估计误差。结果非常有意思冗余度从1提升到2估计误差下降非常显著差不多能下降50%再从2提升到3误差下降幅度变缓大概再降20%左右。这背后的道理也不难理解。冗余量测提供了更多关于系统状态的约束信息最小二乘估计本质上是多个量测信息的折中量测越多随机噪声被平滑掉的效果越好。但冗余度超过3之后新增量测带来的边际收益越来越小而量测系统的建设和维护成本却直线上升。工程上设计量测配置时要找到性价比平衡点不是量测越多越好。在IEEE 14节点系统上做实验时一个很直观的做法是先只布置最小可观测量测比如各节点电压幅值平衡节点和若干支路功率跑一遍WLS记录误差然后逐步添加支路功率量测、注入量测、PMU相角量测每次添加后重新运行记录误差变化。你会看到一个清晰曲线这也是理解量测冗余度最好的实验。4.4 不良数据检测的处理思路实际系统的量测数据难免出现坏数据比如采集模块故障发出一个特别离谱的值或者通信传输偶发丢包导致数据错乱。如果不做处理一个坏数据就足以让WLS估计结果整体偏移污染整个状态估计质量。这就是工程上为什么状态估计必须配不良数据检测模块。最常见的方法是残差分析法。WLS求解收敛后可以计算每个量测的归一化残差r_i (z_i - h_i(x_hat)) / sqrt(R_ii)在正态噪声假设下如果量测正常r_i应当服从标准正态分布大概率落在[-3, 3]区间内。如果某个量测的归一化残差显著超过3基本可以判定它是坏数据。工程实现的时候比较耗时的是剔除坏数据后的重新估计。一种策略是先把所有超限的量测全部清除然后重新跑WLS再检查残差另一种策略是每次只剔除残差最大的那个量测重跑后再检查迭代执行。后者的计算量更大但更稳健因为同时剔除多个量测可能误杀正常数据。这个项目里可以在量测集里刻意注入一两个坏数据比如把某个P_ij量测增大50%然后观察不良数据检测机制能否准确把它揪出来。这个实验做完你对状态估计的工程细节会有一个质的理解提升。5. 从仿真到工程现场的延伸思考5.1 混合量测状态估计的工程意义前面所有工作都是在一个仿真环境里完成的。但仿真项目的价值不只是交一份报告更重要的是理解工程实现背后的逻辑。当前电网调度中心实际使用的状态估计器绝大多数仍然以SCADA量测为主PMU量测作为一种高精度补充逐步融入但还没有做到全覆盖。在这个背景下混合量测状态估计的结构也就是传统量测加PMU量测的WLS框架会长期存在。从维护角度看是否纳入PMU量测以及如何设置PMU权重往往决定了一套状态估计软件能不能在现场稳定运行。有的调度中心遇到过PMU数据跳变导致状态估计结果剧烈波动的案例最后定位下来问题出在PMU相角参考不稳而不是算法本身。处理方式是给PMU量测加上数据质量标签只有数据质量正常的PMU量测才参与估计有效解决了误动问题。5.2 后续还可以扩展的方向这个项目做完后如果还想深入学习有几个清晰的扩展方向可以做。第一个方向是抗差状态估计把WLS目标函数换成M估计的Huber函数或IGG权函数观察在不良数据比例较高时两种算法的性能对比。第二个方向是动态状态估计利用PMU的高速率测量数据结合卡尔曼滤波框架估计系统状态的动态轨迹。第三个方向是分布式状态估计适合解决大规模分区电网的在线状态感知问题。不管往哪个方向扩展Matlab的仿真环境都是绝佳的试验田。你可以快速改动算法核心对比不同策略的效果不用考虑现场硬件限制和数据通信问题这种自由度在工程现场很难获得。从我个人体会来说这个项目最值得投入精力的环节有三个一是逐项推导量测雅可比矩阵并对照数值差分验证这一步做扎实了状态估计就算入门了二是把NR潮流的雅可比矩阵和WLS的量测雅可比矩阵放在一起对比体会两种问题结构的差异三是做PMU配置方案对估计精度影响的敏感性分析它能帮你建立对量测系统设计的直觉。把这些环节啃下来后续无论接触什么状态估计算法你都不会再觉得那是黑盒。