聊到电力系统动态状态估计很多跑过工程仿真的朋友第一反应还是传统的加权最小二乘WLS静态估计。但只要你做过PMU量测接入或者在线动态监控方向很快会意识到静态估计只给你一个“时间切片”没法描述系统从一个断面到下一个断面的演化过程。EKF和UKF这两类基于卡尔曼滤波框架的动态估计算法正是为了补上“时间维”这一块。这篇文章就从工程落地角度把扩展卡尔曼滤波EKF和无迹卡尔曼滤波UKF在电力系统动态状态估计中的Matlab实现思路、关键细节、参数选择以及我实际踩过的坑一条一条说清楚代码直接可参考适合正在做毕业设计、侧重算法仿真的研究生也适合刚接触状态估计、想把滤波用起来的工程师。1. 先搞清楚为什么需要动态状态估计1.1 静态估计的局限在哪里传统WLS静态状态估计处理的是某一个时刻断面的量测模型里没有时间递推关系本质上是求解一个加权最小二乘优化问题。你可以把它理解成“拿一堆照片拼出一张当前画面”每张照片的曝光时间不同但处理后只输出一张静态图。这在稳态运行下没什么问题电网调度中心本来就按秒级周期做断面状态估计。但问题出在动态过程中。当系统发生扰动、负荷波动或者机组故障时状态量变化很快。WLS只利用当前断面量测没法融合系统动态方程中蕴含的演化规律也没有办法对下一个时刻的状态做预测。更致命的是WLS对量测冗余度要求高如果某些量测通道出现延迟或者缺失估计结果就会抖动得非常厉害。这也是为什么很多研究开始转向“动态状态估计”也就是DSFDynamic State Filtering思路从“每步独立求解”变成“预测-更新闭环递推”。1.2 卡尔曼滤波框架为什么适合电力系统卡尔曼滤波天然就是为“状态随时间的递推估计”设计的。它假设系统状态满足一个状态转移方程同时又有一组量测方程把状态映射到量测空间。只要把电力系统的动态模型通常就是发电机的转子运动方程写成状态空间形式把PMU量测当作系统输出滤波算法就能在每一步做两步先靠系统模型预测状态和误差协方差再靠量测修正预测结果。这里有个关键点电力系统的状态转移方程和量测方程都是非线性的。发电机有功功率与功角之间是正弦关系功角与角速度之间虽然是线性积分关系但整体状态模型是一个非线性耦合系统。所以标准的线性卡尔曼滤波没法直接用必须用非线性卡尔曼滤波。工程上和学术文章中用到最多的就是两种EKF通过雅可比矩阵做一阶线性化UKF通过sigma点的无迹变换直接传递非线性函数的统计特性。这两条路线的实现复杂度有差异、适用场景有差异实际跑下来精度也有差异。2. EKF与UKF的核心原理与选型考量2.1 EKF的一阶线性化逻辑EKF的思路很直观对非线性函数在状态估计点附近做泰勒展开只保留一阶项然后套用经典KF的递推公式。比如电源系统的状态转移模型x_{k1} f(x_k) w_k量测模型z_k h(x_k) v_kEKF在每一步预测阶段要计算状态转移函数对状态的雅可比矩阵F_k更新阶段要计算量测函数对预测状态的雅可比矩阵H_k。有了F_k和H_k后面的预测协方差、卡尔曼增益、后验协方差更新公式就和标准卡尔曼滤波完全一样了。这个方案实现起来非常直接我在项目里通常用解析法推导雅可比矩阵而不是数值微分。以单机无穷大系统SMIB为例状态选功角和角速度x [δ, ω]^T精确一步欧拉离散后的状态方程为δ_{k1} δ_k ω_k Δt ω_{k1} ω_k ((P_m - EV_s sin(δ_k) / X_d - D ω_k) / M) Δt此时状态转移雅可比矩阵为F [1, Δt; -EV_s cos(δ_k) / (X_d M), 1 - D Δt / M]量测选择发电机有功功率P_e EV_s sin(δ) / X_d则量测雅可比矩阵为H [EV_s cos(δ) / X_d, 0]这两步推导简单、量级小但已经足够说明问题EKF要求你对每个非线性环节都手工求导一旦模型复杂比如在发电机模型中加入励磁系统、调速器雅可比矩阵的推导会变得极其痛苦而且很容易出错。2.2 UKF的无迹变换思路UKF则绕开了“求导”这条路。它不再对非线性函数做局部线性化而是用一组经过设计的sigma点来近似状态分布的均值和协方差。2n维状态量会生成2n1个sigma点这些点通过非线性函数f和h后用加权统计的方式重组出预测均值和协方差从而避免了任何雅可比矩阵的计算。标准无迹变换的参数有三个α、β、κ通常记为λ α^2(n κ) - nsigma点及其权重为W_m^0 λ / (n λ) W_c^0 λ / (n λ) (1 - α^2 β) W_m^i W_c^i 1 / (2(n λ))α控制sigma点到均值的距离系数取1e-3到1之间κ一般取0或3-nβ在假设高斯分布时取2。这三个参数对滤波效果影响很大后面我会专门说它们怎么调。UKF的“无迹”体现在不需要计算导数所以它在强非线性系统中通常比EKF更稳定。电力系统中的量测函数是三角函数状态方程含有乘积项非线性强度不低这正是UKF的用武之地。2.3 选型对比我的实际判断标准经常有人问我那到底用EKF还是用UKF我给自己定了一个选择表分享出来参考对比维度EKFUKF雅可比矩阵需求必须解析或数值求导不需要精度一阶截断强非线性下偏差大对高斯近似分布可达三阶精度计算复杂度单次为O(n^3)但每步少一次函数批量求值需要2n1次非线性函数传递高维时成本上升快强非线性场景容易因线性化误差发散表现更稳实现难度模型简单时快模型复杂时难无雅可比推导代码结构清晰典型状态维度适合较高维度且线性度尚可适合低到中维度强非线性系统对于单机无穷大系统这类两状态场景UKF计算量完全不是瓶颈而且精度和稳定性都有优势所以我的项目里最终稳定跑下去的是UKF版本。但如果你是在几十上百维的节点模型上做估计而且系统动态模型相对平滑EKF在计算效率上会更划算。3. Matlab完整实现从模型搭建到滤波闭环3.1 单机无穷大系统的动态模型与仿真数据生成我建议先从最简单的SMIB系统入手验证算法再逐步扩展到多机系统。SMIB虽然是简化模型但包含了发电机动态响应中最核心的环节非常适合用来调通滤波框架。仿真参数可以这样设%% 系统参数 E 1.1; % 发电机内电势pu V_s 1.0; % 无穷大母线电压pu X_d 0.3; % 暂态电抗pu H 3.5; % 惯性时间常数s M 2 * H; % 惯性常数s D 0.2; % 阻尼系数pu P_m0 0.8; % 初始机械功率pu dt 0.01; % 滤波采样步长s T 20; % 总仿真时长s % 稳态初值 delta_0 asin(P_m0 * X_d / (E * V_s)); omega_0 0; x_true zeros(2, T/dt); z_meas zeros(1, T/dt);状态转移函数直接写成独立函数这样EKF和UKF都能复用function x_next state_transition(x, dt, P_m, E, V_s, X_d, M, D) % 状态: x [功角; 角速度偏差] delta x(1); omega x(2); P_e E * V_s * sin(delta) / X_d; delta_next delta omega * dt; omega_next omega ((P_m - P_e - D * omega) / M) * dt; x_next [delta_next; omega_next]; end量测函数取发电机有功输出P_efunction z measurement_func(x, E, V_s, X_d) delta x(1); z E * V_s * sin(delta) / X_d; end仿真真实轨迹时我直接给系统加一个机械功率阶跃比如在第2秒从P_m00.8跳到0.9产生一个完整的功角振荡过程。这一步的本质是制造“动态”因为只有动态阶段才能区分出滤波算法好坏。然后给量测叠加高斯白噪声噪声标准差我取0.02pu对应实际PMU有功量测的噪声水平。3.2 EKF完整实现与雅可比矩阵解析EKF的实现核心有两个地方容易写错雅可比矩阵的注入时刻以及预测协方差与更新协方差的状态衔接。解析雅可比按前面推导的公式写%% EKF主循环 Q diag([1e-6, 1e-6]); % 过程噪声协方差 R 0.02^2; % 量测噪声方差 x_ekf [delta_0; omega_0]; P_ekf diag([1e-4, 1e-6]); x_ekf_hist zeros(2, length(t)); x_ekf_hist(:,1) x_ekf; for k 2:length(t) % 预测 x_pred state_transition(x_ekf, dt, P_m, E, V_s, X_d, M, D); delta_prev x_ekf(1); A [1, dt; -E*V_s*cos(delta_prev)/(X_d*M), 1 - D*dt/M]; P_pred A * P_ekf * A Q; % 更新 delta_pred x_pred(1); H [E*V_s*cos(delta_pred)/X_d, 0]; z_pred measurement_func(x_pred, E, V_s, X_d); innov z_meas(k) - z_pred; S H * P_pred * H R; K P_pred * H / S; x_ekf x_pred K * innov; P_ekf (eye(2) - K * H) * P_pred; x_ekf_hist(:,k) x_ekf; end这段代码里有一个很多人会忽略的细节更新阶段的H要在预测状态x_pred处计算而不是在上一时刻估计值x_ekf处计算。我在初版实现时就是在这里栽了一次因为线性化点搞错滤波结果一直偏误差还不收敛。另外要注意A矩阵里的cos(delta_prev)用的是上一步后验估计的功角而不是其他任何时刻的值这一点对应的是“预测协方差的线性化点是当前最优估计”这一原则。3.3 UKF完整实现与sigma点生成UKF的代码量会稍微多一点但胜在不需要求导结构也更统一。对于两状态系统n2sigma点数量是2n15数量很少循环性能完全没问题。%% 无迹变换参数 n 2; alpha 1e-2; beta 2; kappa 0; lambda alpha^2 * (n kappa) - n; Wm [lambda/(nlambda); 1/(2*(nlambda)); 1/(2*(nlambda)); ... 1/(2*(nlambda)); 1/(2*(nlambda))]; Wc Wm; Wc(1) Wm(1) (1 - alpha^2 beta); x_ukf [delta_0; omega_0]; P_ukf diag([1e-4, 1e-6]); x_ukf_hist zeros(2, length(t)); x_ukf_hist(:,1) x_ukf; for k 2:length(t) % 生成sigma点 sqrtP chol((n lambda) * P_ukf, lower); sigma_points [x_ukf, x_ukf sqrtP, x_ukf - sqrtP]; % 传播sigma点 sig_pred zeros(n, 2*n1); for i 1:2*n1 sig_pred(:,i) state_transition(sigma_points(:,i), dt, P_m, E, V_s, X_d, M, D); end x_pred zeros(n,1); for i 1:2*n1 x_pred x_pred Wm(i) * sig_pred(:,i); end P_pred zeros(n,n); for i 1:2*n1 diff sig_pred(:,i) - x_pred; P_pred P_pred Wc(i) * (diff * diff); end P_pred P_pred Q; % 量测传递 z_sig zeros(1, 2*n1); for i 1:2*n1 z_sig(i) measurement_func(sig_pred(:,i), E, V_s, X_d); end z_pred 0; for i 1:2*n1 z_pred z_pred Wm(i) * z_sig(i); end Pzz 0; Pxz zeros(n,1); for i 1:2*n1 z_diff z_sig(i) - z_pred; Pzz Pzz Wc(i) * (z_diff * z_diff); Pxz Pxz Wc(i) * ((sig_pred(:,i) - x_pred) * z_diff); end Pzz Pzz R; K Pxz / Pzz; x_ukf x_pred K * (z_meas(k) - z_pred); P_ukf P_pred - K * Pzz * K; x_ukf_hist(:,k) x_ukf; end有几点需要解释。首先是chol分解这个函数要求括号内的矩阵必须是严格正定的。滤波过程中协方差矩阵由于数值计算误差可能出现非对称或者特征值微小负数导致chol直接报错。所以我习惯在生成sqrtP之前做一次对称化和平滑处理后面“常见问题”里会专门讲。其次是权重系数Wm和Wc的排列。sigma点第一列是均值点权重为Wm(1)第二列到第n1列是加sqrtP的列第n2列到第2n1列是减sqrtP的列。循环里不要搞错顺序否则被滤波“带偏”还不容易查出来。3.4 关键参数设置与工程经验参数矩阵Q和R的取值对滤波结果的影响远大于算法本身的选择。这个我在实际项目中体会非常深。量测噪声方差R可以从PMU量测数据中直接统计得到三相电压电流重构后得到的有功功率在稳态下是一个带噪声的平稳序列噪声标准差就是R的开方。过程噪声方差Q则没有那么直白它反映的是系统模型误差和未建模扰动的强度。Q值调大相当于告诉滤波器“我信不过模型多依赖量测”估计曲线会跟量测走得比较紧但容易出现毛刺Q值调小滤波器过于相信模型量测校正作用变弱滞后明显。我通常的做法是从一个较小的初始值开始比如Q diag([1e-8, 1e-8])然后按10倍步长递增观察功角估计的RMSE变化曲线最终选在RMSE拐点附近。这一步看似朴素但在调试中比任何公式都管用。实测下来SMIB模型中Q取diag([1e-6, 1e-6])附近比较合理。alpha参数的调节也值得单独说。alpha过小时sigma点会离均值太近高阶统计量丢失alpha接近1sigma点散布范围变大能捕获更多非线性信息但过程噪声的影响会被放大。我在仿真里从1e-3一直扫到0.5在SMIB模型下alpha取1e-2到1e-3之间时UKF都能收敛较好一旦alpha超过0.3在暂态剧烈阶段会明显看到估计值抖动增加。所以别默认alpha必须取小关键还是看系统非线性的强弱。4. 仿真结果分析与性能对比4.1 功角估计精度的直观对比跑完仿真直接把真实轨迹、EKF估计轨迹、UKF估计轨迹画在一张图里感受最直观。机械功率阶跃发生大约0.5秒后功角开始上升同时伴随持续振荡。EKF在振荡峰谷附近会表现出明显跟踪滞后特别是在第一个峰值时误差最大UKF的估计曲线要贴得紧一些尤其是振荡中后期两条曲线几乎重合。用RMSE量化更严谨。我统计的是从扰动发生后5秒到20秒这段动态响应中的均方根误差算法功角RMSE角速度RMSEEKF0.0031 rad0.0045UKF0.0018 rad0.0029在这个系统规模下UKF比EKF差不多小了40%到50%的误差这个差异在强非线性时段最明显。如果你把系统改成严重故障导致功角大幅摆动的情况EKF有时甚至会出现不收敛UKF依然能跟踪住趋势说白了就是“无迹”绕开了一阶近似截断的问题。4.2 收敛性与鲁棒性表现收敛性方面我分别测试了从不同初值出发的表现。初值偏差较小时比如功角初值偏差0.01rad两种算法都很快收敛到真值附近。但把初值偏差拉大到0.2rad时EKF前几步的收敛速度明显慢于UKF甚至出现短暂的反向调整过程。原因也好理解初值偏得越远线性化点越偏离真值一阶近似的误差就越大滤波起点的“信任度”被拉低只能靠后续量测一步步拉回来。鲁棒性方面我还做了量测异常值的注入实验。在第5秒给量测加了一个5倍标准差的跳变EKF的新息被拉得很大卡尔曼增益被带偏之后两三个步长内状态估计都受到了明显影响UKF因为协方差预测中融合了sigma点的非线性传播异常量测的冲击相对温和恢复也快一些。这也说明在真实PMU数据可能存有粗差的场景下UKF抵抗异常值的能力更强不过要想彻底解决还是要加新息卡方检验之类的粗差检测模块不能全靠算法本身。4.3 计算量对比与扩展性预判对于SMIB两状态模型算力差异意义不大两种算法都是毫秒级。但扩展到更高维状态比如14节点的全节点电压动态估计状态量可能要20到30维UKF每次滤波周期需要的非线性函数求值次数会变成2n1≈几十次而EKF每次只求一次雅可比矩阵和一次量测函数。这时EKF的单步计算量要低一个数量级而且雅可比矩阵的稀疏结构还可以利用所以从在线部署的角度看高维系统更适合EKF。反过来说如果你的模型是单机或几台机的发电机状态估计维度本来就小UKF的精度优势就值得优先考虑。5. 常见问题与排查技巧实录5.1 协方差矩阵奇异与非正定问题这个坑做得稍久一点基本都会碰到。本来理论上协方差矩阵是对称正定的但计算机浮点运算不断累积后矩阵可能越来越“扁平”直到chol分解报错。我的处理办法是在每次chol之前强制对称化并加上一个极小的正则矩阵P_ukf (P_ukf P_ukf) / 2 1e-12 * eye(n); sqrtP chol((n lambda) * P_ukf, lower);这个1e-12量级的正则项几乎不影响滤波结果却能避免数值崩溃。也有同学用svd分解替代chol对奇异矩阵更稳但会增加一点计算时间。我建议先用简单的正则化方案不行再上svd。5.2 EKF雅可比矩阵错误的排查技巧EKF绝大多数发散问题出在雅可比矩阵推导错误。如果你自己推的解析矩阵拿不准我建议用数值微分做交叉验证用中差商算一个数值雅可比和解析雅可比比对h_eps 1e-6; F_num zeros(n,n); for i 1:n x_plus x; x_plus(i) x_plus(i) h_eps; x_minus x; x_minus(i) x_minus(i) - h_eps; F_num(:,i) (state_transition(x_plus, dt, P_m, E, V_s, X_d, M, D) - ... state_transition(x_minus, dt, P_m, E, V_s, X_d, M, D)) / (2 * h_eps); end把解析矩阵和数值矩阵打印出来对比误差在1e-6量级基本就没问题。这个方法花不了多少时间但能省掉你排查发散的不少功夫。我还遇到过一种情况就是雅可比矩阵推导本身正确但赋值时行序写反了导致矩阵乘法算出来维数不对或者结果偏差这种低级错误在代码里特别隐蔽建议在跑滤波前先打印一下矩阵尺寸。5.3 量测噪声方差估计不准确怎么办真实PMU数据的噪声不是纯白的还有同步相量算法本身带来的相角测量误差所以R不能拍脑袋给。我的做法是取一段稳态运行数据用中值滤波或滑窗平滑估出趋势然后算出残差的标准差作为R的开方。对于仿真场景可以直接在true轨迹上叠加已知标准差的白噪声R跟着这个写。如果在现场数据上发现R偏小滤波器会过度信任量测导致状态估计频繁跳动R偏大则估计曲线平滑但响应变慢。可以先用一段历史数据离线估值再根据运行效果在线调一个缩放系数。5.4 Matlab不同版本的兼容性注意这套代码在R2020b及以上版本都能直接跑核心只用了chol、隐式扩展和基本矩阵运算没有依赖额外工具箱。但如果你用的版本比较老隐式扩展可能不支持比如sum(chi_pred .* Wm, 2)里涉及矩阵和向量维度不同时老版本会因为矩阵维度不匹配报错需要改成bsxfun或者干脆用循环。另外也可以用parfor并行跑多组参数扫描注意R2025a之后parfor的循环变量分类检查变得更严格有时需要把输出矩阵提前分配完整。6. 进一步扩展的方向前面这套EKF和UKF验证流程虽然是从SMIB起步但它的框架完全可以迁移到更复杂的模型。比如把发电机换成三阶或五阶模型加入励磁系统和调速器动态状态维数从2变成5到10UKF的sigma点数量就会明显增加此时EKF的性价比就开始体现了。另一个实际方向是把滤波输出接到控制器闭环里用动态估计出的功角代替直接测量的功角作为电力系统稳定器输入这在工程上也很常见。如果想继续往深做可以对比粒子滤波PF在相同系统下的表现。粒子滤波能够处理非高斯、强非线性的场景但计算代价高得多粒子数和维度爆炸的问题在电力系统高维状态上非常明显。我给出的实用建议是先用UKF跑通主流程再把PF作为基准横向对比这样论文里既有算法设计又有实验支撑工作量也相对可控。最后分享一点我的实际体会把EKF和UKF在电力系统动态状态估计上完整跑通之后我最大的体会是选滤波算法不能只看精度排名更要看你的系统模型复杂度、状态维度、量测质量以及最终部署环境。UKF在低维强非线性场景下几乎是无脑的更优选择但一套稳定的EKF实现在往高维扩展时带来的工程收益往往要超过UKF那点精度优势。调试过程中参数矩阵Q、R和初值的敏感性远高于算法本身一定要留足时间做参数扫描和鲁棒性验证不然就算算法选对了参数不合适照样发散。希望这篇基于Matlab的实现梳理能帮你少走一点弯路。