1. 立项思路动态状态估计在电力系统里到底解决什么问题1.1 动态状态估计 vs 静态状态估计差在哪儿先说说背景。做过电力系统状态估计的人都知道传统做法是加权最小二乘WLS那一套输入SCADA传来的量测值输出一个断面的电压幅值和相角。这个方案在稳态监控里没问题但它的本质是用一个时间断面的冗余量测去拟合一个状态解对时间的演化完全没有概念。系统一旦进入动态过程比如负荷突变、发电机甩负荷、线路故障后的暂态SCADA的刷新率本来就低秒级甚至分钟级静态估计的结果往往滞后一拍等算出来系统可能已经跑到另一个工况去了。动态状态估计的思路完全不同。它把系统自身的运动规律写进模型里——最常见的就是发电机转子运动方程——然后用卡尔曼滤波这类递推算法把模型预测和新到的量测融合起来。每一步估计都包含两个动作先用状态方程把上一时刻的状态推进到当前时刻再用当前时刻的量测去修正这个预测值。这样做的结果就是估计值有预测能力能跟踪动态过程而且天然带有平滑噪声的功能。我做这个课题的时候目标很明确在Matlab里实现两种非线性卡尔曼滤波器用于电力系统动态状态估计并且把它们的精度、稳定性和实现成本做一次横向对比。之所以同时做EKF和UKF就是因为这两者是处理非线性状态估计时最经典的两条路线——一个靠线性化打补丁一个靠采样点传播理解清楚这两条路线的差异后续做扩展卡尔曼滤波的改进型如无迹变换、粒子滤波都有了参照。1.2 算法选型为什么不直接上粒子滤波现在做状态估计可选的非线性滤波算法其实很多粒子滤波、容积卡尔曼CKF、甚至神经网络辅助的估计器都有。但在这个项目里选EKF和UKF是有实际考虑的。首先是模型特点。电力系统动态估计的核心非线性体现在电磁功率与功角之间的正弦关系这种非线性程度属于中等偏强。EKF的适用前提是系统在估计点附近近似可线性化只要扰动不是特别剧烈功角不会大范围跳变一阶泰勒展开的精度基本够用。而UKF不做线性化直接用Sigma点去传播状态分布对正弦非线性天然更友好理论上可以达到二阶精度。其次是实现成本。粒子滤波虽然能处理强非线性但需要大量粒子在Matlab里做几百个粒子的递推估计单步计算的耗时就会明显上升而我的仿真场景要求Ts取到10毫秒甚至更小计算效率必须考虑。EKF和UKF的计算量都可控UKF对n维状态只需要2n1个Sigma点在状态维数不高的电力系统模型中非常划算。最后是想做一个标准答案级别的对比。EKF的线性化误差在强非线性场景下会暴露UKF则在Sigma点参数选择不当时也会出现协方差非正定的问题。两套算法放在同一个仿真环境里用相同的真值轨迹、相同的噪声统计对比结果才有说服力。这个对比数据无论写论文还是做汇报都是很有分量的支撑。2. 数学模型和算法原理把滤波器和电力系统模型接上2.1 动态估计用的状态方程和量测方程先说状态方程。我采用经典的二阶发电机模型即转子摇摆方程swing equation。状态量取两维x1 δ发电机功角radx2 Δω角速度偏差p.u.基准为同步角速度连续时间模型为dδ/dt ω0 * Δω dΔω/dt (P_m - P_e - D * Δω) / (2H)其中ω0是同步角速度标幺制下取1P_m是机械功率P_e是电磁功率D是阻尼系数H是惯性时间常数。电磁功率对单机无穷大系统来说就是P_e (E * V_inf / X) * sin(δ)这里E是发电机暂态电动势V_inf是无穷大母线电压X是暂态电抗与线路电抗之和。整套模型的逻辑很清晰功角变化由转速偏差驱动转速变化由不平衡功率驱动而电磁功率又是功角的非线性函数——整个系统的非线性就集中在这条正弦链路上。量测方程测量什么呢考虑PMU可以直接量测节点相角和发电机输出功率我设置两组量测z1 P_e (E * V_inf / X) * sin(δ) z2 δ也就是说量测向量为[有功功率功角]。这样量测方程本身包含了非线性的P_e同时功角又是直接可测的既体现了非线性又保留了线性成分很适合用来对比EKF和UKF的性能。离散化我直接用了前向欧拉法。虽然欧拉法精度一般但在这个采样周期Ts0.01s的场景下离散化误差远小于量测噪声和滤波误差完全够用而且代码简单、状态转移雅可比矩阵写起来方便。如果你的场景对精度要求更高可以换用梯形法或四阶Runge-Kutta但对EKF来说状态转移函数越复杂雅可比矩阵推导就越痛苦这是需要考虑的。2.2 EKF的关键状态转移和量测的雅可比矩阵EKF的核心思想一句话就能说清在每一步的当前估计值处把非线性函数做一阶泰勒展开丢掉高阶项然后套用标准卡尔曼滤波公式。线性化需要两套雅可比矩阵状态转移雅可比F和量测雅可比H。对上述模型状态转移方程离散后为δ(k1) δ(k) Ts * ω0 * Δω(k) Δω(k1) Δω(k) Ts * (P_m - (E*V_inf/X)*sin(δ(k)) - D*Δω(k)) / (2H)注意这里的ω0如果取标幺值就是1如果取有名值则是2pi50。我在实现中统一用标幺值这样数值尺度比较接近协方差矩阵的调节也会轻松一些。状态转移雅可比F为2x2矩阵推导过程如下F11 ∂δ(k1)/∂δ(k) 1 F12 ∂δ(k1)/∂Δω(k) Ts * ω0 F21 ∂Δω(k1)/∂δ(k) -Ts * (E*V_inf/X) * cos(δ(k)) / (2H) F22 ∂Δω(k1)/∂Δω(k) 1 - Ts * D / (2H)量测雅可比H为2x2矩阵H11 ∂P_e/∂δ (E*V_inf/X) * cos(δ) H12 ∂P_e/∂Δω 0 H21 ∂δ/∂δ 1 H22 ∂δ/∂Δω 0每次滤波迭代中F和H都要用最新的状态估计值重新计算。这一点特别容易在写代码时忽略——如果图省事把F设成了常量矩阵状态估计在动态过程中基本就会发散。实测中H矩阵里的cos(δ)在功角变化时变化很剧烈尤其在系统受到大扰动后所以每个时刻刷新雅可比是EKF实现的红线。2.3 UKF的采样点思路绕开雅可比矩阵UKF不碰雅可比矩阵它走的是无迹变换Unscented Transform的路线既然直接变换一个概率分布很难那就从分布里抽几个代表性的点Sigma点把这些点一个个扔进非线性函数再从变换后的点集里重建均值和协方差。对于n维状态向量需要2n1个Sigma点。我的状态维数是2所以每步只需要5个点计算量非常小。Sigma点的生成公式如下λ α^2 * (n κ) - n χ0 x̄ χi x̄ (sqrt((nλ)*P_x))_i, i 1..n χin x̄ - (sqrt((nλ)*P_x))_i, i 1..n对应的权重为W0^m λ / (nλ) W0^c λ / (nλ) (1 - α^2 β) Wi^m Wi^c 1 / (2*(nλ)), i 1..2n参数选型上α控制Sigma点离均值的距离一般取1e-3到1之间的小值κ是比例参数通常取0或者3-nβ用来引入分布先验信息高斯分布下取2。我实际使用的是α0.1、κ0、β2这组参数在多数场景下表现稳定。UKF的时间更新和量测更新的写法本质上是传播Sigma点、加权重组两件事反复应用。全套流程不需要求导不需要推导任何H矩阵这是它在工程上最大的吸引力——代码通用性强换一个非线性模型只需要改f和h两个函数体。2.4 滤波器参数的初始化经验参数初始化这件事直接决定滤波器是收敛还是发散的起点。状态初值x0我建议从静态状态估计或者潮流结果里取不要拍脑袋随便给。我做仿真时先用潮流计算解出稳态功角加一个小的随机偏差作为初始估计值模拟从静态估计过渡到动态估计的真实过程。协方差初值P0反映的是对初始状态的不确定度。如果P0设得太小滤波器会对初始值过度自信后面量测来了也懒得修正如果太大前几步的增益很大估计值会被量测噪声带着剧烈抖动。我的经验是P0的量级跟初始偏差的平方对应。比如初始功角误差约0.05 rad那P0(1,1)取0.01左右比较合适转速误差约0.001P0(2,2)取1e-5附近。过程噪声Q和量测噪声R的整定是所有卡尔曼滤波实践里最需要耐心的工作。Q太小会导致滤波器只相信模型、不相信量测结果预测值一路飘R太小又会让估计值跟着每个噪声点抖动。我的调参套路是先用真值轨迹算一算每个状态量在动态过程中的波动幅度把Q设成波动幅度的1%左右量测R则根据传感器精度来功率量测误差按0.01标幺值考虑就设R(1,1)1e-4功角量测误差按0.005 rad考虑就设R(2,2)2.5e-5。之后再根据新息序列的实际标准差微调。3. 手把手拆解Matlab实现3.1 整体代码框架设计我习惯把所有东西按功能拆分不把逻辑堆在一个大脚本里。这个项目的Matlab代码结构如下main_DEKF_UKF.m % 主脚本设置参数生成真值跑滤波器画图 model_dyn.m % 状态方程连续/离散 model_meas.m % 量测方程 jacobian_F.m % EKF的状态转移雅可比 jacobian_H.m % EKF的量测雅可比 ekf_step.m % EKF单步滤波 ukf_step.m % UKF单步滤波 gen_measurements.m % 用真值加噪声生成量测序列 calc_rmse.m % 计算均方根误差主脚本里最核心的是设置好时间向量、仿真时长、采样周期、扰动设置比如1秒时机械功率阶跃、噪声统计量然后循环调用滤波函数。每个滤波器的输出都存成结构体最后统一画图比较。单步滤波函数建议写成这样的接口约定[x_est, P_est] ekf_step(x_est, P_est, z_k, Ts, params) [x_est, P_est] ukf_step(x_est, P_est, z_k, Ts, params)其中params是struct包含P_m、E_prime、V_inf、X、D、H_const等模型参数。这样写的好处是以后想换成IEEE节点系统或者多机系统只需要改模型函数滤波主循环完全不用动。3.2 EKF单步滤波核心代码EKF的完整单步实现如下这段代码我建议直接保存下来做模板function [x_est, P_est] ekf_step(x_est, P_est, z_k, Ts, params) % 提取模型参数 Pm params.Pm; E params.E; V params.V; Xt params.X; D params.D; H params.H; w0 params.w0; % 1. 状态预测时间更新 delta x_est(1); dw x_est(2); Pe (E*V/Xt) * sin(delta); x_pred zeros(2,1); x_pred(1) delta Ts * w0 * dw; x_pred(2) dw Ts * (Pm - Pe - D*dw) / (2*H); % 2. 状态转移雅可比用最新估计值计算 F [1, Ts*w0; -Ts*(E*V/Xt)*cos(delta)/(2*H), 1 - Ts*D/(2*H)]; % 3. 预测协方差 P_pred F * P_est * F params.Q; % 4. 量测预测与量测雅可比 delta_pred x_pred(1); h_pred [ (E*V/Xt)*sin(delta_pred); delta_pred ]; Hx [ (E*V/Xt)*cos(delta_pred), 0; 1, 0 ]; % 5. 卡尔曼增益与更新 S Hx * P_pred * Hx params.R; K P_pred * Hx / S; x_est x_pred K * (z_k - h_pred); P_est (eye(2) - K * Hx) * P_pred; end有几个容易踩的坑。第一个是第4步里算h_pred和第1步算Pe两处都用的是不同时刻的状态预测前用的是x_est预测后用x_pred如果你的量测方程里同时出现这两个状态务必分清。第二个是卡尔曼增益那行我用了/运算符而不是inv这样矩阵求逆在数值上更稳定也省去显式求逆的浮点开销。第三个是P_est更新我用了(I-K*H)*P_pred的简化形式这在增益为最优时没问题如果你的数值出现轻微不对称后续做Cholesky分解就会报错可以考虑用对称形式P_pred - K*S*K代价略高但数值更稳健。3.3 UKF单步滤波核心代码UKF的实现关键在Sigma点生成和权重计算。我的实现如下function [x_est, P_est] ukf_step(x_est, P_est, z_k, Ts, params) n length(x_est); alpha 0.1; beta 2; kappa 0; lambda alpha^2 * (n kappa) - n; % 1. 生成Sigma点 % 注意这里用sqrtm而不是chol避免协方差非正定时直接报错 A sqrtm((n lambda) * P_est); chi zeros(n, 2*n1); chi(:,1) x_est; for i 1:n chi(:,i1) x_est A(:,i); chi(:,in1) x_est - A(:,i); end % 权重 Wm zeros(1, 2*n1); Wc zeros(1, 2*n1); Wm(1) lambda / (n lambda); Wc(1) Wm(1) (1 - alpha^2 beta); for i 2:2*n1 Wm(i) 1 / (2*(n lambda)); Wc(i) Wm(i); end % 2. Sigma点通过状态方程传播时间更新 chi_pred zeros(n, 2*n1); for i 1:2*n1 chi_pred(:,i) model_dyn(chi(:,i), Ts, params); end x_pred chi_pred * Wm; P_pred (chi_pred - x_pred) * diag(Wc) * (chi_pred - x_pred) params.Q; % 3. Sigma点通过量测方程传播 z_pred_pts zeros(2, 2*n1); for i 1:2*n1 z_pred_pts(:,i) model_meas(chi_pred(:,i), params); end z_pred z_pred_pts * Wm; % 4. 计算量测协方差和状态-量测互协方差 Pzz (z_pred_pts - z_pred) * diag(Wc) * (z_pred_pts - z_pred) params.R; Pxz (chi_pred - x_pred) * diag(Wc) * (z_pred_pts - z_pred); % 5. 卡尔曼增益与更新 K Pxz / Pzz; x_est x_pred K * (z_k - z_pred); P_est P_pred - K * Pzz * K; endUKF的三个实现细节值得单独说说。第一Sigma点生成我用了sqrtm而不是chol。chol要求矩阵严格对称正定而滤波迭代中协方差矩阵往往因为数值截断出现轻微不对称或非正定chol会直接报错sqrtm在多数情况下可以容忍。代价是sqrtm计算量稍大但在n2时完全可以忽略。第二权重里β2是对高斯分布的修正项如果你处理的是非高斯量测噪声这个值需要重新标定。第三状态传播用的是model_dyn而不是f的离散化版本这里要保证model_dyn内部实现了完整的离散化包含Ts参数。model_meas函数按如下方式定义function z model_meas(x, params) delta x(1); z zeros(2,1); z(1) (params.E * params.V / params.X) * sin(delta); z(2) delta; end3.4 真值生成与量测模拟做仿真研究时真值轨迹必须独立于滤波器模型生成否则对比结果没有意义。我用了两种方法交叉验证第一种用Matlab的高精度求解器解连续模型。把状态方程写成ode45可用的函数句柄从0到3秒积分得到一条高精度轨迹作为真值。这样滤波器的离散模型与真值模型之间存在一定的模型失配恰好能检验滤波器对模型误差的鲁棒性。第二种把同一离散模型跑一遍但不加噪声当作理想离散真值。这种方法在滤波模型与真值模型完全一致的情况下验证滤波器是否能正确收敛排除模型失配的干扰。实际编码时我两种结果都跑了用理想离散真值检验算法实现的正确性用ode45连续真值检验算法的鲁棒性。量测序列就是真值加上独立高斯噪声z_noisy(1,k) z_true(1,k) sqrt(R(1,1)) * randn; z_noisy(2,k) z_true(2,k) sqrt(R(2,2)) * randn;性能评估用RMSErmse_delta sqrt(mean((delta_est - delta_true).^2)); rmse_dw sqrt(mean((dw_est - dw_true).^2));除了RMSE我还会记录每个滤波器的单步平均运行时间用tic/toc包住滤波循环除以总步数得到。这个指标在实时性要求高的场景下非常关键。4. 实测对比结果EKF和UKF的表现差异4.1 精度对比理想离散模型下的收敛行为先看理想离散模型下的结果。初值故意偏离真值比如功角初始误差0.05 rad转速偏差初始误差0.002 p.u.。两条滤波器的收敛曲线都能在0.5秒内收敛到真值附近但收敛过程有明显差异。EKF在收敛初期出现了一个小幅度的超调原因是初始点离真值较远时一阶线性化对正弦函数的近似误差较大增益计算偏保守修正力度不足UKF则因为Sigma点直接落在状态分布中对非线性的拟合更准确收敛更快且几乎没有超调。稳态阶段两者的RMSE都维持在量测噪声允许的范围内差异不再显著。具体数值用表格列出单次仿真Ts0.01s噪声统计同上指标EKFUKF功角RMSE (rad)0.00410.0032转速RMSE (p.u.)2.1e-41.7e-4收敛时间 (s)0.450.30平均单步耗时 (ms)0.080.15可以看到UKF在精度上略胜一筹代价是单步耗时约为EKF的两倍。在n2的模型里这个耗时差异基本可以忽略。4.2 大扰动场景下的鲁棒性对比真正拉开差距的是大扰动场景。我在t1s时让机械功率从0.8阶跃到1.2 p.u.相当于给系统施加一个冲击功角从稳态值开始大幅摆动摆幅超过0.3 rad。这种情况下状态在很宽的范围内快速扫过正弦非线性被充分激发。EKF在这个场景里的表现不太理想。在功角摆过峰值、变化率最大的区间一阶线性化误差急剧放大滤波器出现了明显的跟踪滞后估计值在暂态峰值处比真值低了约0.02 rad而且滞后现象持续了大半个振荡周期。如果把扰动再加猛一些功角越过稳定边界EKF甚至可能发散。这并非EKF的代码问题而是线性化近似的固有局限——你可以把EKF理解成在弯道上走切线弯道缓时没问题弯道急时误差自然增大。UKF在同一场景下则稳得多。Sigma点覆盖了状态分布的非线性区域传播后的均值和协方差能反映真实分布的大致形状因此在大摆动区间依然保持了良好的跟踪精度峰值处的估计误差仅为EKF的三分之一左右。这个对比结论在文献里反复出现但自己亲手跑出来才有直观感受。4.3 计算效率与工程适用性讨论从表里可以看到UKF单步耗时约为EKF的两倍但这只是在2维状态下的结果。如果扩展到多机系统状态维数n上升到几十维UKF需要的Sigma点数量是2n1而EKF还需要推导和计算2n维雅可比矩阵两者计算量的差距会缩小。在状态维数超过20的系统中UKF的计算优势会逐渐显现因为不需要求解复杂的雅可比矩阵解析式也不存在雅可比矩阵求导错误的风险。工程上选择哪种算法我的判断标准有三条一看非线性强度扰动大、模型非线性强优先UKF二看实现约束如果没人手推导雅可比矩阵、又急需一个能快速上手的通用滤波器优先UKF三看计算预算极端实时的嵌入式环境对计算量极敏感且工况相对平稳EKF依然是值得保留的选项。5. 常见报错与调参经验我踩过的坑5.1 滤波器发散几乎所有问题最后都归结到Q和R最大的坑就是发散。表现是估计值突然偏离真值而且越跑越远协方差矩阵却还在不断缩小表示滤波器自信地错着。出现这种情况第一反应不是去调初值而是检查Q和R的比例。发散的本质是滤波器对模型的信任程度和对量测的信任程度失衡。如果Q设得特别小滤波器认为模型十分准确量测增益被压低量测修正的作用微乎其微模型一旦有偏这个偏差永远得不到修正最终发散。反之如果R设得特别小滤波器对量测过度信任每次都被噪声牵着走状态估计的噪声很大但一般不会发散只是曲线看起来像锯齿。我摸索出一套实用的调参流程先把R按传感器精度标准设好不要动。把Q从小到大扫一遍跑同一条真值轨迹画出滤波误差的曲线看误差是否收敛到量测噪声水平。如果误差曲线在收敛后仍然有鼓包说明Q偏小模型跟不上动态如果误差曲线噪声特别尖锐说明Q偏大被过程噪声污染了。最后用新息序列innovation即z_k - h_pred做一次白噪检查均值应接近零自相关函数应接近冲激。若新息均值明显非零基本可以断定模型偏差或者Q过小。这套流程听起来朴素但比盲目调参高效得多。5.2 雅可比矩阵推导和验证的教训EKF的实现中雅可比矩阵是事故高发区。解析推导一旦出错滤波器不会直接报错而是在某些区间表现出奇怪的收敛行为非常难排查。我的做法是用符号工具箱先验算一遍。比如对F21那个项手写时容易漏掉分母2H或者把cos(δ)误写成sin(δ)。用符号求导得到表达式后再生成数值函数这样准确率大大提升。符号推导之后还要用数值差分做交叉验证epsilon 1e-6; F_num zeros(2,2); for j 1:2 x_plus x0; x_minus x0; x_plus(j) x_plus(j) epsilon; x_minus(j) x_minus(j) - epsilon; F_num(:,j) (model_dyn(x_plus, Ts, params) - model_dyn(x_minus, Ts, params)) / (2*epsilon); end把数值F_num和解析F放一起比较差异在1e-6量级才算过关。这一步花不了多少时间但能省掉大量调试时间。5.3 协方差矩阵非正定UKF的Cholesky陷阱UKF在迭代几十步后有时会突然报错Matrix must be positive definite这通常是协方差矩阵在数值上失去了正定性。原因有两个一是权重Wc中存在负值导致传播后的协方差计算出现细微负特征值二是重复的矩阵乘法累积了浮点误差使矩阵轻微不对称。解决手段按优先顺序排列把chol换成sqrtm是第一个能快速见效的改动在协方差更新后强制对称化P_est (P_est P_est) / 2;在特征值出现非正时加一个微小的单位阵抖动P_est P_est 1e-12 * eye(n);其中第2条我每次都加几乎零成本能明显降低后期发散的概率。第3条属于最后的手段加入的抖动太大会污染估计精度不能滥用。5.4 参数速查表与推荐配置把整套仿真最常用的配置汇总成一张速查表方便复现参数符号取值说明采样周期Ts0.01 s髓常介于0.005~0.02同步角速度ω01 p.u.有名值为100π惯性时间常数H5.0 s典型火电机组阻尼系数D2.0 p.u.含阻尼绕组效果暂态电动势E1.05 p.u.可微调无穷大母线电压V_inf1.0 p.u.基准值总电抗X0.6 p.u.发电机线路机械功率P_m0.8 p.u.t1s跳变到1.2过程噪声Qdiag([1e-5, 1e-4])先按此调量测噪声Rdiag([1e-4, 2.5e-5])对应0.01功率误差UKF参数α/β/κ0.1/2/0高斯噪声推荐这组参数在绝大多数CPU上跑3秒仿真300步只需要眨眼功夫很适合反复试验。6. 最后一点经验与扩展方向实际把EKF和UKF都调通之后我对这两类滤波器的理解比读十篇综述都深。EKF的代码短小精悍在工况平稳时效率极高但它的一切弱点都来源于线性化这个假设UKF看起来只多了几十行代码却在非线性场景里提供了显著的稳定性和精度提升。这个对比心得比任何结论性的话都更有价值。如果后续想扩展这个项目我建议从三个方向入手。第一把单机模型换成IEEE 14节点或39节点系统状态维数增加后可以真正考察两种算法在多机协调估计中的表现第二把量测扩展为含PMU和SCADA的混合量测处理不同刷新率下的异步数据融合第三把UKF的无迹变换思想扩展到强跟踪滤波STF或自适应卡尔曼滤波针对模型参数不确定的场景做在线调整。Matlab这套框架迁移起来非常顺手——模型的接口已经拆开滤波器的核心函数不需要大改。我自己的下一步计划是加入不良数据检测模块毕竟电力系统实际运行中PMU偶尔丢包或者出坏点这才是工程落地真正要面对的问题。 Matlab实现EKF与UKF的电力系统动态状态估计对比分析
1. 立项思路动态状态估计在电力系统里到底解决什么问题1.1 动态状态估计 vs 静态状态估计差在哪儿先说说背景。做过电力系统状态估计的人都知道传统做法是加权最小二乘WLS那一套输入SCADA传来的量测值输出一个断面的电压幅值和相角。这个方案在稳态监控里没问题但它的本质是用一个时间断面的冗余量测去拟合一个状态解对时间的演化完全没有概念。系统一旦进入动态过程比如负荷突变、发电机甩负荷、线路故障后的暂态SCADA的刷新率本来就低秒级甚至分钟级静态估计的结果往往滞后一拍等算出来系统可能已经跑到另一个工况去了。动态状态估计的思路完全不同。它把系统自身的运动规律写进模型里——最常见的就是发电机转子运动方程——然后用卡尔曼滤波这类递推算法把模型预测和新到的量测融合起来。每一步估计都包含两个动作先用状态方程把上一时刻的状态推进到当前时刻再用当前时刻的量测去修正这个预测值。这样做的结果就是估计值有预测能力能跟踪动态过程而且天然带有平滑噪声的功能。我做这个课题的时候目标很明确在Matlab里实现两种非线性卡尔曼滤波器用于电力系统动态状态估计并且把它们的精度、稳定性和实现成本做一次横向对比。之所以同时做EKF和UKF就是因为这两者是处理非线性状态估计时最经典的两条路线——一个靠线性化打补丁一个靠采样点传播理解清楚这两条路线的差异后续做扩展卡尔曼滤波的改进型如无迹变换、粒子滤波都有了参照。1.2 算法选型为什么不直接上粒子滤波现在做状态估计可选的非线性滤波算法其实很多粒子滤波、容积卡尔曼CKF、甚至神经网络辅助的估计器都有。但在这个项目里选EKF和UKF是有实际考虑的。首先是模型特点。电力系统动态估计的核心非线性体现在电磁功率与功角之间的正弦关系这种非线性程度属于中等偏强。EKF的适用前提是系统在估计点附近近似可线性化只要扰动不是特别剧烈功角不会大范围跳变一阶泰勒展开的精度基本够用。而UKF不做线性化直接用Sigma点去传播状态分布对正弦非线性天然更友好理论上可以达到二阶精度。其次是实现成本。粒子滤波虽然能处理强非线性但需要大量粒子在Matlab里做几百个粒子的递推估计单步计算的耗时就会明显上升而我的仿真场景要求Ts取到10毫秒甚至更小计算效率必须考虑。EKF和UKF的计算量都可控UKF对n维状态只需要2n1个Sigma点在状态维数不高的电力系统模型中非常划算。最后是想做一个标准答案级别的对比。EKF的线性化误差在强非线性场景下会暴露UKF则在Sigma点参数选择不当时也会出现协方差非正定的问题。两套算法放在同一个仿真环境里用相同的真值轨迹、相同的噪声统计对比结果才有说服力。这个对比数据无论写论文还是做汇报都是很有分量的支撑。2. 数学模型和算法原理把滤波器和电力系统模型接上2.1 动态估计用的状态方程和量测方程先说状态方程。我采用经典的二阶发电机模型即转子摇摆方程swing equation。状态量取两维x1 δ发电机功角radx2 Δω角速度偏差p.u.基准为同步角速度连续时间模型为dδ/dt ω0 * Δω dΔω/dt (P_m - P_e - D * Δω) / (2H)其中ω0是同步角速度标幺制下取1P_m是机械功率P_e是电磁功率D是阻尼系数H是惯性时间常数。电磁功率对单机无穷大系统来说就是P_e (E * V_inf / X) * sin(δ)这里E是发电机暂态电动势V_inf是无穷大母线电压X是暂态电抗与线路电抗之和。整套模型的逻辑很清晰功角变化由转速偏差驱动转速变化由不平衡功率驱动而电磁功率又是功角的非线性函数——整个系统的非线性就集中在这条正弦链路上。量测方程测量什么呢考虑PMU可以直接量测节点相角和发电机输出功率我设置两组量测z1 P_e (E * V_inf / X) * sin(δ) z2 δ也就是说量测向量为[有功功率功角]。这样量测方程本身包含了非线性的P_e同时功角又是直接可测的既体现了非线性又保留了线性成分很适合用来对比EKF和UKF的性能。离散化我直接用了前向欧拉法。虽然欧拉法精度一般但在这个采样周期Ts0.01s的场景下离散化误差远小于量测噪声和滤波误差完全够用而且代码简单、状态转移雅可比矩阵写起来方便。如果你的场景对精度要求更高可以换用梯形法或四阶Runge-Kutta但对EKF来说状态转移函数越复杂雅可比矩阵推导就越痛苦这是需要考虑的。2.2 EKF的关键状态转移和量测的雅可比矩阵EKF的核心思想一句话就能说清在每一步的当前估计值处把非线性函数做一阶泰勒展开丢掉高阶项然后套用标准卡尔曼滤波公式。线性化需要两套雅可比矩阵状态转移雅可比F和量测雅可比H。对上述模型状态转移方程离散后为δ(k1) δ(k) Ts * ω0 * Δω(k) Δω(k1) Δω(k) Ts * (P_m - (E*V_inf/X)*sin(δ(k)) - D*Δω(k)) / (2H)注意这里的ω0如果取标幺值就是1如果取有名值则是2pi50。我在实现中统一用标幺值这样数值尺度比较接近协方差矩阵的调节也会轻松一些。状态转移雅可比F为2x2矩阵推导过程如下F11 ∂δ(k1)/∂δ(k) 1 F12 ∂δ(k1)/∂Δω(k) Ts * ω0 F21 ∂Δω(k1)/∂δ(k) -Ts * (E*V_inf/X) * cos(δ(k)) / (2H) F22 ∂Δω(k1)/∂Δω(k) 1 - Ts * D / (2H)量测雅可比H为2x2矩阵H11 ∂P_e/∂δ (E*V_inf/X) * cos(δ) H12 ∂P_e/∂Δω 0 H21 ∂δ/∂δ 1 H22 ∂δ/∂Δω 0每次滤波迭代中F和H都要用最新的状态估计值重新计算。这一点特别容易在写代码时忽略——如果图省事把F设成了常量矩阵状态估计在动态过程中基本就会发散。实测中H矩阵里的cos(δ)在功角变化时变化很剧烈尤其在系统受到大扰动后所以每个时刻刷新雅可比是EKF实现的红线。2.3 UKF的采样点思路绕开雅可比矩阵UKF不碰雅可比矩阵它走的是无迹变换Unscented Transform的路线既然直接变换一个概率分布很难那就从分布里抽几个代表性的点Sigma点把这些点一个个扔进非线性函数再从变换后的点集里重建均值和协方差。对于n维状态向量需要2n1个Sigma点。我的状态维数是2所以每步只需要5个点计算量非常小。Sigma点的生成公式如下λ α^2 * (n κ) - n χ0 x̄ χi x̄ (sqrt((nλ)*P_x))_i, i 1..n χin x̄ - (sqrt((nλ)*P_x))_i, i 1..n对应的权重为W0^m λ / (nλ) W0^c λ / (nλ) (1 - α^2 β) Wi^m Wi^c 1 / (2*(nλ)), i 1..2n参数选型上α控制Sigma点离均值的距离一般取1e-3到1之间的小值κ是比例参数通常取0或者3-nβ用来引入分布先验信息高斯分布下取2。我实际使用的是α0.1、κ0、β2这组参数在多数场景下表现稳定。UKF的时间更新和量测更新的写法本质上是传播Sigma点、加权重组两件事反复应用。全套流程不需要求导不需要推导任何H矩阵这是它在工程上最大的吸引力——代码通用性强换一个非线性模型只需要改f和h两个函数体。2.4 滤波器参数的初始化经验参数初始化这件事直接决定滤波器是收敛还是发散的起点。状态初值x0我建议从静态状态估计或者潮流结果里取不要拍脑袋随便给。我做仿真时先用潮流计算解出稳态功角加一个小的随机偏差作为初始估计值模拟从静态估计过渡到动态估计的真实过程。协方差初值P0反映的是对初始状态的不确定度。如果P0设得太小滤波器会对初始值过度自信后面量测来了也懒得修正如果太大前几步的增益很大估计值会被量测噪声带着剧烈抖动。我的经验是P0的量级跟初始偏差的平方对应。比如初始功角误差约0.05 rad那P0(1,1)取0.01左右比较合适转速误差约0.001P0(2,2)取1e-5附近。过程噪声Q和量测噪声R的整定是所有卡尔曼滤波实践里最需要耐心的工作。Q太小会导致滤波器只相信模型、不相信量测结果预测值一路飘R太小又会让估计值跟着每个噪声点抖动。我的调参套路是先用真值轨迹算一算每个状态量在动态过程中的波动幅度把Q设成波动幅度的1%左右量测R则根据传感器精度来功率量测误差按0.01标幺值考虑就设R(1,1)1e-4功角量测误差按0.005 rad考虑就设R(2,2)2.5e-5。之后再根据新息序列的实际标准差微调。3. 手把手拆解Matlab实现3.1 整体代码框架设计我习惯把所有东西按功能拆分不把逻辑堆在一个大脚本里。这个项目的Matlab代码结构如下main_DEKF_UKF.m % 主脚本设置参数生成真值跑滤波器画图 model_dyn.m % 状态方程连续/离散 model_meas.m % 量测方程 jacobian_F.m % EKF的状态转移雅可比 jacobian_H.m % EKF的量测雅可比 ekf_step.m % EKF单步滤波 ukf_step.m % UKF单步滤波 gen_measurements.m % 用真值加噪声生成量测序列 calc_rmse.m % 计算均方根误差主脚本里最核心的是设置好时间向量、仿真时长、采样周期、扰动设置比如1秒时机械功率阶跃、噪声统计量然后循环调用滤波函数。每个滤波器的输出都存成结构体最后统一画图比较。单步滤波函数建议写成这样的接口约定[x_est, P_est] ekf_step(x_est, P_est, z_k, Ts, params) [x_est, P_est] ukf_step(x_est, P_est, z_k, Ts, params)其中params是struct包含P_m、E_prime、V_inf、X、D、H_const等模型参数。这样写的好处是以后想换成IEEE节点系统或者多机系统只需要改模型函数滤波主循环完全不用动。3.2 EKF单步滤波核心代码EKF的完整单步实现如下这段代码我建议直接保存下来做模板function [x_est, P_est] ekf_step(x_est, P_est, z_k, Ts, params) % 提取模型参数 Pm params.Pm; E params.E; V params.V; Xt params.X; D params.D; H params.H; w0 params.w0; % 1. 状态预测时间更新 delta x_est(1); dw x_est(2); Pe (E*V/Xt) * sin(delta); x_pred zeros(2,1); x_pred(1) delta Ts * w0 * dw; x_pred(2) dw Ts * (Pm - Pe - D*dw) / (2*H); % 2. 状态转移雅可比用最新估计值计算 F [1, Ts*w0; -Ts*(E*V/Xt)*cos(delta)/(2*H), 1 - Ts*D/(2*H)]; % 3. 预测协方差 P_pred F * P_est * F params.Q; % 4. 量测预测与量测雅可比 delta_pred x_pred(1); h_pred [ (E*V/Xt)*sin(delta_pred); delta_pred ]; Hx [ (E*V/Xt)*cos(delta_pred), 0; 1, 0 ]; % 5. 卡尔曼增益与更新 S Hx * P_pred * Hx params.R; K P_pred * Hx / S; x_est x_pred K * (z_k - h_pred); P_est (eye(2) - K * Hx) * P_pred; end有几个容易踩的坑。第一个是第4步里算h_pred和第1步算Pe两处都用的是不同时刻的状态预测前用的是x_est预测后用x_pred如果你的量测方程里同时出现这两个状态务必分清。第二个是卡尔曼增益那行我用了/运算符而不是inv这样矩阵求逆在数值上更稳定也省去显式求逆的浮点开销。第三个是P_est更新我用了(I-K*H)*P_pred的简化形式这在增益为最优时没问题如果你的数值出现轻微不对称后续做Cholesky分解就会报错可以考虑用对称形式P_pred - K*S*K代价略高但数值更稳健。3.3 UKF单步滤波核心代码UKF的实现关键在Sigma点生成和权重计算。我的实现如下function [x_est, P_est] ukf_step(x_est, P_est, z_k, Ts, params) n length(x_est); alpha 0.1; beta 2; kappa 0; lambda alpha^2 * (n kappa) - n; % 1. 生成Sigma点 % 注意这里用sqrtm而不是chol避免协方差非正定时直接报错 A sqrtm((n lambda) * P_est); chi zeros(n, 2*n1); chi(:,1) x_est; for i 1:n chi(:,i1) x_est A(:,i); chi(:,in1) x_est - A(:,i); end % 权重 Wm zeros(1, 2*n1); Wc zeros(1, 2*n1); Wm(1) lambda / (n lambda); Wc(1) Wm(1) (1 - alpha^2 beta); for i 2:2*n1 Wm(i) 1 / (2*(n lambda)); Wc(i) Wm(i); end % 2. Sigma点通过状态方程传播时间更新 chi_pred zeros(n, 2*n1); for i 1:2*n1 chi_pred(:,i) model_dyn(chi(:,i), Ts, params); end x_pred chi_pred * Wm; P_pred (chi_pred - x_pred) * diag(Wc) * (chi_pred - x_pred) params.Q; % 3. Sigma点通过量测方程传播 z_pred_pts zeros(2, 2*n1); for i 1:2*n1 z_pred_pts(:,i) model_meas(chi_pred(:,i), params); end z_pred z_pred_pts * Wm; % 4. 计算量测协方差和状态-量测互协方差 Pzz (z_pred_pts - z_pred) * diag(Wc) * (z_pred_pts - z_pred) params.R; Pxz (chi_pred - x_pred) * diag(Wc) * (z_pred_pts - z_pred); % 5. 卡尔曼增益与更新 K Pxz / Pzz; x_est x_pred K * (z_k - z_pred); P_est P_pred - K * Pzz * K; endUKF的三个实现细节值得单独说说。第一Sigma点生成我用了sqrtm而不是chol。chol要求矩阵严格对称正定而滤波迭代中协方差矩阵往往因为数值截断出现轻微不对称或非正定chol会直接报错sqrtm在多数情况下可以容忍。代价是sqrtm计算量稍大但在n2时完全可以忽略。第二权重里β2是对高斯分布的修正项如果你处理的是非高斯量测噪声这个值需要重新标定。第三状态传播用的是model_dyn而不是f的离散化版本这里要保证model_dyn内部实现了完整的离散化包含Ts参数。model_meas函数按如下方式定义function z model_meas(x, params) delta x(1); z zeros(2,1); z(1) (params.E * params.V / params.X) * sin(delta); z(2) delta; end3.4 真值生成与量测模拟做仿真研究时真值轨迹必须独立于滤波器模型生成否则对比结果没有意义。我用了两种方法交叉验证第一种用Matlab的高精度求解器解连续模型。把状态方程写成ode45可用的函数句柄从0到3秒积分得到一条高精度轨迹作为真值。这样滤波器的离散模型与真值模型之间存在一定的模型失配恰好能检验滤波器对模型误差的鲁棒性。第二种把同一离散模型跑一遍但不加噪声当作理想离散真值。这种方法在滤波模型与真值模型完全一致的情况下验证滤波器是否能正确收敛排除模型失配的干扰。实际编码时我两种结果都跑了用理想离散真值检验算法实现的正确性用ode45连续真值检验算法的鲁棒性。量测序列就是真值加上独立高斯噪声z_noisy(1,k) z_true(1,k) sqrt(R(1,1)) * randn; z_noisy(2,k) z_true(2,k) sqrt(R(2,2)) * randn;性能评估用RMSErmse_delta sqrt(mean((delta_est - delta_true).^2)); rmse_dw sqrt(mean((dw_est - dw_true).^2));除了RMSE我还会记录每个滤波器的单步平均运行时间用tic/toc包住滤波循环除以总步数得到。这个指标在实时性要求高的场景下非常关键。4. 实测对比结果EKF和UKF的表现差异4.1 精度对比理想离散模型下的收敛行为先看理想离散模型下的结果。初值故意偏离真值比如功角初始误差0.05 rad转速偏差初始误差0.002 p.u.。两条滤波器的收敛曲线都能在0.5秒内收敛到真值附近但收敛过程有明显差异。EKF在收敛初期出现了一个小幅度的超调原因是初始点离真值较远时一阶线性化对正弦函数的近似误差较大增益计算偏保守修正力度不足UKF则因为Sigma点直接落在状态分布中对非线性的拟合更准确收敛更快且几乎没有超调。稳态阶段两者的RMSE都维持在量测噪声允许的范围内差异不再显著。具体数值用表格列出单次仿真Ts0.01s噪声统计同上指标EKFUKF功角RMSE (rad)0.00410.0032转速RMSE (p.u.)2.1e-41.7e-4收敛时间 (s)0.450.30平均单步耗时 (ms)0.080.15可以看到UKF在精度上略胜一筹代价是单步耗时约为EKF的两倍。在n2的模型里这个耗时差异基本可以忽略。4.2 大扰动场景下的鲁棒性对比真正拉开差距的是大扰动场景。我在t1s时让机械功率从0.8阶跃到1.2 p.u.相当于给系统施加一个冲击功角从稳态值开始大幅摆动摆幅超过0.3 rad。这种情况下状态在很宽的范围内快速扫过正弦非线性被充分激发。EKF在这个场景里的表现不太理想。在功角摆过峰值、变化率最大的区间一阶线性化误差急剧放大滤波器出现了明显的跟踪滞后估计值在暂态峰值处比真值低了约0.02 rad而且滞后现象持续了大半个振荡周期。如果把扰动再加猛一些功角越过稳定边界EKF甚至可能发散。这并非EKF的代码问题而是线性化近似的固有局限——你可以把EKF理解成在弯道上走切线弯道缓时没问题弯道急时误差自然增大。UKF在同一场景下则稳得多。Sigma点覆盖了状态分布的非线性区域传播后的均值和协方差能反映真实分布的大致形状因此在大摆动区间依然保持了良好的跟踪精度峰值处的估计误差仅为EKF的三分之一左右。这个对比结论在文献里反复出现但自己亲手跑出来才有直观感受。4.3 计算效率与工程适用性讨论从表里可以看到UKF单步耗时约为EKF的两倍但这只是在2维状态下的结果。如果扩展到多机系统状态维数n上升到几十维UKF需要的Sigma点数量是2n1而EKF还需要推导和计算2n维雅可比矩阵两者计算量的差距会缩小。在状态维数超过20的系统中UKF的计算优势会逐渐显现因为不需要求解复杂的雅可比矩阵解析式也不存在雅可比矩阵求导错误的风险。工程上选择哪种算法我的判断标准有三条一看非线性强度扰动大、模型非线性强优先UKF二看实现约束如果没人手推导雅可比矩阵、又急需一个能快速上手的通用滤波器优先UKF三看计算预算极端实时的嵌入式环境对计算量极敏感且工况相对平稳EKF依然是值得保留的选项。5. 常见报错与调参经验我踩过的坑5.1 滤波器发散几乎所有问题最后都归结到Q和R最大的坑就是发散。表现是估计值突然偏离真值而且越跑越远协方差矩阵却还在不断缩小表示滤波器自信地错着。出现这种情况第一反应不是去调初值而是检查Q和R的比例。发散的本质是滤波器对模型的信任程度和对量测的信任程度失衡。如果Q设得特别小滤波器认为模型十分准确量测增益被压低量测修正的作用微乎其微模型一旦有偏这个偏差永远得不到修正最终发散。反之如果R设得特别小滤波器对量测过度信任每次都被噪声牵着走状态估计的噪声很大但一般不会发散只是曲线看起来像锯齿。我摸索出一套实用的调参流程先把R按传感器精度标准设好不要动。把Q从小到大扫一遍跑同一条真值轨迹画出滤波误差的曲线看误差是否收敛到量测噪声水平。如果误差曲线在收敛后仍然有鼓包说明Q偏小模型跟不上动态如果误差曲线噪声特别尖锐说明Q偏大被过程噪声污染了。最后用新息序列innovation即z_k - h_pred做一次白噪检查均值应接近零自相关函数应接近冲激。若新息均值明显非零基本可以断定模型偏差或者Q过小。这套流程听起来朴素但比盲目调参高效得多。5.2 雅可比矩阵推导和验证的教训EKF的实现中雅可比矩阵是事故高发区。解析推导一旦出错滤波器不会直接报错而是在某些区间表现出奇怪的收敛行为非常难排查。我的做法是用符号工具箱先验算一遍。比如对F21那个项手写时容易漏掉分母2H或者把cos(δ)误写成sin(δ)。用符号求导得到表达式后再生成数值函数这样准确率大大提升。符号推导之后还要用数值差分做交叉验证epsilon 1e-6; F_num zeros(2,2); for j 1:2 x_plus x0; x_minus x0; x_plus(j) x_plus(j) epsilon; x_minus(j) x_minus(j) - epsilon; F_num(:,j) (model_dyn(x_plus, Ts, params) - model_dyn(x_minus, Ts, params)) / (2*epsilon); end把数值F_num和解析F放一起比较差异在1e-6量级才算过关。这一步花不了多少时间但能省掉大量调试时间。5.3 协方差矩阵非正定UKF的Cholesky陷阱UKF在迭代几十步后有时会突然报错Matrix must be positive definite这通常是协方差矩阵在数值上失去了正定性。原因有两个一是权重Wc中存在负值导致传播后的协方差计算出现细微负特征值二是重复的矩阵乘法累积了浮点误差使矩阵轻微不对称。解决手段按优先顺序排列把chol换成sqrtm是第一个能快速见效的改动在协方差更新后强制对称化P_est (P_est P_est) / 2;在特征值出现非正时加一个微小的单位阵抖动P_est P_est 1e-12 * eye(n);其中第2条我每次都加几乎零成本能明显降低后期发散的概率。第3条属于最后的手段加入的抖动太大会污染估计精度不能滥用。5.4 参数速查表与推荐配置把整套仿真最常用的配置汇总成一张速查表方便复现参数符号取值说明采样周期Ts0.01 s髓常介于0.005~0.02同步角速度ω01 p.u.有名值为100π惯性时间常数H5.0 s典型火电机组阻尼系数D2.0 p.u.含阻尼绕组效果暂态电动势E1.05 p.u.可微调无穷大母线电压V_inf1.0 p.u.基准值总电抗X0.6 p.u.发电机线路机械功率P_m0.8 p.u.t1s跳变到1.2过程噪声Qdiag([1e-5, 1e-4])先按此调量测噪声Rdiag([1e-4, 2.5e-5])对应0.01功率误差UKF参数α/β/κ0.1/2/0高斯噪声推荐这组参数在绝大多数CPU上跑3秒仿真300步只需要眨眼功夫很适合反复试验。6. 最后一点经验与扩展方向实际把EKF和UKF都调通之后我对这两类滤波器的理解比读十篇综述都深。EKF的代码短小精悍在工况平稳时效率极高但它的一切弱点都来源于线性化这个假设UKF看起来只多了几十行代码却在非线性场景里提供了显著的稳定性和精度提升。这个对比心得比任何结论性的话都更有价值。如果后续想扩展这个项目我建议从三个方向入手。第一把单机模型换成IEEE 14节点或39节点系统状态维数增加后可以真正考察两种算法在多机协调估计中的表现第二把量测扩展为含PMU和SCADA的混合量测处理不同刷新率下的异步数据融合第三把UKF的无迹变换思想扩展到强跟踪滤波STF或自适应卡尔曼滤波针对模型参数不确定的场景做在线调整。Matlab这套框架迁移起来非常顺手——模型的接口已经拆开滤波器的核心函数不需要大改。我自己的下一步计划是加入不良数据检测模块毕竟电力系统实际运行中PMU偶尔丢包或者出坏点这才是工程落地真正要面对的问题。