做电力系统动态状态估计DSE课题我猜你迟早会撞上 EKF 和 UKF 这两个名字——它们几乎成了非线性滤波的代名词。我在 Matlab 里把两套滤波器的代码完整落地跑通仿真、做对比、调参数过程中踩了不少坑也把每个坑背后的原因摸清楚了大半。这篇就把整个项目的设计思路、算法原理、Matlab 实现细节和排错经验一次讲透无论你是正在写课程代码、准备毕业论文还是刚接手状态估计课题想快速跑通一个基线方案都可以直接照着做。先交代一下这个项目在干什么电力系统动态状态估计简单说就是在 PMU 同步相量量测的基础上结合发电机机电暂态模型实时递推估计发电机功角、转速这类动态状态。整个项目用 Matlab 实现核心是两套非线性滤波器扩展卡尔曼滤波 EKF 和无迹卡尔曼滤波 UKF。除了把算法跑通我还做了精度对比、参数敏感性分析和发散问题排查下面的内容都是我实际验证过的不是教科书式的空谈。1. 项目拆解动态状态估计与两种滤波器的定位1.1 电力系统动态状态估计在解决什么问题传统电力系统状态估计大多是静态的基于 SCADA 系统的量测数据按秒级甚至分钟级周期计算全网状态输出给 EMS 做在线分析。但这些年新能源比例上来以后系统的动态行为明显变复杂频率波动、暂态过程、次同步振荡这类现象越来越频繁秒级静态估计根本跟不上机电暂态过程的变化速度。动态状态估计的意义就在这里它在状态方程里引入发电机的转子运动方程用 PMU 提供的毫秒级同步相量量测实时跟踪功角、转速这类变量的变化轨迹。从控制理论角度看这就是一个典型的状态观测问题。对象是发电机的动态模型输入是机械功率和励磁电压输出是 PMU 量测到的电压、电流相量状态则是功角、角速度、暂态电动势。PMU 本身有量测噪声发电机模型也有不确定性所以要把贝叶斯递推框架套上去——预测加更新用最新量测不断修正模型推算结果。这个问题的本质是系统模型和量测模型都非线性经典卡尔曼滤波 KF 只适用于线性高斯系统所以必须上非线性滤波方法。1.2 为什么偏偏是 EKF 和 UKF非线性滤波的选项其实不少粒子滤波、H∞滤波、平方根滤波都能做但 EKF 和 UKF 是入门和研究最常用的两个原因也很实际粒子滤波需要大量粒子才能保证精度计算量对在线电力系统应用来说通常不可接受EKF 和 UKF 的计算量都在毫秒级以下适合实时递推。EKF 的思路是线性化在每个时间点把非线性状态方程和量测方程做一阶泰勒展开得到雅可比矩阵然后套用经典 KF 的公式。它的优点是原理直观、计算量最小缺点也明显——强非线性场景下线性化误差会累积而且手推雅可比矩阵非常容易出错。UKF 完全是另一条路不去线性化函数而是用无迹变换通过一组精心构造的 Sigma 点来捕获状态分布的一阶矩和二阶矩。对高斯分布来说理论上精度能达到三阶而且不需要求任何导数实现流程是机械式的。在电力系统这种量测方程包含功率潮流计算的场景里UKF 往往比 EKF 稳健因为它不会因为某个局部线性化点选得不好而出大错。把这两个放一起做成了完整的对比实验一是为了论文里有理有据地分析精度差异二是想真正理解非线性滤波的两大流派逼近函数还是逼近分布。1.3 整体方案怎么分层项目代码我没有写成一个大脚本堆到底而是分了四层系统模型层发电机动态方程、网络量测方程统一封装成状态转移函数 f 和量测函数 h量测生成层根据真实状态轨迹叠加高斯噪声生成模拟 PMU 量测序列滤波算法层EKF 和 UKF 两个独立实现接口完全一致评估层计算 RMSE、绘制状态估计曲线和误差曲线这样分层的好处很明显以后换系统模型比如从单机无穷大系统扩展到多机系统滤波算法层的代码完全不用动想新加一种滤波器比如容积卡尔曼滤波 CKF只需要照着同一套接口写实现就行。我在做这个项目的过程中反复体会到代码解耦省下来的调试时间比写代码本身多得多。2. 算法核心细节与实操要点2.1 EKF线性化的代价和边界EKF 的递推公式不复杂核心就五个式子。设状态预测值为 x_pred协方差预测为 P_pred滤波器更新后的状态为 x_est协方差为 P_estx_pred f(x_prev) P_pred F * P_prev * F Q K P_pred * H / (H * P_pred * H R) x_est x_pred K * (z - h(x_pred)) P_est (I - K * H) * P_pred这里的 F 和 H 分别是状态转移函数和量测函数对状态的雅可比矩阵。真正麻烦的地方也在这里。我用一个很经典的二阶发电机模型来解释状态取功角 δ 和角速度偏差 Δω状态方程离散化后是δ(k1) δ(k) ωn * Δω(k) * dt Δω(k1) Δω(k) (Pm - Pmax*sin(δ(k)) - D*Δω(k)) / (2H) * dt其中 ωn 是同步转速 2π50Pmax 是发电机的最大电磁功率标幺值D 是阻尼系数H 是惯性时间常数。量测可以取功角 δ 和电磁功率 Pe Pmaxsin(δ)h1 δ h2 Pmax*sin(δ)对应的雅可比矩阵手推起来还算轻松F [1, ωn*dt; -Pmax*cos(δ)*dt/(2H), 1 - D*dt/(2H)] H [1, 0; Pmax*cos(δ), 0]但这里藏着一个特别容易踩的坑H 矩阵第二行是对 δ 求偏导有的同学会把 Pmaxcos(δ) 写成 Pmaxsin(δ) 或者丢掉前面系数一错就全错。而且一旦换成三阶模型加上暂态电动势 E q 那一维Pe 的表达式会包含 Id、Iq 这些与 δ 耦合的中间变量雅可比推导的复杂度会直线上升。我的建议是手推完雅可比以后一定用数值差分交叉验证一下——Matlab 里写一个两点差分函数比如 (f(xh)-f(x-h))/(2h)对比两者是否一致。这个验证步骤替我省下了至少两天调试时间。EKF 的另一个边界条件是线性化点必须在真实状态附近。如果初始状态给得很偏或者过程噪声 Q 设置得太小导致模型预测误差一直压不下来线性化点就会跑偏滤波结果自然越来越差。这是 EKF 在很多实际问题里“莫名发散”的根源之一。2.2 UKF无迹变换与三个关键参数UKF 的核心是无迹变换。对 n 维状态 x均值 x̄协方差 P构造 2n1 个 Sigma 点每个点配一个权重然后让这些点通过非线性函数传播再用传播后的点重新计算均值和协方差。公式长这样λ α^2*(nκ) - n χ0 x̄ χi x̄ sqrt((nλ)*P)_i, i 1,...,n χ(in) x̄ - sqrt((nλ)*P)_i, i 1,...,n权重W0m λ/(nλ) W0c λ/(nλ) (1-α^2β) Wim Wic 1/(2*(nλ))参数选择有固定套路。α 控制 Sigma 点围绕均值的散布程度一般取 1e-3 到 1 之间β 和状态分布的先验有关高斯分布取 2 最优κ 通常取 0 或 3-n对高斯分布来说影响很小。这里我提醒一个很多人没注意到的问题λ 不能取 -n否则 (nλ) 乘完协方差矩阵之后为零矩阵Sigma 点全部坍缩到均值上滤波器直接失效。我习惯设 α0.1、κ0这样对 n3 的三阶系统nλ 算出来是正的小量不会出幺蛾子。无迹变换的思想用一句生活化的话说如果直接算一条复杂曲线的积分很难那就多取几个样本点用样本的分布去逼近真实的分布。Sigma 点就像是精心布置的采样点它们经过非线性函数后携带的信息比单点线性化丰富得多。2.3 两种算法流程的统一结构EKF 和 UKF 表面差很多但骨子里都是同一个贝叶斯递推框架预测阶段用状态方程推先验更新阶段用量测方程把先验修正成后验。区别只在于非线性怎么处理。我在代码里把两个算法都写成了统一的函数接口输入是上一时刻的状态、协方差、当前量测、输入量和系统结构体 sys输出是当前时刻的状态和协方差估计。系统结构体里封装了 f、h、Q、R、雅可比函数等所有依赖。主程序里切换 EKF 和 UKF 只需要改一行函数名对比实验做起来非常舒服。下面这个表是我做完整轮测试后整理的可以直接拿去用对比项EKFUKF非线性处理方式一阶泰勒展开线性化无迹变换逼近状态分布是否需要雅可比矩阵需要手推易错不需要理论精度一阶精度高斯分布下三阶精度计算量单次预测加更新2n1 次传播随维度线性增长实现难度模型简单时低复杂时雅可比繁琐流程固定机械性强典型失效模式线性化误差累积导致发散协方差非正定或参数选择不当3. Matlab 代码实现与参数整定实战3.1 代码架构与模型定义项目的代码我按模块拆好了目录结构大致这样main.m 主脚本 model/omib_model.m 单机无穷大系统模型 filters/ekf_predict_update.m EKF单步滤波 filters/ukf_predict_update.m UKF单步滤波 utils/generate_measurements.m 量测生成 utils/rmse_calc.m RMSE计算 scripts/plot_results.m 绘图脚本模型定义我放在一个函数里返回系统结构体。以二阶发电机模型为例function sys omib_model() % 系统参数 sys.w0 2*pi*50; % 同步角速度 rad/s sys.H 5.0; % 惯性时间常数 s sys.D 1.0; % 阻尼系数 pu sys.Pmax 1.5; % 最大电磁功率 pu sys.n 2; % 状态维度 sys.m 2; % 量测维度 sys.dt 0.01; % 采样步长 s % 状态转移函数 f(x,u,dt) % 状态 x [delta; domega] sys.f (x, u, dt) [x(1) sys.w0 * x(2) * dt; x(2) (u(1) - sys.Pmax*sin(x(1)) - sys.D*x(2)) / (2*sys.H) * dt]; % 量测函数 h(x,u) sys.h (x, u) [x(1); sys.Pmax*sin(x(1))]; % EKF 用到的雅可比矩阵 sys.F_jacobian (x, u, dt) [1, sys.w0*dt; -sys.Pmax*cos(x(1))*dt/(2*sys.H), 1 - sys.D*dt/(2*sys.H)]; sys.H_jacobian (x, u) [1, 0; sys.Pmax*cos(x(1)), 0]; % 噪声协方差 sys.Q diag([1e-5, 1e-5]); % 过程噪声 sys.R diag([5e-5, 1e-4]); % 量测噪声 end这里把状态转移函数和量测函数都写成匿名函数后面的滤波算法只跟 sys 打交道不关心具体模型。换模型时只需要改这个文件。3.2 EKF 主循环代码与实现要点EKF 的单步滤波函数不长核心就十几行function [x_est, P_est] ekf_predict_update(x_prev, P_prev, z, u, sys) % 预测 x_pred sys.f(x_prev, u, sys.dt); F sys.F_jacobian(x_prev, u, sys.dt); P_pred F * P_prev * F sys.Q; % 更新 h_pred sys.h(x_pred, u); H sys.H_jacobian(x_pred, u); S H * P_pred * H sys.R; K P_pred * H / S; % 用右除代替 inv数值稳定性更好 x_est x_pred K * (z - h_pred); P_est (eye(sys.n) - K * H) * P_pred; end这里有一个细节我要多说一句雅可比矩阵在预测阶段是在 x_prev 处求导在更新阶段是在 x_pred 处求导两处的线性化点不一样。我在最初写的版本里图省事用了同一个点结果滤波轨迹在扰动后出现了明显的滞后偏差排查了好久才注意到这个细节。另一个细节是 K 的计算我用的是P_pred * H / S这是 Matlab 里的右除等价于乘以 S 的逆矩阵但数值上更稳定、效率也更高。除非是为了教学演示否则不建议写inv(S)。3.3 UKF 主循环代码与实现要点UKF 的单步滤波函数代码量比 EKF 大但是结构非常固定function [x_est, P_est] ukf_predict_update(x_prev, P_prev, z, u, sys) n sys.n; alpha 0.1; beta 2; kappa 0; lambda alpha^2 * (n kappa) - n; % 生成 Sigma 点 sqrtP chol((n lambda) * P_prev, lower); chi [x_prev, x_prev sqrtP, x_prev - sqrtP]; % n x (2n1) % 权重 Wm zeros(1, 2*n1); Wc zeros(1, 2*n1); Wm(1) lambda / (n lambda); Wc(1) lambda / (n lambda) (1 - alpha^2 beta); Wm(2:end) 1 / (2 * (n lambda)); Wc(2:end) 1 / (2 * (n lambda)); % 预测Sigma 点通过状态方程传播 chi_pred zeros(n, 2*n1); for i 1:2*n1 chi_pred(:,i) sys.f(chi(:,i), u, sys.dt); end x_pred chi_pred * Wm; P_pred zeros(n, n); for i 1:2*n1 dx chi_pred(:,i) - x_pred; P_pred P_pred Wc(i) * (dx * dx); end P_pred P_pred sys.Q; % 更新预测 Sigma 点通过量测方程传播 z_pred_pts zeros(sys.m, 2*n1); for i 1:2*n1 z_pred_pts(:,i) sys.h(chi_pred(:,i), u); end z_pred z_pred_pts * Wm; % 计算协方差和增益 Pz zeros(sys.m, sys.m); Pxz zeros(n, sys.m); for i 1:2*n1 dz z_pred_pts(:,i) - z_pred; dx chi_pred(:,i) - x_pred; Pz Pz Wc(i) * (dz * dz); Pxz Pxz Wc(i) * (dx * dz); end Pz Pz sys.R; K Pxz / Pz; % 状态更新 x_est x_pred K * (z - z_pred); P_est P_pred - K * Pz * K; end写 UKF 时最需要注意的是chol分解前的矩阵必须是正定的。我遇到过协方差矩阵因为数值误差变得轻微不对称chol 直接报错的情况。解决办法是先对称化再尝试分解P_sym (P_prev P_prev) / 2; sqrtP chol((n lambda) * P_sym, lower);如果这样还是报错说明滤波器已经处于崩溃边缘需要回到 4.2 节的排查流程了。3.4 仿真场景搭建与参数初始化仿真场景我建议从单机无穷大系统开始状态维度低物理概念清晰调参直观算法跑顺了再往多机系统扩展也不迟。仿真时长设为 10 秒步长 0.01 秒共 1000 个滤波周期。扰动设置成机械功率在 t1 秒时从 0.8 pu 阶跃到 1.2 pu这样功角会经历明显的摆动过程正好考验滤波器的动态跟踪能力。量测序列的生成逻辑先用确定性模型算出真实状态轨迹再在量测方程的输出上叠加高斯噪声。为了让 EKF 和 UKF 的对比公平生成噪声前固定随机种子rng(42);滤波初值设置上有一点要提醒初始状态不能直接给真实值否则滤波器一开始就“开挂”了没法评估收敛性能。我是这么设的x_true_0 [asin(0.8/1.5); 0]; % 真实初始状态 x_est_0 x_true_0 [0.02; 0.005]; % 滤波初值稍微偏离真实值 P_0 diag([1e-2, 1e-3]); % 初始协方差反映对初值的不确定程度过程噪声 Q 和量测噪声 R 的选择是整场调试里最费时间的环节。我的经验是Q 反映的是模型误差和输入扰动的强度R 反映的是量测设备精度。PMU 的电压相量精度通常在 0.1% 左右功角换算后的误差约 0.01 弧度量级所以 R 的对角元素按物理精度去估计是靠谱的Q 反而更难定只能根据“滤波器会不会跟踪太慢”和“会不会噪声太大”来回试。最终我用的是 Q diag([1e-5, 1e-5])、R diag([5e-5, 1e-4])这个组合在多种扰动工况下都表现得比较稳。3.5 结果评估与对比指标评估指标我用的是 RMSE对每个状态分量分别计算function rmse rmse_calc(x_true, x_est) diff x_true - x_est; rmse sqrt(mean(diff.^2, 2)); end最后画三张图第一张是真值、EKF 估计值、UKF 估计值的功角轨迹对比第二张是角速度轨迹对比第三张是两个滤波器的误差曲线。误差曲线比轨迹更能看出问题——误差如果不是围绕零附近波动而是出现单向漂移基本可以断定滤波器有系统性偏差。从我的测试结果看在噪声不太大的工况下EKF 和 UKF 的精度差距并不悬殊UKF 在扰动后的前几百毫秒内跟踪更紧稳态阶段的 RMSE 则比较接近。但如果量测噪声加大或者模型参数失配UKF 的优势会越来越明显。这也符合理论线性化误差在强非线性或者大噪声场景下被放大得更厉害。4. 常见问题与排查技巧实录4.1 滤波器发散最常见的三个原因滤波器发散是最折磨人的问题症状一般是估计误差越来越大或者估计值直接飘走。我总结下来最常见的三个原因按出现频率排序过程噪声 Q 和量测噪声 R 的比例严重失衡。Q 太小滤波器会过分相信模型模型一有偏差预测就跟不上量测又修正不动Q 太大则会把量测里的噪声也当信息吸收进来估计曲线毛刺密布。调试的时候我会先固定 R把 Q 从极小值开始逐渐增大观察误差曲线的变化趋势。初始协方差 P0 设置不合理。P0 反映的是“我对初值有多不信任”。P0 太大滤波初期会剧烈跳动P0 太小滤波器收敛速度会非常慢甚至一直没反应过来。一般取初始状态不确定度的平方量级比较合适。模型和量测的不一致。比如状态方程里用了简化功率表达式但量测生成时用了更精细的网络方程这会导致量测残差里包含确定性偏差滤波器再怎么调都消除不掉。遇到这种问题先把“真值加噪”的量测生成方式改成“从滤波器自己的模型推出来的量测”如果残差立刻变小那就是模型的锅不是算法的锅。4.2 chol 分解报错与数值不稳定UKF 里 chol 报 Matrix must be positive definite 是非常经典的报错。原因通常是协方差矩阵 P 失去了正定性可能是 PMU 量测的数值尺度差异过大也可能是连续几次更新后数值误差逐渐累积。我的处理流程是先检查 P 矩阵有没有明显不对称不对称就做对称化处理 然后看看有没有负的特征值有负特征值说明滤波器已经崩溃单靠对称化救不回来 回退到上一时刻的良好状态或者增大 Q 重新初始化让协方差矩阵“膨胀”回正定区域。还有一个经验如果在 UKF 里用到高阶模型比如 n10 以上建议改用平方根 UKF 或者加一个简单的协方差映射策略比如 P P 1e-12 * eye(n) 这种极小的正则化项能避免大部分数值问题。4.3 EKF 和 UKF 公平对比的三个注意点做对比实验时如果方法不严谨结论很容易失真。我总结出三个必须注意的点固定随机种子。两套滤波器必须在完全相同的量测序列上跑否则每一次随机噪声都不一样对比结果毫无意义。Matlab 里在生成量测前写一句 rng(42) 就能解决。用 Monte Carlo 多次仿真取统计指标。单次仿真结果受噪声实现影响很大一次 UKF 恰好比 EKF 好下一次可能就反过来了。我一般至少要跑 50 次每次重新生成噪声序列但保持随机种子族一致最后用平均 RMSE 和中位数 RMSE 一起报告。分别调优再对比。我看到一些论文里 UKF 胜出是因为 EKF 的雅可比写错了或者 Q 没调好这属于无效对比。公平的做法是先让两套滤波器分别在典型场景下达到自己最好的性能再统一比较。我自己的经验是在二阶模型这种线性化程度不高的例子里EKF 调好之后精度并不逊色太多所以不要为了突出 UKF 而故意“黑”EKF。4.4 调试滤波器的几个小技巧分享几个常规文档里不会写的调试技巧都是我实测下来效率最高的。看卡尔曼增益 K 的变化趋势。K 的值代表滤波器对量测的信任程度。如果 K 持续趋近于零说明量测几乎不起作用滤波器变成了纯模型预测此时要么 Q 太小要么 R 太大如果 K 异常大说明滤波器基本抛弃了模型完全跟着量测跑此时量测噪声很可能被低估了。打印预测残差 z - h(x_pred)。残差的均值如果明显不为零说明模型里有系统性偏差残差如果剧烈跳动则通常是噪声设置问题。这个残差序列是最早暴露问题的信号。用“极端参数测试”定位问题范围。把 R 设成极小值如果滤波结果仍然很差问题大概率在模型或算法实现把 Q 设成极小值如果跟踪仍然不好问题大概率在量测方程或状态转移函数。这个二分类测试能快速缩小排查范围。仿真步长不要太大。欧拉法离散在步长 0.01 秒下对发电机慢动态没有问题但如果把步长拉到 0.1 秒部分工况下离散误差会大到让滤波器出现虚假振荡。我建议滤波步长保持在 0.01 到 0.02 秒之间既符合 PMU 的采样频率也能保证离散精度。做完这个项目最深的感触是滤波公式本身并不复杂真正花时间的永远是模型构建和参数整定。EKF 看起来实现简单但雅可比矩阵手推错一处结果就全错UKF 实现起来代码量更大但流程机械、不容易错对复杂模型的扩展也更友好。如果要我给新手排一个优先级我会说先把物理模型写对再调 Q 和 R最后才纠结选 EKF 还是 UKF——因为选错滤波器最多精度差点模型错了再花哨的算法也白搭。最后再分享一个我自己还在用的小技巧把状态方程和量测函数写成独立模块并在内部加一个日志开关每次调用记录输入输出调试时直接看每个时刻的预测残差几乎能立刻定位到是模型问题还是滤波参数问题。这个习惯让后面的多机系统扩展少走了很多弯路。