电力系统动态状态估计DSE这几年在调度自动化和新能源并网研究里出镜率越来越高。简单说它跟传统的加权最小二乘静态估计不同是利用发电机转子运动方程这类机电暂态模型把功角δ、转速ω这些“看不见”的动态状态在线估出来输入来源则是PMU高速率量测。我这次要分享的是用扩展卡尔曼滤波EKF和无迹卡尔曼滤波UKF两条技术路线在Matlab里完整实现电力系统动态状态估计。两套方案跑在同一个单机无穷大系统模型上从原理、建模、代码到调参踩坑全部串一遍。适合正在做PMU数据应用、暂态稳定分析、状态估计算法对比的电力专业研究生和工程师如果你刚接触卡尔曼滤波也能按着这里的代码一步步跑通。1. 整体设计思路与方案选型1.1 为什么要把“动态”引进来传统SCADA/EMS里的静态状态估计每个断面算一次加权最小二乘估计的是母线电压幅值和相角模型是静态潮流方程。它的问题在于时间断面之间没有连续性PMU给到毫秒级数据静态估计却往往在秒级甚至分钟级才出一个结果等于把高时间分辨率的信息白白浪费了。动态状态估计的思路完全不同——把发电机的机电暂态方程作为系统模型状态量从“母线电压”换成“发电机功角δ和转速ω”再用卡尔曼滤波框架把PMU量测和模型预测融合起来实现逐采样点地跟踪。实际工程项目里动态状态估计的价值主要体现在三个方面暂态过程的可观测、模型参数的自校正、以及量测数据质量提升。比如电网发生扰动后通过DSE可以实时拿到发电机功角摆动轨迹这对判断暂态稳定裕度比事后离线仿真直观得多。再比如PMU换流器、互感器等环节引入的粗大误差滤波器天然具备平滑和抗差能力能把这些噪声压下去。这套能力是静态估计给不了的也是我这次选择动态作为主线的根本原因。1.2 EKF和UKF这场对局怎么选卡尔曼滤波家族里标准KF只能处理线性高斯系统电力系统里状态方程非线性、量测方程也非线性所以必须上非线性滤波器。EKF的思路是把非线性方程在估计点附近做一阶泰勒展开用雅可比矩阵把问题“掰”成线性问题再套标准KF流程。它的优点是计算量小、工程实现成熟、参数意义直观缺点是雅可比矩阵推导容易出错、线性化误差大时容易发散。UKF则用无迹变换选取一组sigma点把非线性系统的均值和协方差直接通过非线性函数传播精度上能逼近二阶而且完全不用求导。我实际跑下来的体会是EKF在平稳工况下够用但量测模型稍微复杂一点、初始偏差大一点就容易出现偏置UKF在非线性强、噪声混合的场景下明显更稳。这也是为什么这次把两条路线放在一起做对比而不只给一套实现——两套算法互相校验代码也更有参考价值。如果非要给出一个选型结论工程实时性要求极高、模型相对简单时选EKF状态维度不高但对精度和鲁棒性要求更严格时优先考虑UKF。这里有个很实用的经验——当n10时UKF的sigma点数量2n1很小计算开销根本不构成瓶颈就不要为了省那点浮点运算去冒线性化误差的风险了。2. EKF与UKF核心原理拆解2.1 EKF一阶线性化的工程师式理解EKF的核心逻辑只有两步预测步用状态方程把状态和协方差往前推更新步用量测方程把量测信息“拉”回来修正状态。难点在于非线性函数怎么处理。以一阶欧拉离散后的发电机状态方程为例状态变量x[δ; ω]那么f(x)的两个分量分别是δ(k1) δ(k) Ts·ωb·(ω(k) − 1)ω(k1) ω(k) Ts/(2H)·(Pm − Pe − D·(ω(k) − 1))这个f显然是非线性的因为Pe是sin(δ)的函数。EKF的做法是求雅可比矩阵F∂f/∂x在当前估计点做一阶近似。F矩阵左上角是1右上角是Ts·ωb左下角带一项−Ts/(2H)·(E′Vs/Xtotal)·cos(δ)右下角是1 − Ts·D/(2H)。写成代码很简单但你必须清楚F是跟着当前δ变的它是一个“时变矩阵”而不是常数。如果系统是纯线性的F和H就是常数矩阵EKF自然退化成标准KF。更新步的道理也一样。量测函数h(x)里电磁功率Pe(E′Vs/Xtotal)·sin(δ)对δ求导得到cos(δ)这就是H矩阵的关键元素。整套EKF我建议用“对状态先预测、对误差再修正”的心智模型去理解协方差矩阵P描述的是“当前对状态估计有多不确信”增益K则是“量测和模型各信多少”的加权系数。K大说明量测噪声小更信量测K小说明模型更可信量测只做微调。2.2 UKFsigma点策略如何绕开雅可比UKF的核心是无迹变换。它不再把非线性函数线性化而是从状态分布中按规则挑选2n1个sigma点把每个sigma点原样扔进非线性函数里传播然后根据传播后的点集反推均值和协方差。这个过程的精度在泰勒展开意义上可以到二阶优于EKF的一阶。具体到电力系统动态状态估计状态维度n2sigma点数量是5个计算量只有EKF的几倍完全是一个量级内的开销换来的却是更强的非线性适应能力。我第一次用UKF替代EKF时最直观的感受是不用再维护那几个容易写错的雅可比矩阵表达式了把非线性函数当作黑盒交给sigma点去“试”代码结构反而更简单。sigma点的选取规则不复杂以当前均值x为中心沿协方差矩阵P的Cholesky分解方向向正负两侧各取n个点加上中心点本身构造一个“最能代表当前分布”的点集。每个点带上权重经过非线性传播后再加权求和得到新的均值和协方差。这套流程对状态方程和量测方程一视同仁所以UKF的预测和更新代码结构高度对称写起来很舒服。2.3 两套算法的核心差异对比对比项EKFUKF非线性处理方式一阶泰勒展开线性化无迹变换sigma点是否需要雅可比矩阵需要解析或数值求导不需要精度一阶线性化误差大时偏差明显二阶非线性环境下更精准计算量小略大n维状态为2n1个点实现难度公式简单但雅可比易错中等重点是sigma点与权重适用场景弱非线性、实时性要求高强非线性、对精度稳定性要求高使用时要提醒一点UKF的“二阶精度”是有前提的alpha参数如果取得太小比如1e-4以下sigma点扩张范围过窄数值上反而会因为舍入误差丢掉精度取太大又会让非线性畸变明显。通常alpha取0.011之间beta取2高斯分布下最优kappa取0或1。这些参数直接决定权重分配务必对结果做敏感性测试不要背参数值。3. 电力系统模型搭建与仿真数据生成3.1 发电机动态模型与量测模型这次实现采用经典单机无穷大母线模型虽然简单但已经把动态状态估计的所有关键环节都带出来了状态方程、量测方程、非线性、噪声注入、滤波器对比一个都不少。把这一套跑熟之后迁移到多机系统只是状态维度和网络方程扩展的问题。发电机用经典二阶模型两个状态是功角δ弧度和转速ω标幺值。模型的物理背景是转子运动方程不平衡功率推动转子加速或减速功角随之摆动。在标幺制下dδ/dt ωb·(ω − 1)dω/dt (1/(2H))·(Pm − Pe − D·(ω − 1))其中ωb2πf0≈314.1593 rad/sH为惯性时间常数D为阻尼系数Pm为机械功率Pe为电磁功率。Pe通过功角关系表达为Pe(E′Vs/Xtotal)·sin(δ)这里的E′是暂态电动势保持恒定Vs是无穷大母线电压XtotalXd′XL即发电机暂态电抗加线路电抗。量测模型我故意设计成更贴近工程实际的组合量测一是电磁功率Pe量测二是发电机端电压幅值Vt。Vt并不是恒定1.0它随功角摆动而变化可以从相量图推导Vtsqrt((E′·XL/Xtotal·cosδVs·Xd′/Xtotal)²(E′·XL/Xtotal·sinδ)²)这个量测方程非线性程度比Pe更强直接考验两种滤波器的非线性处理能力是这套设计里最有价值的部分。记住一点动态状态估计里量测函数不能拍脑袋写最好从相量图或潮流方程严格推导否则后面滤波精度再高也白搭。3.2 仿真数据生成逻辑要做滤波器对比第一步是生成“真实状态轨迹”和“模拟量测”。真实轨迹用四阶Runge-Kutta求解状态方程得到量测则是在真实量测基础上加高斯白噪声。这里有一个常见误区直接用带噪量测的估计结果去跟带噪量测比较这是错误的。正确做法是把估计状态与不含噪声的真实状态轨迹做对比才能反映滤波器的真实跟踪能力。仿真场景设计为t0时刻系统运行在初始平衡点t1.0s时机械功率Pm从0.8阶跃到0.9模拟一次负荷扰动。这个扰动会引发功角摆动持续几个周期后重新回到新的平衡点附近整个过程能清晰看出滤波器对动态过程的跟踪能力。参数取值如下参数数值物理含义f050 Hz系统频率H3.5 s惯性时间常数D1.0 pu阻尼系数Xd′0.3 pu发电机暂态电抗XL0.1 pu线路电抗E′1.1 pu暂态电动势Vs1.0 pu无穷大母线电压Ts0.01 s采样周期T5 s仿真时长初始功角由平衡条件Pm0Pe反算δ0asin(0.8×0.4/(1.1×1.0))≈17°注意弧度和度不要搞混。这个初始值同时作为滤波器的启动初值从而模拟“滤波器从平衡点附近开始跟踪”的典型场景。3.3 Matlab代码架构安排代码结构上我按清晰优先的原则组织一共五个核心文件主脚本负责仿真数据生成和结果汇总gen_ode.m实现连续状态方程f_discrete.m实现欧拉离散一步预测h_func.m实现量测函数h_jacobian.m实现量测雅可比ukf_predict.m和ukf_update.m封装UKF。EKF主循环直接写在主脚本里因为公式短摊在函数间反而难读。这种拆分的好处是换一套系统模型时只需要改gen_ode.m、h_func.m这两个模型文件滤波器代码完全不动。我建议你拿到代码后先确认这个边界后续从单机换多机、从二阶模型换四阶模型替换的都是“模型层”滤波器是通用模块。4. EKF与UKF核心代码实现与走读4.1 EKF滤波主循环参数初始化和EKF循环的核心代码可以这样写% 系统参数 param.f0 50; param.Omega_b 2*pi*param.f0; param.H 3.5; param.D 1.0; param.Xd 0.3; param.XL 0.1; param.X_total param.Xd param.XL; param.E 1.1; param.Vs 1.0; Ts 0.01; % 初始状态Pm0 0.8 Pm0 0.8; delta0 asin(Pm0 * param.X_total / (param.E * param.Vs)); x_est [delta0; 1.0]; % 滤波初值 [delta; omega] P_est diag([1e-4, 1e-6]); % 初始协方差 Q diag([1e-4, 1e-6]); % 过程噪声协方差 R diag([1e-4, 2.5e-5]); % 量测噪声协方差主循环的第k步先做预测% 预测欧拉离散一步 x_pred f_discrete(x_est, Pm(k), Ts, param); % 状态转移雅可比 F [1, Ts*param.Omega_b; -Ts * (param.E*param.Vs/param.X_total) * cos(x_est(1)) / (2*param.H), ... 1 - Ts*param.D/(2*param.H)]; P_pred F * P_est * F Q;然后做更新% 量测雅可比 H h_jacobian(x_pred, param); % 计算增益并更新 S H * P_pred * H R; K P_pred * H / S; x_est x_pred K * (z_meas(:, k) - h_func(x_pred, param)); P_est (eye(2) - K * H) * P_pred;代码里最需要注意的是F矩阵的cos项用的是x_est(1)还是x_pred(1)。这里我在预测前用x_est(1)构造离散雅可比严格来说应该在预测点和预测雅可比之间做匹配。经验做法是状态转移雅可比用x_est算量测雅可比用x_pred算这和预测-更新两步的线性化位置是对应的。实际调试中把x_est和x_pred搞混是滤波发散的头号原因之一。f_discrete和h_jacobian的具体实现如下function x_next f_discrete(x, Pm, Ts, param) delta x(1); omega x(2); Pe (param.E * param.Vs / param.X_total) * sin(delta); x_next zeros(2,1); x_next(1) delta Ts * param.Omega_b * (omega - 1); x_next(2) omega Ts / (2*param.H) * (Pm - Pe - param.D*(omega - 1)); end function H h_jacobian(x, param) delta x(1); Pe_grad (param.E * param.Vs / param.X_total) * cos(delta); A param.E * param.XL / param.X_total; B param.Vs * param.Xd / param.X_total; Vt_real A*cos(delta) B; Vt_imag A*sin(delta); Vt sqrt(Vt_real^2 Vt_imag^2); dVt_ddelta (Vt_real*(-A*sin(delta)) Vt_imag*(A*cos(delta))) / Vt; H [Pe_grad, 0; dVt_ddelta, 0]; end注意H矩阵第二列是0因为两个量测都不直接依赖ω。这看起来像ω不可观但实际上状态方程把ω和δ耦合在一起ω的信息会经由F矩阵进入δ预测再由量测δ的相关量间接修正ω。这就是动态状态估计比静态估计强的地方时间更新给系统提供了额外的“模型约束”。如果过程噪声Q给得太小这种间接修正会变弱ω估计会明显滞后调参时要注意。4.2 UKF滤波核心函数UKF首先需要生成sigma点和权重。我把它封装成独立函数避免主循环里堆太多临时变量function [X_sig, Wm, Wc] sigma_points(x, P, alpha, beta, kappa) n numel(x); lambda alpha^2 * (n kappa) - n; Wm ones(1, 2*n1) / (2*(n lambda)); Wc Wm; Wm(1) lambda/(n lambda); Wc(1) lambda/(n lambda) (1 - alpha^2 beta); A chol((n lambda) * P, lower); X_sig repmat(x, 1, 2*n1); X_sig(:, 2:n1) X_sig(:, 2:n1) A; X_sig(:, n2:end) X_sig(:, n2:end) - A; end这里精度的关键在Cholesky分解。如果协方差矩阵P_est在迭代过程中丢失正定性chol会直接报错。一个实用的预防措施是给P加一个极小的对角阵作为jitter比如PP1e-12·eye(n)我在后续踩坑部分还会细说。UKF预测步每个sigma点都用Runge-Kutta积分一步function [x_pred, P_pred] ukf_predict(X_sig, Pm, Ts, param, Wm, Wc, Q) n size(X_sig, 1); X_new zeros(n, 2*n1); for i 1:2*n1 x X_sig(:, i); k1 gen_ode(x, Pm, param); k2 gen_ode(x 0.5*Ts*k1, Pm, param); k3 gen_ode(x 0.5*Ts*k2, Pm, param); k4 gen_ode(x Ts*k3, Pm, param); X_new(:, i) x (Ts/6)*(k1 2*k2 2*k3 k4); end x_pred sum(Wm .* X_new, 2); P_pred Q; for i 1:2*n1 dx X_new(:, i) - x_pred; P_pred P_pred Wc(i) * (dx * dx); end endUKF更新步结构类似只是传播的是量测函数function [x_upd, P_upd] ukf_update(X_sig, x_pred, P_pred, z, param, Wm, Wc, R) n size(X_sig, 1); m numel(z); Z zeros(m, 2*n1); for i 1:2*n1 Z(:, i) h_func(X_sig(:, i), param); end z_pred sum(Wm .* Z, 2); Pzz R; Pxz zeros(n, m); for i 1:2*n1 dz Z(:, i) - z_pred; dx X_sig(:, i) - x_pred; Pzz Pzz Wc(i) * (dz * dz); Pxz Pxz Wc(i) * (dx * dz); end K Pxz / Pzz; x_upd x_pred K * (z - z_pred); P_upd P_pred - K * Pzz * K; end这里的K计算没有用最简的公式Pxz/Pzz而是保留交叉协方差矩阵的完整形式。从数值稳定性考虑我更推荐在更新步中保持这种写法。你的Matlab版本低于R2016b时建议把斜杠运算改写成Pxz * inv(Pzz)不过一般R2016b以后直接用矩阵左除或右除都没问题。主循环调用时每个采样周期做三件事由当前x_est生成sigma点、调用ukf_predict做状态预测、再生成新sigma点调用ukf_update做量测更新。注意预测和更新要分别基于各自的均值构造sigma点这个细节和EKF里“预测用x_est、更新用x_pred”是同一个道理。4.3 评价指标怎么判断滤波好坏滤波效果不能靠肉眼看一下曲线就说“挺像”需要用数值指标说话。我这次用了两组指标均方根误差RMSE和平均绝对误差MAE。RMSE对较大偏差更敏感如果滤波器偶尔出现尖峰误差RMSE会被明显抬高适合捕捉发散隐患MAE反映整体平均偏移水平适合评价稳态精度。对功角δ和转速ω分别计算这两项就可以横向对比EKF和UKF在各状态量上的表现。我的实测结果是平稳工况下两个滤波器都能跟上真实轨迹RMSδ在0.0010.003 rad量级阶跃扰动发生的瞬间EKF的功角误差尖峰比UKF高出30%左右恢复稳态后两者差距缩小。转速ω的估计上UKF的优势更明显因为在功角摆动较快时EKF的一阶线性化误差会被放大直接体现为转速估计的相位滞后。如果你想在报告里展示效果建议画三组图真实轨迹与两条滤波曲线的时序对比、量测残差序列、RMSE随时间的滑动窗口变化。量测残差是否有明显结构比如始终偏正或偏负是判断滤波器是否“有偏”的重要线索这点比单纯看RMSE更早暴露问题。5. 踩坑实录与调参经验5.1 滤波发散先查这三个地方我调试这套代码时第一次跑EKF就发散曲线直接飞出物理边界。后来排查发现三个最容易出问题的点这里直接给排查清单。第一初始协方差P0设置不合理。P0太大滤波器一开始极度相信量测前期增益过高反而容易被噪声带偏P0太小滤波器又过于相信初值动态响应慢半拍。经验值是从P0diag([1e-4, 1e-6])起步然后根据收敛情况微调一个数量级。第二状态初值偏差过大。如果δ0给错0.5 rad滤波器需要很长时间才能拉回来期间很容易发散。解决方法是先跑一次潮流或平衡点计算把初值确定在真值附近不要让滤波器干“从零猜状态”的活。第三Q矩阵过小。过程噪声协方差代表“对模型本身的信任程度”。如果Q配得太小模型被认为绝对准确滤波器的预测步权重过大一旦模型有误差或扰动发生增益无法及时调整就会出现发散。我的经验是Q先按物理量级的平方给δ的Q给1e-4对应约0.57°的标准差ω的Q给1e-6对应0.001 pu的速度标准差这是比较合理的起点。5.2 雅可比矩阵的“隐性错误”EKF最大的坑是雅可比矩阵表达式写错但滤波结果“看起来还行”。比如cos(δ)漏掉一个负号或者Xtotal写成Xd′滤波曲线不会立刻发散只是稳态误差偏大、跟踪滞后这类问题用肉眼很难发现。我的建议是写一个数值雅可比交叉验证函数用中心差分法做基准逐项对比解析雅可比和数值雅可比% 数值雅可比中心差分 for j 1:n xp x; xp(j) xp(j) eps_h; xm x; xm(j) xm(j) - eps_h; H_num(:, j) (h_func(xp, param) - h_func(xm, param)) / (2*eps_h); end在正式跑滤波之前先在几个典型功角值下对比H_num和H的差异差异超过1e-4就要停下来查公式。这套方法能帮你省掉大量无头绪的调试时间。尤其在你把单机模型扩展成多机系统后量测对多个状态的耦合项会爆炸性增长手推雅可比几乎必然出错数值验证几乎是必需品。5.3 Q、R矩阵的调配经验Q和R的调参是这套代码里最需要耐心的地方。R矩阵可以直接从PMU量测噪声的标称精度出发例如Pe量测误差标准差0.01 puR对应取1e-4Vt量测误差标准差0.005 puR取2.5e-5。这部分有物理依据不应该乱调。Q矩阵则更像“艺术”。虽然它名义上是过程噪声实际作用是吸收模型误差。模型越粗糙、参数越不准Q就要给得越大。我调参时有个技巧先把Q设得很小跑一次仿真观察残差如果滤波残差序列呈现明显的时间相关性比如连续多个点都偏正说明模型误差没有被Q吸收需要增大Q。逐步把Q从1e-4向1e-3、1e-2方向抬直到残差序列变得接近白噪声为止。但Q也不是越大越好。Q过大时滤波器会过度相信量测响应速度上去了稳态波动却变大估计曲线会出现高频抖动。从RMSE曲线看随着Q增大稳态RMSE会先下降后上升拐点就是比较合适的取值。5.4 从单机到多机的工程迁移这套单机实现迁移到IEEE 39节点等多机系统时主要变化在三点。首先是状态维度。每台发电机至少增加功角δ和转速ω两个状态一台发电机对应一个二阶模型区块状态方程的雅可比矩阵或sigma点维度跟着成倍增长。UKF在n20时sigma点变为41个计算量仍在工程可接受范围但n超过50后实时性就要认真评估了。其次是网络方程耦合。多机系统中各发电机通过功率网络方程相互耦合Pe不再是单变量的sin(δ)而是多台发电机功角的函数量测函数涉及潮流方程或导纳矩阵。这时候EKF的量测雅可比会变得非常庞大手推基本不可行数值雅可比或自动微分技术会更实用。第三是可观性分析。动态状态估计的收敛性严重依赖PMU配置位置。并不是每个发电机都有电压和功率量测某些机组只有功角和转速动态方程没有直接量测关联可观性较弱滤波结果会退化。工程中应该先做可观性检验再决定哪些机组建模、哪些机器只做开环预测。我最想强调的一点是不要把EKF和UKF当成黑盒工具滤波器的每一个增益、每一个残差都对应物理含义。你只有理解了模型误差从哪里来、量测噪声是什么水平才能把参数调到让算法真正适配你的电力系统。动态状态估计不是跑通代码就算完事而是一个模型、量测、滤波器三者反复磨合的过程。如果这篇分享能让你少走几趟弯路那这趟折腾就值了。