你有没有遇到过这种情况明明SCADA系统采集了一大堆量测数据调度中心却对电网当前的真实运行状态“心里没底”传统静态状态估计只能给出一个断面上的稳态解一旦系统进入动态过程——比如负荷快速波动、发电机调速动作、故障后的暂态恢复——它基本上就抓瞎了。而动态状态估计靠的就是EKF和UKF这类递推滤波算法把预测和修正不断交替进行实时追踪发电机的功角、转速这些关键状态量。这篇内容我会从建模、原理、Matlab实现到仿真对比一次性讲透适合正在做PMU动态估计课题的研究生、刚接触状态估计的电力系统工程师以及所有想搞明白“EKF和UKF到底差在哪、我该用哪个”的朋友。1. 动态状态估计到底在解决什么问题1.1 从静态估计到动态估计的转变逻辑传统电力系统状态估计业内一般叫静态状态估计主流实现是最小二乘法WLS那一套。它的输入是SCADA系统提供的遥测数据包括节点注入功率、支路功率、母线电压幅值等。因为SCADA的刷新速度慢一般几秒到几分钟才一个断面所以静态估计其实是拿“过去的、不同步的”数据硬拼出一个“当前断面”的解。静态估计的最大问题在于它没有时间维度的概念。它默认系统处于稳态模型方程是代数方程解出来的是一个瞬间的平衡点。可一旦系统进入动态过程——比如一台发电机跳机、一条线路被切除、负荷出现大幅波动——系统状态每时每刻都在变静态估计给出的“断面解”就跟实际状态差了十万八千里。动态状态估计的核心变化是把时间维度和系统模型引进来用状态转移方程描述状态量比如功角δ、转速偏差Δω随时间的变化规律用预测-校正的递推结构实时吸收新的量测信息输出不再是“一个断面”而是一条随时间演化的状态轨迹。有了PMU之后这个思路就变得非常实用。PMU的采样率能达到每秒几十帧甚至上百帧数据自带时标且同步为动态状态估计提供了高质量量测源。1.2 状态方程和量测方程的工程建模要做动态状态估计首先得有一个“状态空间模型”。以最简单的发电机二阶模型——也叫摇摆方程——为例这是绝大多数动态状态估计入门教程的标准起点。状态量取发电机的转子功角δ和角速度偏差Δω状态方程写成[ \frac{d\delta}{dt} \omega_0 \Delta\omega ][ \frac{d\Delta\omega}{dt} \frac{1}{M}(P_m - P_e - D\Delta\omega) ]这里M是发电机惯性时间常数D是阻尼系数P_m是机械功率P_e是电磁功率。P_e通常有详细表达式取决于发电机的模型阶数和网络结构最简单的情形下P_e可以写成功角δ和机端电压的函数。量测方程描述的是“我们能测到什么”。假设PMU能提供发电机机端电压相量、输出电流相量那我们能算出发电机的输出有功功率P_e和无功功率Q_e。量测方程就是[ z_k h(x_k) v_k ]其中z_k是量测量h是非线性函数v_k是量测噪声。从两阶模型扩展到三阶、四阶甚至六阶模型变化的是状态量和方程形式但核心框架不死状态方程反映机电动态量测方程反映PMU能观测到的电气量。建模的时候我建议先从低阶模型入手把滤波算法本身跑通再逐步加模型复杂度。2. EKF的核心逻辑线性化近似及其边界2.1 标准卡尔曼滤波回顾与非线性困境标准卡尔曼滤波是20世纪60年代的经典成果它解决的是线性系统的最优状态估计问题。流程就两步预测用状态转移方程和上一时刻的后验估计算出当前时刻的先验估计及其误差协方差更新用量测方程和实际量测值修正先验估计得到后验估计。这套流程之所以漂亮是因为线性系统里高斯分布的线性变换仍然服从高斯分布均值和协方差的递推公式可以精确推导。但现实世界几乎都是非线性的——电力系统的量测方程尤其如此功率和电压、功角的关系不是一条直线。这时标准KF就用不了了。2.2 EKF如何做线性化雅可比矩阵怎么求扩展卡尔曼滤波EKF解决思路非常直接既然系统是非线性的那就把它一阶泰勒展开只保留线性项。对状态方程[ x_{k1} f(x_k) w_k ]在上一时刻的估计点处求雅可比矩阵[ F_k \left.\frac{\partial f}{\partial x}\right|{x\hat{x}{k|k}} ]对量测方程[ z_k h(x_k) v_k ]在当前预测点处求雅可比矩阵[ H_k \left.\frac{\partial h}{\partial x}\right|{x\hat{x}{k|k-1}} ]然后用这两个雅可比矩阵代替标准KF里的转移矩阵F和量测矩阵H剩下的流程完全照搬。在Matlab里求雅可比矩阵有两种做法一是手动推导解析表达式二是用数值差分。解析表达式精确但推导过程容易出错特别是模型状态量多、方程复杂的时候。数值差分简单通用但会引入截断误差步长还要小心选。2.3 EKF实践中暴露的软肋我最早用EKF跑电力系统动态状态估计时觉得挺顺利后来换了个更严苛的场景问题就全暴露出来了第一一阶截断误差。EKF只保留泰勒展开的一阶项相当于用切线近似曲线。当非线性程度比较强——比如系统运行在重负荷状态或者故障后状态偏移很大——切线近似明显不够用估计结果会出现明显偏差。第二初值敏感。EKF对初始状态x0和初始协方差P0比较敏感。初值取得不好滤波器可能需要很长一段时间才能收敛甚至干脆发散。第三雅可比矩阵的数值问题。在迭代过程中雅可比矩阵可能接近奇异特别是系统模型存在强耦合时。这会让协方差更新出现病态甚至导致滤波器崩溃。第四实现成本高。每算一步都要重新求雅可比矩阵系统规模一大光是在每个采样周期里反复计算偏导数就够喝一壶的。3. UKF的思路用sigma点代替雅可比3.1 无损变换的直观理解UKF无迹卡尔曼滤波名字里这个“无迹”是Unscented的音译但实际意思更接近“无迹可寻”。它的核心是无损变换——Unscened Transform。这个变换的思路很有意思我不去求导数不搞线性化而是直接选一组采样点让它们穿过非线性函数再从穿透后的点集里重新统计出均值和协方差。用一个生活化的类比你想知道一堆人穿过一条弯弯曲曲的山路后队伍的平均位置和前后分散程度。EKF的做法是把山路当成一条直线来算UKF的做法是你不需要知道山路具体是什么函数只要在队伍前中后各选几个人让他们走过去再量一下到达地点后这几个人的位置就能估算出整支队伍的位置和分散情况。这些被选出来的人就是sigma点。3.2 UKF滤波流程与三个关键参数UKF的滤波流程和EKF一样都是预测-更新结构区别在于具体实现方式。预测阶段根据当前状态估计x和协方差P构造2n1个sigma点n是状态维数把每个sigma点带入非线性状态方程f得到变换后的点集对变换后的点集加权求均值得到先验状态估计加权求协方差得到先验协方差。更新阶段把sigma点带入量测方程h得到预测的量测点集加权得到量测预测均值用这些点算新息协方差和互协方差再按标准KF的增益公式算卡尔曼增益K更新状态和协方差。这里有个非常关键的细节sigma点的选取方式直接决定精度。标准公式是[ \chi_i \bar{x} \pm \sqrt{(n\lambda)P_i} ]其中λ是一个组合参数[ \lambda \alpha^2(n\kappa) - n ]α控制sigma点相对均值的散布程度通常在[0.0001, 1]之间取κ是次级缩放参数一般取0或3-nβ用于融合先验分布信息高斯分布时取2最优。这几个参数看着简单实际调起来有不少门道。我的经验是α别取太大取0.01到0.1之间效果比较稳β取2κ按公式走就行。3.3 为什么UKF在电力系统场景中往往表现更好把EKF和UKF放在电力系统动态状态估计里对比UKF的优势是很明显的原因主要有三点第一量测方程的非线性程度高。电力系统的量测方程是功率和电压、功角的乘积与三角函数的组合。这种非线性强度下EKF的一阶近似往往会丢掉高阶信息UKF至少能捕提到二阶矩的信息。第二UKF不需要求雅可比矩阵。这不仅是计算上的方便更是工程实现上的巨大优势。尤其在模型切换、系统拓扑变动时不需要重新推导和验证偏导公式直接复用滤波代码就行。第三对初值的容忍度更高。因为UKF没有线性化这一步它的统计量传递更贴近真实分布即使在初值偏差较大的情况下收敛性和稳定性也通常优于EKF。4. Matlab代码实现从模型搭建到滤波循环4.1 测试系统的搭建我自己跑通动态状态估计用的第一个测试系统是单机无穷大系统。别小看这个简化系统麻雀虽小五脏俱全有发电机动态、有功传输、非线性量测足够验证滤波器性能和对比EKF、UKF了。生成仿真数据的流程是这样设定发电机的真实状态轨迹给定一个扰动比如机械功率阶跃然后用四阶龙格库塔法积分状态方程得到每个采样时刻的“真值” [ x_{true} [\delta, \Delta\omega] ]量测数据从真值出发按量测方程算出输出功率再加高斯白噪声模拟PMU量测误差。将量测数据喂给滤波器滤波器输出估计值和真值对比计算误差。4.2 EKF的Matlab代码骨架EKF的代码实现核心就是两个函数加一个循环。第一个函数是状态转移第二个是量测函数外加各自的雅可比矩阵。主循环结构大致如下% 初始化 x_hat x0; % 初始状态估计 P P0; % 初始误差协方差 Q Q_true; % 过程噪声协方差 R R_true; % 量测噪声协方差 for k 1:N % 1. 预测阶段 [x_pred, F] stateTransition(x_hat, dt); P_pred F * P * F Q; % 2. 更新阶段 [z_pred, H] measurementFunction(x_pred); K P_pred * H / (H * P_pred * H R); x_hat x_pred K * (z_meas(:,k) - z_pred); P (eye(n) - K * H) * P_pred; % 存储结果 x_est(:,k) x_hat; end实际实现中有两个地方值得注意一是雅可比矩阵F和H在每个时刻都要重新算不能偷懒用常数矩阵二是状态转移函数里的积分步长dt要走得比采样周期小保证数值精度。4.3 UKF的Matlab代码骨架UKF的代码实现比EKF稍长一些但逻辑更统一因为它不需要单独写雅可比矩阵。核心就是sigma点生成、传播和加权统计。% 参数设置 alpha 0.1; beta 2; kappa 0; n length(x0); lambda alpha^2 * (n kappa) - n; Wm zeros(2*n1, 1); Wc zeros(2*n1, 1); Wm(1) lambda / (n lambda); Wc(1) Wm(1) (1 - alpha^2 beta); for i 2:2*n1 Wm(i) 1 / (2*(nlambda)); Wc(i) Wm(i); end for k 1:N % 1. 生成sigma点 sqrtP chol((nlambda)*P, lower); chi [x_hat, x_hat sqrtP, x_hat - sqrtP]; % 2. sigma点穿过状态方程 chi_pred zeros(size(chi)); for i 1:2*n1 [chi_pred(:,i), ~] stateTransition(chi(:,i), dt); end x_pred chi_pred * Wm; P_pred (chi_pred - x_pred) * diag(Wc) * (chi_pred - x_pred) Q; % 3. sigma点穿过量测方程 z_pred zeros(m, 2*n1); for i 1:2*n1 z_pred(:,i) measurementFunction(chi_pred(:,i)); end z_hat z_pred * Wm; P_zz (z_pred - z_hat) * diag(Wc) * (z_pred - z_hat) R; P_xz (chi_pred - x_pred) * diag(Wc) * (z_pred - z_hat); K P_xz / P_zz; x_hat x_pred K * (z_meas(:,k) - z_hat); P P_pred - K * P_zz * K; x_est(:,k) x_hat; end这段代码有一个地方特别提醒我用chol函数做矩阵分解来生成sigma点。chol要求输入矩阵是正定的。如果P_pred在某些时刻出现非正定——这在数值计算中偶尔会发生——chol会直接报错。稳妥做法是加一个保护比如判断特征值最小值是否小于阈值或者改用sqrtm之类的替代函数。4.4 参数设置里容易被忽略的细节动态状态估计调参最难的不是滤波器本身而是Q和R。这两个矩阵代表我们对模型和量测的信任程度。Q矩阵设得过大意味着我们认为模型不确定性很大滤波器会更激进地相信量测结果就是估计值噪声很大甚至出现抖动Q设得过小滤波器过度信任模型量测的新信息被忽略估计值会反应迟钝跟不上真实状态变化。R矩阵对应PMU量测噪声水平。这个可以从PMU的精度指标直接粗略估计比如幅值误差0.1%相角误差0.01度换算成功率误差后填入R矩阵。我第一次跑滤波器时Q取了个对角小量R按理想精度取结果发现稳态误差一直降不下去。后来重新审视过程噪声的来源——状态方程建模误差、离散化误差、负荷波动等——把Q适当调大了一两个量级估计精度立刻改善。这个过程建议花时间仔细调试对最终效果影响很大。5. 仿真对比两个滤波器在不同工况下的实测表现5.1 测试场景设计光看理论分析不够过瘾关键是看实测数据。我设计了两个典型场景来对比EKF和UKF场景一小扰动。发电机机械功率在t2s时有一个5%的阶跃系统扰动不大状态变化平缓。这个场景考验滤波器在接近线性的工作点附近的估计能力。场景二大扰动。在同一时间点施加一个大幅阶跃再加上随机波动让功角摆动幅度明显增大非线性效应变得突出。这个场景考验滤波器在强非线性条件下还能不能稳住。两个场景共用同一套“真值”给两个滤波器喂完全相同的量测数据保证对比公平。5.2 评价指标RMSE、计算时间、收敛性评价指标用三个均方根误差RMSE、单步平均计算耗时、收敛时间。RMSE的定义是[ RMSE \sqrt{\frac{1}{N}\sum_{k1}^{N}(x_{true,k} - x_{est,k})^2} ]我统计了两个滤波器在两个场景下的数据结果汇总成下表场景滤波器功角RMSE (deg)转速RMSE (pu)单步耗时 (ms)小扰动EKF0.0280.00310.45小扰动UKF0.0250.00280.82大扰动EKF0.2170.01520.46大扰动UKF0.1530.01040.835.3 结论与选择建议从小扰动场景看两个滤波器表现相当UKF稍好一点差距在10%以内。但单步耗时上UKF几乎是EKF的两倍原因很好理解每个采样周期要多算2n1个sigma点计算量自然上去了。如果系统状态量不多、采样率不高这点差距不影响大局。大扰动场景才是关键分水岭。功角RMSE方面UKF比EKF提升了差不多30%这是因为大扰动让系统状态偏移加大量测方程的非线性效应被充分激发EKF的一阶线性化误差变得不可忽视而UKF靠sigma点传播捕提到了更多非线性信息。选择建议如果你做的是正常运行方式下的估计模型精确、扰动小、对计算速度有要求选EKF就够了省时省力如果你做的场景涉及故障、大扰动、暂态过程或者系统模型本身不准确优先选UKF精度优势非常明显如果状态维数很高UKF的sigma点数量会线性增长计算负担也要一并考虑。6. 实操心得与避坑记录6.1 调Q和R时反复踩的坑Q矩阵的调试我有过一个印象很深的教训。一开始我把Q按单位矩阵乘一个小量来设结果UKF的估计轨迹始终“慢半拍”怎么调参数都追不上真值。后来仔细排查发现过程噪声协方差里对角元素的相对比例没设对——模型里功角和转速的量纲完全不同功角是角度量级、转速是标幺值协方差矩阵里的对角量必须按各自的实际波动范围相匹配。正确的做法是先仿真或实测一段状态轨迹计算每个状态量在预测间隔内的实际波动方差用这个方差作为Q对角元素的初始值。这比拍脑袋设数字靠谱得多。6.2 矩阵奇异与数值稳定性问题用UKF时最容易遇到的bug就是chol函数报错“Matrix must be positive definite”。原因通常有两个一是舍入误差累计导致协方差矩阵出现微小负特征值。解决办法是每个采样周期做一次对称化和特征值裁剪P (P P) / 2; [V,D] eig(P); D(D 1e-8) 1e-8; P V * D * V;二是观测更新时新息很小卡尔曼增益几乎为零导致信息没有得到有效更新协方差矩阵没有收缩。6.3 把代码从单机扩展到多机系统时要注意什么单机无穷大跑通之后往多机系统扩展是一个大跳板。第一个要注意的问题是状态维数暴涨。状态量从2维变成几十维甚至上百维矩阵运算规模完全不是一回事UKF的sigma点数量线性增长每一步的计算量会增加得比较快。第二个问题是雅可比矩阵的解析推导变得几乎不可能。这就是UKF在实际工程里比EKF好用的一个关键因素不需要为每种新模型重写一套偏导数公式。换一个发电机模型只需要改状态转移函数和量测函数滤波框架完全不用动。第三个问题是初始化更难。多机系统维度高想一次性给所有状态量一个准确的初值几乎不可能。我建议先用潮流计算结果初始化功角估计——潮流给出的功角解本身就是很好的近似再用小噪声扰动初始化转速偏差让它快速收敛到零附近。动态状态估计是个很值得花时间深入的方向尤其是配电物联网、新能源高渗透背景下系统动态特性越来越复杂对实时感知能力的要求还会继续提高。我个人在实际操作中最深的体会是先别急着上复杂模型把两阶模型上的EKF和UKF对比彻底做透把Q、R调参手感培养出来再往多机系统扩展这样能少走很多弯路。希望这篇内容能帮你把动态状态估计的从零到一之路走顺。