最近在帮几个师弟调电力系统动态状态估计的仿真代码发现大家卡住的位置都差不多——模型都能写出来卡尔曼增益公式也背得熟但真到Matlab里把EKF和UKF跑起来各种发散、矩阵奇异、量测噪声调不动的问题全冒出来了。这篇文章就围绕EKF和UKF在电力系统动态状态估计中的Matlab实现把数学模型、滤波原理、代码细节和调试经验完整梳理一遍。适合电力系统方向的研究生、做广域测量系统动态分析或新能源并网研究的工程师参考也适合刚接触卡尔曼滤波、想用Matlab做非线性状态估计的同学直接抄作业。1. 项目概述与总体设计思路1.1 动态状态估计到底在解决什么问题传统电力系统状态估计用的是加权最小二乘WLS那一套静态方法基于SCADA量测每隔几秒甚至几分钟出一个断面。但现代电力系统里新能源出力波动、负荷快速变化、功角振荡这些动态过程越来越频繁静态断面已经跟不上实际运行节奏。PMU同步相量测量单元的普及改变了这个局面它能以几十到几百赫兹的频率输出带时标的相量数据这就让动态状态估计DSE变得可行。动态状态估计的本质是把电力系统的动态模型写成状态空间形式用滤波器从含噪声的PMU量测里实时递推状态量。这里的“状态”通常指发电机功角、角速度、暂态电动势等。它的输出是“每一时刻”的状态估计值而不是某一个断面的静态解所以能跟踪系统动态轨迹为后续的稳定性分析、保护控制提供连续的输入。我用的是最常用的发电机二阶模型状态向量取功角δ和角速度ω。量测取PMU可提供的功角、角速度以及电磁功率三者都叠加独立的高斯白噪声。这样既保留了电力系统的非线性特征又不至于让模型复杂到难以调试适合作为EKF和UKF的对比研究对象。1.2 为什么偏偏选EKF和UKF非线性滤波的算法其实不少粒子滤波、H∞滤波、滑模观测器都有研究应用。但EKF和UKF在电力系统DSE里的地位很特殊主要原因是它在精度和计算量之间取得了最实用的平衡。EKF的思想最直白在估计点附近把非线性函数做一阶泰勒展开用雅可比矩阵近似线性化然后套用经典卡尔曼滤波框架。它的实现门槛低计算量小在系统运行点变化平缓时表现很好。缺点也很明确一阶线性化会丢掉高阶信息遇到强非线性场景比如故障后的功角大摆误差会变大甚至发散。UKF是在这个思路上做了升级。它不显式求雅可比而是取一组确定性采样点sigma点让这些点通过非线性函数传播再用加权统计量近似状态的均值和协方差。相当于用“抽样逼近”代替“线性化逼近”理论上能捕获到二阶甚至三阶的统计特性。代价是计算量比EKF大一些但不需要推导雅可比矩阵这在大规模系统和复杂量测方程里是很大的优势。另外还有一个非常实际的原因Matlab里实现这两种算法都不需要额外工具箱核心代码量能控制在三四百行以内非常适合教学和科研验证。1.3 我对这个项目的工程假设在做仿真前我先把假设条件定死省得后面调试的时候变量太多系统采用标幺值频率基准50Hzω02π×50。发电机用经典二阶模型E、Xd、Vt均视为恒定常数。过程噪声和量测噪声均为零均值高斯白噪声且彼此不相关。采样周期Ts固定为0.01s对应100Hz的PMU上报频率。滤波器的初始状态在真实值附近加一个较小的偏差初始协方差按经验设定。这些假设和实际PMU工程场景基本能对上又不会引入太多参数干扰算法本身的对比。仿真里我让系统在0.5s的时刻发生一次机械功率阶跃扰动用来检验动态跟踪能力。2. 两种滤波算法的数学原理与递推流程2.1 EKF的线性化思想与核心递推式假设系统为x_{k} f(x_{k-1}, u_{k-1}) w_{k-1}z_{k} h(x_{k}) v_{k}其中w和v分别是过程噪声和量测噪声。EKF在每一步先计算雅可比矩阵F_{k-1} ∂f/∂x|{x̂{k-1}}H_{k} ∂h/∂x|{x̂{k|k-1}}然后走标准预测-更新两步。预测步x̂_{k|k-1} f(x̂_{k-1|k-1}, u_{k-1})P_{k|k-1} F_{k-1} P_{k-1|k-1} F_{k-1}^T Q更新步K_{k} P_{k|k-1} H_{k}^T (H_{k} P_{k|k-1} H_{k}^T R)^{-1}x̂_{k|k} x̂_{k|k-1} K_{k} (z_{k} - h(x̂_{k|k-1}))P_{k|k} (I - K_{k} H_{k}) P_{k|k-1}这里最考验人的就是雅可比矩阵。解析推导要小心求导数值求则要注意步长。我的实践结论是对于二阶发电机模型这种维数不高、非线性度适中的系统数值雅可比完全够用解析表达反而容易因为推导错误引入低级bug。2.2 UKF的sigma点传播思路UKF回避了雅可比它的核心是构造2n1个sigma点n是状态维数。设状态均值x̄、协方差P使用比例对称采样λ α²(n κ) - nS chol((n λ)P)取Cholesky分解的下三角矩阵sigma点集合为χ₀ x̄χ_i x̄ S_iχ_{in} x̄ - S_i其中S_i是S的第i列。权重按下面方式分配W₀^m λ/(n λ)W₀^c λ/(n λ) (1 - α² β)其余点的权值均为1/[2(n λ)]α控制sigma点离均值多远β在高斯分布时取2最优κ是比例参数通常取0或3-n。sigma点生成之后每个点都通过状态方程和量测方程传播然后带权组合出预测均值、协方差和互协方差x̂_{k|k-1} Σ W_i^m χ_{i,k|k-1}P_{k|k-1} Σ W_i^c (χ_{i,k|k-1} - x̂)(...)ᵀ Qẑ_{k|k-1} Σ W_i^m ζ_{i,k|k-1}P_{zz} Σ W_i^c (ζ_i - ẑ)(...)ᵀ RP_{xz} Σ W_i^c (χ_i - x̂)(ζ_i - ẑ)ᵀ然后增益K P_{xz} P_{zz}^{-1}状态和协方差更新。从实现角度看UKF比EKF多的工作主要集中在sigma点生成和两组权重数组管理上核心计算量约是EKF的三倍左右。但换来的是不需要任何导数信息在量测函数很复杂的时候优势极其明显。2.3 两者在电力系统DSE场景下的对比维度EKFUKF非线性处理一阶泰勒展开sigma点统计逼近雅可比矩阵需要不需要计算量小约为EKF的2-3倍强非线性下精度可能恶化通常更好实现难度低但要仔细求雅可比中要处理数稳细节典型适用弱非线性、实时性要求高的场景强非线性、量测复杂场景我在这次仿真里的实测结论是单机系统正常工况下两者精度差异不大EKF甚至因为数值噪声更小显得更稳。但把扰动加大、功角摆动幅度拉高后UKF的跟踪优势就开始体现出来。如果以后扩展到大电网多机系统我会优先考虑UKF。3. 电力系统动态模型构建与离散化3.1 发电机二阶状态方程的建立单机无穷大系统的经典二阶模型写成连续形式dδ/dt ω - ω0dω/dt (Pm - Pe - D(ω - ω0)) / M其中M 2H/ω0H为惯性时间常数D为阻尼系数Pm为机械功率Pe为电磁功率。在经典模型中Pe由暂态电势E、机端电压Vt和直轴暂态电抗Xd决定Pe (E Vt / Xd) sin δ虽然这个模型在电力系统教科书里很基础但作为动态状态估计的验证平台非常合适因为它保留了sin(δ)这个关键非线性项能真实检验滤波器对非线性的处理能力。我用Matlab把这些参数放到结构体里params.E 1.05; % 暂态电动势 pu params.Vt 1.0; % 机端电压 pu params.Xd 0.3; % 直轴暂态电抗 pu params.Pm 0.8; % 机械功率 pu params.H 5; % 惯性时间常数 s params.D 2; % 阻尼系数 pu params.Ts 0.01; % 采样周期 s params.omega0 2*pi*50; params.M 2*params.H/params.omega0;这里有一个常见单位坑功角δ是弧度角速度ω一般也用带ω0的量纲rad/s。为了和标幺体系兼容我让所有计算按“pu”统一进行最后画图时再换算成有名值或角度值。3.2 量测方程的设计我设计的量测向量包含三个量z [δ, ω, Pe]ᵀ其中前两个是状态量本身的直接量测第三个是通过电气量间接计算得到的电磁功率。量测方程写成δ_m δ v₁ω_m ω v₂P_em (E Vt / Xd) sin δ v₃把第三项加进去量测方程就变成非线性的了EKF在更新步就必须要算h对x的雅可比这个设计能更加真实地考验UKF的值逼近能力。量测噪声标准差我设成pa.delta0.005rad、pa.omega0.005pu、pa.Pe0.02pu也就是量测里功角和转速的噪声很小电磁功率噪声稍大。3.3 真实轨迹与仿真量测的生成做滤波仿真之前要先造一套“真实值”用来评价算法好坏。我的做法是直接用状态方程对真实初始状态做开环递推把结果当成系统真实轨迹再在量测方程上叠加高斯噪声形成PMU量测。x_true zeros(2, N); z_meas zeros(3, N); x_true(:,1) [0.2; omega0]; % 初始功角0.2rad for k 1:N-1 x_true(:, k1) f_func(x_true(:,k), params); end % 量测生成 for k 1:N h_val h_func(x_true(:,k), params); z_meas(:,k) h_val sigma_v .* randn(3,1); end这样生成的量测序列同时包含动态信息和噪声后面跑滤波器时用z_meas作为输入再把状态估计值和x_true做对比计算均方根误差RMSE这类量化指标。4. Matlab代码实现的关键环节4.1 状态方程与量测方程的函数封装我习惯把非线性函数写成独立的M函数文件这样EKF和UKF可以共用同一套模型避免模型不一致导致的对比失真。function x_next f_func(x, params) delta x(1); omega x(2); Pe params.E * params.Vt / params.Xd * sin(delta); x_next [ delta (omega - params.omega0) * params.Ts; omega (params.Pm - Pe - params.D*(omega - params.omega0))/params.M * params.Ts ]; end function z h_func(x, params) delta x(1); omega x(2); Pe params.E * params.Vt / params.Xd * sin(delta); z [delta; omega; Pe]; end这里有个容易被忽视的细节状态方程的欧拉离散需要在delta的递推中使用当前的omega而不是使用上一步的omega。写成上面的形式就是正确的时序关系。4.2 EKF的雅可比矩阵数值计算解析雅可比需要手推偏导在状态方程里又要对sin(δ)求导又要对ω变量求导虽然这个模型不算复杂但一旦以后扩展成四阶、六阶发电机模型解析推导的工作量会直线上升。所以我在主程序里统一用中心差分法求数值雅可比代码十分简洁function J numerical_jacobian(fun, x, params, dx) if nargin 4 dx 1e-6; end n length(x); f0 fun(x, params); m length(f0); J zeros(m, n); for i 1:n xp x; xp(i) x(i) dx; xm x; xm(i) x(i) - dx; J(:, i) (fun(xp, params) - fun(xm, params)) / (2*dx); end end中心差分比单边差分精度高步长dx我控制在1e-6到1e-7之间。dx太大会让线性化误差占主导太小又容易受浮点舍入噪声影响。EKF主循环代码x_ekf(:,1) x_init; P_ekf P_init; for k 1:N-1 % 预测 x_pred f_func(x_ekf(:,k), params); F_k numerical_jacobian(f_func, x_ekf(:,k), params); P_pred F_k * P_ekf * F_k Q; % 更新 H_k numerical_jacobian(h_func, x_pred, params); z_pred h_func(x_pred, params); S_k H_k * P_pred * H_k R; K_k P_pred * H_k / S_k; x_ekf(:,k1) x_pred K_k * (z_meas(:,k1) - z_pred); P_ekf (eye(2) - K_k * H_k) * P_pred; end这里我用的是K P_pred * H / S_kMatlab里“/”等价于右除乘以S_k的逆数值上更推荐用右除而不是显式求逆能有效避免条件数很大的情况。4.3 UKF的sigma点生成与主循环实现UKF实现里最容易出错的地方是Cholesky分解的方向。Matlab默认chol返回的是上三角矩阵如果直接用它去加减列向量列的对应关系会乱。我建议加lower参数或者用chol(...)做转置。L chol((n lambda) * P_ukf, lower); X_sigma [x_ukf, x_ukf L, x_ukf - L]; % 2n1列每列是一个sigma点权重数组提前算好n 2; alpha 1e-2; beta 2; kappa 0; lambda alpha^2 * (n kappa) - n; Wm ones(2*n1, 1) / (2*(n lambda)); Wc Wm; Wm(1) lambda / (n lambda); Wc(1) lambda / (n lambda) (1 - alpha^2 beta);UKF主循环for k 1:N-1 % sigma点通过状态方程 X_sigma_pred zeros(n, 2*n1); for j 1:2*n1 X_sigma_pred(:,j) f_func(X_sigma(:,j), params); end x_pred X_sigma_pred * Wm; P_pred (X_sigma_pred - x_pred) * diag(Wc) * (X_sigma_pred - x_pred) Q; % sigma点通过量测方程 Z_sigma_pred zeros(3, 2*n1); for j 1:2*n1 Z_sigma_pred(:,j) h_func(X_sigma_pred(:,j), params); end z_pred Z_sigma_pred * Wm; Pzz (Z_sigma_pred - z_pred) * diag(Wc) * (Z_sigma_pred - z_pred) R; Pxz (X_sigma_pred - x_pred) * diag(Wc) * (Z_sigma_pred - z_pred); K Pxz / Pzz; x_ukf(:,k1) x_pred K * (z_meas(:,k1) - z_pred); P_ukf P_pred - K * Pzz * K; P_ukf (P_ukf P_ukf) / 2; % 强制对称 end这个实现用清晰的分布式写法牺牲了一点向量化效率但好处是逻辑透明初学者对照公式能一步步看懂。如果追求速度可以把sigma点循环改成矩阵运算不过对DSE仿真来说这点差距无所谓。4.4 结果可视化与评估指标仿真结束以后我画三张图功角跟踪曲线、角速度跟踪曲线、电磁功率量测与预测对比。计算两个指标功角估计的RMSE和角速度估计的RMSE。rmse_delta sqrt(mean((x_ekf(1,:) - x_true(1,:)).^2)); rmse_omega sqrt(mean((x_ekf(2,:) - x_true(2,:)).^2));RMSE别把整个仿真段都算进去。前几十个采样点滤波器还在收敛状态估计偏差很大会严重拉高RMSE。我会跳过前10%的暂态段再计算稳态段RMSE这样反映的才是滤波器真正的跟踪精度。5. 参数选取策略与仿真结果分析5.1 过程噪声Q和量测噪声R的调参心得噪声协方差矩阵是卡尔曼滤波里最敏感的参数也是很多人调起来最没头绪的地方。对于量测噪声R我的原则非常直接既然量测是真实PMU数据R就应该由传感器精度决定。我这里直接在生成量测时用已知噪声标准差来构造R等于给了滤波器一个“标准答案”。如果使用实际PMU数据最稳妥的办法是统计一段静态工况下量测的时间序列方差。对于过程噪声Q麻烦一些。Q描述的是状态方程未能反映的动态随机性包括模型误差、参数漂移、外部扰动等没法直接测量。我采取的做法是先固定R然后把Q当做一个标量q乘以单位阵去调Q q × [1, 0; 0, 1]从q1e-4开始依次试到q1e-8观察滤波曲线的平滑程度与延迟。q大了滤波器会过度信任量测跟着噪声走曲线毛刺多q小了滤波器过度信任模型真实轨迹发生突变时跟踪滞后明显。最后我在单机系统里取q1e-6功角和角速度的跟踪曲线既平滑又能跟上扰动。5.2 滤波器初始化状态估计的初值我取真实值的邻域进行偏离比如真实初始功角是0.2rad我设初始估计为0.25rad角速度初始估计为ω0 0.1人为制造一个启动误差这样能检验滤波器收敛能力。初始协方差P_init设为单位阵乘以0.1左右。P_init表示对初始状态的信任程度取得太大会让滤波器在最初阶段出现大幅修正取得太小则容易陷在初始误差里。我在实际调试中体会是P_init稍微给大一点点让滤波器“自动”收敛到正确状态比反复试调初值省心得多。5.3 EKF与UKF的仿真结果对比机械功率在0.5s时从0.8阶跃到0.9相当于给系统施加了一个阶跃扰动。从仿真结果来看EKF和UKF的功角估计最终都收敛到真实轨迹附近稳态RMSE差别不大。但在扰动发生后的前0.2sEKF的跟踪曲线出现了一个较为明显的振荡UKF则更平滑地跟上了过渡过程。这个结果完全符合理论预期——阶跃扰动让功角轨迹的局部非线性增强一阶线性化的EKF在这个时刻出现精度损失而sigma点传播的UKF对非线性分布拟合得更好。角速度估计上两种滤波器的差距相对小一些因为ω的动态方程在扰动后很快回到线性区域。这张表可以直观看到两种算法在相同参数下的精度数据模拟结果算法功角稳态RMSE (rad)角速度稳态RMSE (pu)单步耗时 (μs)EKF0.00850.00316.2UKF0.00680.002615.8单步耗时是Matlab里用tic/toc粗略测的绝对数值和机器配置关系很大但相对差距基本稳定UKF大约是EKF的2.5倍计算量换来了约20%-30%的精度提升。这个性价比在单机系统里不算突出但在量测更复杂、非线性更强的场合会明显放大。6. 常见问题与排查技巧实录6.1 滤波发散的典型场景我调试过程中遇到过三次发散三次的原因各不相同都很有代表性。第一次是因为初始协方差P_init设得太小。初始状态估计偏差明明不小P_init却只有1e-8的级别滤波器觉得自己“很确定”卡尔曼增益被压得很低状态估计长时间不更新误差越来越大看起来像发散。解法是把P_init放大到0.01以上让滤波器在前期敢于修正。第二次是数值雅可比步长dx取值太离谱。我一开始用了1e-4对二阶模型来说线性化误差偏大在强非线性段导致增益计算错误最终滤波轨迹飞出天际。把dx改成1e-6之后立刻恢复正常。第三次源于量测序列里的噪声设定与R不一致。我生成量测时噪声加倍了R还是按原来较小的噪声方差设置滤波器过度信任被污染的测量数据结果在扰动时刻前后反复震荡。这提醒我仿真里的R必须和实际量测噪声统计特性严格一致。6.2 协方差矩阵非对称与非正定问题EKF和UKF在递推过程中协方差矩阵P容易出现非对称甚至非正定的情况尤其在长时间递推或非线性较强时。我给出的处理办法很粗暴有效每次更新之后强制对称化P (P P) / 2如果协方差矩阵仍然出现非正定的疑难情况比如某些特征值变成负的多半是数值精度或者sigma点分布出了问题。我会检查cholesky分解能否成功如果chol报错就说明P确实不是正定矩阵需要回到Q和R的调参上找原因。还有个极为实用的细节是生成sigma点时用的Cholesky分解如果P有微弱非正定倾向就会直接报错。我见过不少同学在这里被卡住。一个稳妥的补丁方法是在分解前给P对角加一个极小的扰动P_ukf P_ukf 1e-9 * eye(n);这个微小的正定化操作不影响滤波精度但能让程序稳健运行。6.3 计算效率与批量仿真意识单次仿真跑几十万步时EKF与UKF的速度差距才会真正体现出来。如果要做蒙特卡洛实验或者参数扫描建议注意两点第一提前把Q、R、sigma点权重这些不随迭代变化的量计算好别放在循环里重复生成。第二把最内层的sigma点传播循环尽量向量化。我的经验是对n2的简单模型向量化能快5倍左右从代码可读性角度还是保持循环形式但一旦做批量仿真就切换成向量化版本。第三能预分配的大数组x_ekf、x_true、z_meas一定提前预分配否则Matlab在循环里反复扩容会让速度拖慢一个数量级。N 2000; x_ekf zeros(2, N); P_ekf P_init;这样的小改动对长仿真的提升极其明显。6.4 实用避坑清单做这类仿真踩过太多次坑整理成一张速查清单放在这里比我前面啰嗦的所有内容都实用功角单位坑画图时记得把弧度转角度不然看着像发散实际只是单位问题。量测噪声坑randn默认标准差为1生成特定噪声方差时记得乘上对应倍率。扰动时刻坑仿真中设置阶跃扰动时滤波器的预测值和量测值会出现短暂失配这是正常现象不要一看到波动就怀疑算法错误。数据段划分坑RMSE的计算一定要区分启动段和稳态段否则启动误差会掩盖真实精度差异。代码版本坑Matlab的chol函数在旧版本中默认返回上三角新版本支持lower参数。写代码时务必注意你的Matlab版本否则sigma点生成全是错的。我个人在实际调试中最深的一个体会是EKF和UKF的差距远没有网上很多文章渲染得那么大。在单机电力系统模型里只要Q和R调得合理EKF的表现其实非常能打。UKF的真正价值在于省去了雅可比推导当系统规模变大、量测方程变得复杂甚至出现不连续环节的时候UKF的工程便利性才真正凸显出来。所以刚上手做这个方向的同学我建议先把EKF老老实实写通理解了线性化和增益更新是怎么回事再切换到UKF你会发现后面这个无非就是“换了一种算概率分布的方式”架构上并没有本质障碍。