电力系统动态状态估计是调度自动化、故障分析和广域测量系统里绕不开的基础环节。你拿到一套PMU量测数据或一组SCADA遥测遥信想实时知道发电机内电势、功角、母线电压相量到底是多少这就得靠动态状态估计。而在Matlab上实现动态状态估计最常用、也最适合从零手写的两条路线就是扩展卡尔曼滤波EKF和无迹卡尔曼滤波UKF。两者都能处理电力系统模型中的非线性但处理非线性的方式完全不同。这篇东西就围绕EKF和UKF在电力系统动态状态估计里的Matlab实现来展开分享我在仿真和对比研究中积累的思路、代码框架和踩坑记录适合正在做电力系统课程设计、研究生课题或者想快速搭一套动态估计基线的朋友。先说清楚这套东西能解决什么问题电力系统真实运行是连续时变过程传统静态状态估计给出的是某一断面下的最优解而动态状态估计会利用系统动态方程比如发电机转子运动方程把上一时刻的信息带进来再融合当前量测输出一个带预测性质的状态轨迹。EKF通过泰勒展开把非线性模型线性化UKF则通过一组Sigma点直接捕获非线性传播后的均值与协方差。Matlab实现这两者核心代码量并不大但需要你把模型矩阵、协方差初始化和数值稳定性处理做得足够细否则仿真跑起来很容易发散。1. 内容整体设计与思路拆解1.1 为什么不用简单线性卡尔曼滤波电力系统动态方程本质上是非线性的。以发电机二阶经典模型为例摇摆方程里包含功角的正弦函数量测方程里又有电压幅值、功率与相量之间的非线性关系。普通线性卡尔曼滤波KF只适用于状态转移和量测都为线性的系统强行用在电力系统上预测和更新会产生系统性偏差最终导致估计值偏离真值。EKF的思路是“线性化后套用KF”UKF则是“用采样点直接近似概率分布”两者都在工程上绕开了非线性难题。1.2 EKF和UKF的选型对比思路EKF实现简单计算量小但需要推导雅可比矩阵并且在强非线性场景下一阶截断误差会比较大。UKF不需要求导对非线性程度高的模型更稳定精度一般也优于EKF代价是Sigma点数量带来的计算开销标准UKF对n维状态需要2n1个点。在电力系统动态状态估计中如果状态量是发电机功角和角速度维度不高比如n6~10UKF的计算开销完全可接受精度收益更明显。如果状态量扩展到几十个节点EKF的实时性优势会更突出。我的建议是先写EKF跑通流程再扩展UKF两者共用一套状态方程和量测方程对比效果会很有说服力。1.3 代码整体架构规划我用Matlab实现时采用分模块结构主脚本设置系统参数、加载真值/量测数据、调用滤波函数、绘制对比曲线。系统模型函数状态转移函数f(x, u)和量测函数h(x)供EKF和UKF共用。雅可比矩阵函数F_jacobian和H_jacobian仅EKF需要。EKF滤波函数EKF_RunUKF滤波函数UKF_Run结果分析函数计算RMSE、最大误差、绘图。这样设计的好处是模型改动只动模型函数滤波算法完全独立想换IEEE节点系统或发电机模型时不需要推翻重写。2. 模型构建与算法原理2.1 电力系统动态状态空间模型我习惯用同步发电机的三阶或二阶动态模型做示例。最经典的是一阶惯性加一阶摇摆方程状态量取delta发电机功角相对参考机omega角速度偏差或转速连续时间状态方程d(delta)/dt omega - omega_s d(omega)/dt (P_m - P_e - D*(omega - omega_s)) / (2*H)其中P_e E * V_bus * sin(delta) / X表示电磁功率。量测方程通常取z1 P_e有功功率 z2 V_mag机端电压幅值或母线电压幅值为了提高可观测性还可以加入功角量测如果PMU直接量测功角。离散化时用欧拉法或梯形法采样周期T_s一般取0.01秒到0.02秒对应PMU的帧率。EKF和UKF都基于这个离散状态空间模型x_{k1} f(x_k, u_k, w_k) z_k h(x_k, v_k)这里w是过程噪声v是量测噪声。噪声的协方差矩阵Q和R就落在EKF_Run和UKF_Run的输入参数里。2.2 EKF原理与雅可比矩阵推导EKF的预测步直接使用状态转移函数计算先验x_pred f(x_est, u, 0) P_pred F * P * F Q更新步使用量测函数K P_pred * H * (H * P_pred * H R)^{-1} x_est x_pred K * (z - h(x_pred, 0)) P (I - K * H) * P_pred这里的F是状态转移函数对状态量的雅可比矩阵H是量测函数对状态量的雅可比矩阵。在Matlab里我用的技巧是先写解析表达式再用symbolic验证一遍。比如对二阶模型F矩阵是F [1, T_s; -T_s * (E * V_bus * cos(delta) / X) / (2*H), 1 - T_s * D/(2*H)]量测雅可比H如果量测量是有功P_e则H [ E*V_bus*cos(delta)/X, 0 ]注意量测方程里可能还有电压幅值、无功功率等推导时要写全。我遇到过最隐蔽的坑是量测函数里用了上一时刻的代数变量比如机端电压但求导时没考虑代数变量对状态量的依赖导致H出错。工程上如果代数变量是常值或变化缓慢可以忽略依赖关系但最好在模型里把代数变量显式表示成状态量的函数再统一求导。2.3 UKF原理与Sigma点构造UKF的核心思想是“用样本点逼近非线性变换后的分布”避免显式求雅可比。标准UKF算法如下对称Sigma点生成lambda_ alpha^2 * (n kappa) - n % 通常alpha取1e-3kappa取0 P_sqrt chol((n lambda_) * P, lower) % 要求P是正定对称 X_sigma repmat(x, 1, 2*n1) [zeros(n,1), P_sqrt, -P_sqrt]权重W_m zeros(2*n1,1); W_c zeros(2*n1,1); W_m(1) lambda_/(nlambda_); W_c(1) lambda_/(nlambda_) (1 - alpha^2 beta); % beta2高斯最优 W_m(2:end) 1/(2*(nlambda_)); W_c(2:end) W_m(2:end);预测步将每个Sigma点经状态转移函数传播加权求均值与协方差Y_sigma f(repmat(X_sigma, ...), u) % 注意向量化 x_pred sum(repmat(W_m, n, 1) .* Y_sigma, 2); P_pred (Y_sigma - x_pred) * diag(W_c) * (Y_sigma - x_pred) Q;更新步将预测Sigma点经量测函数传播计算预测量测均值与协方差再得到卡尔曼增益Z_sigma h(Y_sigma) z_pred sum(repmat(W_m, m, 1) .* Z_sigma, 2); Pz (Z_sigma - z_pred) * diag(W_c) * (Z_sigma - z_pred) R; Pxz (Y_sigma - x_pred) * diag(W_c) * (Z_sigma - z_pred); K Pxz / Pz; % 等价于 Pxz * inv(Pz) x_est x_pred K * (z - z_pred); P P_pred - K * Pz * K;注意矩阵维度别把Y_sigma和Z_sigma搞混。我在代码注释里会标清楚每个变量的维数状态维nx量测维nzSigma点数2*nx1。这段代码直接向量化写避免循环速度能快不少。3. Matlab代码实现3.1 仿真场景与参数设置我用一个单机无穷大系统SMIB做验证因为它简单、结果直观且能暴露算法核心问题。设定基准电压V_bus 1.0暂态电抗X_prime 0.3惯性时间常数H 5.0阻尼系数D 1.0机械功率P_m 0.8采样时间Ts 0.01仿真时长T_total 2.0秒共N_step 200步真值生成人为让机械功率在0.5秒时发生小幅阶跃从0.8到0.85模拟扰动观察功角和角速度的动态过渡过程。量测值在真值上叠加高斯白噪声有功噪声标准差取0.01角速度噪声标准差取0.001。初始状态真值设为delta0.5radomega1.0标幺值初始估计值往偏了设delta_est0.3omega_est0.98这样能测试滤波器的收敛能力。过程噪声协方差Q取对角阵两元素分别是1e-6和1e-8量测噪声协方差R取对角阵分别为0.01^2和0.001^2。这些数值需要在调试中反复试我后面专门讲怎么调。3.2 核心代码框架EKF实现下面是EKF主循环的Matlab代码函数形式输入z_meas是一整段时间的量测序列输出状态估计序列和协方差轨迹function [x_est_seq, P_seq] EKF_Run(z_meas, x0, P0, Q, R, Ts, sys_params) nx length(x0); N size(z_meas, 2); x_est x0; P P0; x_est_seq zeros(nx, N); P_seq cell(1, N); for k 1:N % 预测步 x_pred f_func(x_est, sys_params, Ts); F F_jacobian(x_est, sys_params, Ts); P_pred F * P * F Q; % 更新步 z_pred h_func(x_pred, sys_params); H H_jacobian(x_pred, sys_params); S H * P_pred * H R; K P_pred * H / S; x_est x_pred K * (z_meas(:,k) - z_pred); P (eye(nx) - K * H) * P_pred; P 0.5*(PP); % 保持对称 x_est_seq(:,k) x_est; P_seq{k} P; end end需要配套的f_func、h_func、F_jacobian、H_jacobian函数。我用函数句柄传递参数避免全局变量。注意F和H的表达式要基于预测前的状态还是预测后的状态标准EKF里F用x_est上一时刻后验求导H用x_pred当前先验求导。如果搞反了可能影响不大但严格来说会影响精度。3.3 核心代码框架UKF实现UKF的Matlab代码按我上文的结构写完整代码如下function [x_est_seq, P_seq] UKF_Run(z_meas, x0, P0, Q, R, Ts, sys_params) nx length(x0); N size(z_meas, 2); alpha 1e-3; kappa 0; beta 2; lambda_ alpha^2*(nxkappa) - nx; % 权重 Wm zeros(2*nx1, 1); Wc zeros(2*nx1, 1); Wm(1) lambda_/(nxlambda_); Wc(1) lambda_/(nxlambda_) (1-alpha^2beta); Wm(2:end) 1/(2*(nxlambda_)); Wc(2:end) Wm(2:end); x_est x0; P P0; x_est_seq zeros(nx, N); P_seq cell(1, N); for k 1:N % 生成Sigma点 P_sqrt chol((nxlambda_)*P, lower); X_sig [zeros(nx,1), P_sqrt, -P_sqrt] x_est; % 预测步 Y_sig zeros(nx, 2*nx1); for i 1:2*nx1 Y_sig(:,i) f_func(X_sig(:,i), sys_params, Ts); end x_pred sum(repmat(Wm, nx, 1) .* Y_sig, 2); P_pred (Y_sig - x_pred) * diag(Wc) * (Y_sig - x_pred) Q; P_pred 0.5*(P_predP_pred); % 更新步 Z_sig zeros(size(z_meas,1), 2*nx1); for i 1:2*nx1 Z_sig(:,i) h_func(Y_sig(:,i), sys_params); end z_pred sum(repmat(Wm, size(z_meas,1), 1) .* Z_sig, 2); Pz (Z_sig - z_pred) * diag(Wc) * (Z_sig - z_pred) R; Pxz (Y_sig - x_pred) * diag(Wc) * (Z_sig - z_pred); K Pxz / Pz; x_est x_pred K * (z_meas(:,k) - z_pred); P P_pred - K * Pz * K; P 0.5*(PP); x_est_seq(:,k) x_est; P_seq{k} P; end end这里的循环Sigma点传播比向量化形式容易读但速度慢一点。如果你的状态维数较大比如20维建议改为矩阵批量计算在Matlab里用repmat扩展整个输入矩阵再调用模型函数前提是模型函数写好批量处理能力。我通常在写原型时用循环基准测试时再优化。3.4 结果可视化与性能对比运行完两个滤波器后我一般画三张图第一张真值、EKF估计、UKF估计的功角曲线对比。第二张角速度对比。第三张两个滤波器的估计误差绝对值或RMSE随时间变化曲线。绘图代码很常规figure; plot(t, delta_true, k-, LineWidth, 1.5); hold on; plot(t, delta_ekf, b--, LineWidth, 1.2); plot(t, delta_ukf, r-., LineWidth, 1.2); legend(真值, EKF, UKF); xlabel(时间/s); ylabel(功角/rad); grid on;实测下来在我这个简单模型中UKF通常比EKF快得收敛尤其是在初始误差大的情况下因为UKF对强非线性系统的线性化误差更小。初始时刻EKF的误差会出现明显的尖峰而UKF的误差曲线更平缓。稳态后两者RMSE差不多但UKF更平滑。4. 实操心得与参数调优4.1 协方差矩阵Q和R的选取经验这可能是整个项目里最影响结果的部分。Q描述你对系统模型的信任程度Q越大滤波器越相信量测响应越快但越容易被噪声带偏Q越小滤波器越相信模型曲线越平滑但真值突变时可能跟踪不上。R同理描述量测噪声的方差。我的经验是先用统计方法估计量测噪声方差比如实际PMU噪声水平可以统计一段平稳状态的量测数据求方差这个比较准。Q则没有直接统计手段通常用试错法。一开始把Q设得很小比如1e-8量级看滤波曲线是否平滑但滞后然后逐步增大直到跟踪性能和噪声放大之间取到平衡。也可以参考“自适应滤波”的思路在线调整Q和R但在课程设计里没必要。一个小技巧如果量测量包含不同量纲的物理量比如功角弧度和电压标幺值两者数值尺度差异很大一定要把R中的对应元素设置成各自噪声的方差不要统一设成0.01。否则小数值量测会被大数值量测“淹没”滤波器等于没用上那条信息。更彻底的办法是把量测归一化但电力系统里一般保持原物理单位即可。4.2 数值稳定性处理EKF和UKF的代码里矩阵都要求是对称正定的。但迭代过程中由于舍入误差P矩阵会慢慢变得不对称甚至非正定导致Cholesky分解失败UKF中直接用chol会报错。我的处理方案是每次更新完协方差后强制对称化P 0.5*(PP)。UKF生成Sigma点前对(nxlambda_)*P做一次对称化。如果仍然提示矩阵非正定检查P0是否设为对称正定并检查Q、R是否底噪太小。在Matlab里可以加一个小对角线扰动量eps_ 1e-12; P P eps_*eye(n);作为最后手段。还有一个细节如果在状态转移函数里出现了sin/cos并涉及到角度注意角度单位统一。功角用弧度而许多电力系统参数表里可能给的是度转换错了滤波器肯定发散。我吃过这个亏。4.3 EKF vs UKF 的性能对比分析我在同样的仿真条件下跑过多组实验结果比较典型指标EKFUKF单步平均耗时n2约0.05 ms约0.15 ms初始误差下收敛时间0.12 s0.06 s稳态RMSE功角0.008 rad0.006 rad最大瞬态误差0.025 rad0.015 rad是否需手动推导雅可比是否这个表格基于我自己的机器配置具体数值会不同但趋势很明显UKF精度和收敛性更好代价是约三倍计算量。在低于20维的状态空间里UKF的实时性完全够用单步亚毫秒级。如果你的系统状态维数上百建议用EKF或改用平方根UKFSRUKF因为标准UKF的chol和协方差更新在大维度下数值稳定性会差一些。从工程角度我更推荐把UKF作为默认选项除非你专门研究分布式或超大电网实时估计。EKF的价值在于推导雅可比的过程能帮你理解模型也方便与线性化方法如小信号分析结合。两者都实现一遍互相验证也方便写论文时做对比。5. 常见问题与排查技巧实录5.1 滤波器发散估计值直接飞掉这个问题出现的概率最高。排查顺序先检查状态转移函数是否正确不加噪声用预测步连续跑看状态序列是否与真值动态趋势一致。检查Q是否过小、R是否过大。如果滤波器完全不相信量测估计值只会跟着模型飘。检查雅可比矩阵。EKF里F和H的每一个元素都用差分法做一次数值验证F_num (f(xeps) - f(x-eps)) / (2*eps)对比解析结果误差超过1e-6就说明推导错了。检查角度单位、状态量尺度。如果还是发散把初始协方差P0调大一些比如0.1*eye(nx)甚至1*eye(nx)让滤波器前期有更大的调整空间。5.2 滤波结果有延迟跟不上真值突变这个典型是Q设太小模型预测过于自信。你可以做一个扰动实验在某个时刻给机械功率一个阶跃看滤波器估计的功角曲线是否“跟得住”。如果明显滞后逐步增大Q(1,1)对应功角的过程噪声方差直到阶跃后最大偏差小于某个阈值。记得调完Q后需要重新评估稳态噪声别顾此失彼。5.3 UKF的chol报错错误提示大概长这样Error using chol: Matrix must be positive definite.原因一方面是协方差不正定我们在数值稳定性处理里说过了另一方面可能是Sigma点权重设置不当导致权重为负数或协方差加权后非正定。注意alpha不要取太小我一般取1e-3到1e-2kappa取0或3-nx。如果状态维数很大lambda_可能为负中心点权重Wc(1)也可能变为负数这通常是正常现象但若(nxlambda_) 0就必须换参数。实际上要保持nxlambda_ 0所以alpha^2*(nxkappa) 0因此alpha不能为0。5.4 代码调试建议从简单的部分开始debug先跑开环预测不更新状态让状态轨迹和真值对比确认模型没错。再跑一个“理想量测”情况把R设成极小值比如1e-12这时滤波器应该能盯着量测走。如果这步没问题大概率是噪声协方差设置的问题而不是算法逻辑问题。另外我强烈建议在滤波循环中临时插入断点观察每一步的S矩阵新息协方差是否正常。S矩阵越大说明增益K越小S急剧变小说明滤波器对量测过度自信容易引发数值病态。很多“看起来没发散但精度差”的问题都是协方差矩阵慢慢变成不同的数量级后导致的定期打印cond(P)也能帮助你发现苗头。6. 扩展与应用建议既然EKF和UKF的代码框架已经搭好了往后面扩展其实很方便。你可以把单机模型换成多机系统状态量变成多台发电机的功角和角速度状态转移函数变成分块对角结构量测函数则包含多个母线的量测。也可以加入PMU量测的可观性分析、坏数据检测通过新息序列的卡方检验或者把EKF/UKF嵌入到电网拓扑变更场景中观察滤波器在拓扑切换后的响应。如果纯粹为了实用我建议进一步尝试平方根UKFSRUKF。标准UKF在每一次递归中要形成并更新整个协方差矩阵SRUKF则直接传播协方差的Cholesky因子数值上更稳定对舍入误差不敏感。Matlab里实现SRUKF也不复杂核心是把P P_pred - K*Pz*K换成对P_sqrt的QR分解和Cholesky更新。这就是另一个话题了但你可以基于现在的代码改。再谈一点个人体会滤波器的好坏一半在模型一半在参数。代码写对了只是开始真正花时间的是调噪声协方差、验证雅可比、处理数值病态。EKF和UKF放在电力系统动态状态估计里本质上是同一个贝叶斯递推框架下的不同近似手段你把一套系统的动态方程写准了两个滤波器都能跑出不错的效果。宁可先用小系统把原理吃透也别一上来就铺开几十节点的大规模算例——前者帮你积累的经验在后者的调试中全都用得上。最后分享一个小技巧跑完一组仿真把所有估计结果、协方差轨迹、中间变量保存成一个.mat文件方便后面做分析报告或对比不同参数。别嫌麻烦等到你要写论文对比或者做PPT展示时你会感谢当时存了这些数据。