做电力系统动态状态估计DSE这件事最早吸引我入坑的原因是调度中心拿到的PMU/SCADA量测永远带着噪声你没法直接把数值拿去做控制或保护判断。传统静态估计只给一个时间断面的快照但电力系统是连续变化的尤其故障、负荷突变那几秒状态变化的剧烈程度跟稳态完全不在一个量级。这时候就需要递推滤波也就是动态状态估计。讲到非线性系统的递推滤波EKF和UKF是最绕不开的两条路也是这个项目标题的核心。这篇博文围绕基于EKF和UKF的电力系统动态状态估计完整梳理我从状态建模、滤波原理、Matlab代码实现到调参避坑的整个过程适合刚接触DSE、或者已经跑通静态状态估计想进阶到动态估计的同学参考目标是让你看完能自己把代码拼起来跑出结果。我会把几个关键决策讲透为什么用发电机二阶摆动方程做状态模型、Jacobian矩阵为什么要分解析和数值两种写法、UKF的sigma点参数怎么取才能保证协方差正定、以及实测里哪些坑会让你一晚上白干。文中所有Matlab代码片段都是围绕核心逻辑写的方向明确照抄拼装即可运行但更重要的是看懂每一步到底在算什么。1. 为什么电力系统需要动态状态估计从断面快照到实时追踪1.1 静态状态估计的局限性传统加权最小二乘WLS状态估计做的事情是在一个时间断面内利用冗余量测求解一个使加权残差平方和最小的状态向量。它只用一个时刻的数据没有任何时序信息。这个方案在稳态工况下完全够用但有三类问题很要命对坏数据敏感一个通道的尖峰毛刺可能把整个断面的解拉偏事后要靠残差检测去揪。不能预测只能回答现在最可能是多少回答不了下一秒大概是多少。没有动态模型扰动期间的状态轨迹只能靠事后离线分析没法在线跟随。而动态状态估计的思路是把问题放进递推框架里用上一时刻的状态预测当前时刻再用当前时刻的量测修正预测。这就是卡尔曼滤波家族的看家本领。1.2 动态状态估计的数学框架标准离散时间模型长这样状态转移方程x_{k1} f(x_k) w_k其中 w_k ~ N(0, Q)量测方程z_k h(x_k) v_k其中 v_k ~ N(0, R)在电力系统DSE里最常见的状态量选取是发电机的转子运动状态。经典二阶模型摆动方程连续时间形式dδ/dt ω - ω0dω/dt (Pm - Pe - D(ω - ω0)) / M其中 δ 是发电机功角相对参考机或惯量中心ω 是角频率Pm 是机械功率Pe 是电磁功率M 是惯性时间常数D 是阻尼系数。这里值得多说一句要不要把励磁系统、调速器动态也建模进去我的建议是初学阶段先不要。把Pm当作已知常量或慢变输入把模型没刻画的那部分差异交给过程噪声Q去吸收。先把EKF/UKF跑通后面再逐步往状态方程里加励磁动态否则你区分不清滤波发散是因为算法写错还是模型没建对。离散化用一阶欧拉法最简单δ_{k1} δ_k Δt(ω_k - ω0)ω_{k1} ω_k Δt((Pm - Pe_k - D(ω_k - ω0)) / M)Pe 对 δ 的函数关系在网络里是非线性的。用单机无穷大母线SMIB模型举例Pe C * sin(δ)C E * V / x_total比如取 E1.0 p.u.V1.0 p.u.x_total0.5 p.u.则 C2.0Pe 2.0 sin(δ)。这就是一个非常清晰的非线性来源。多机系统里Pe要通过网络方程和潮流联立求解量测函数更复杂但滤波逻辑完全一致。1.3 量测配置直接量测与非线性量测PMU能给的量测包括功角、角频率、电压幅值、有功功率、无功功率等。在DSE里可以把量测分成两类直接量测δ、ω对应 h(x) 是线性的非线性量测有功功率 P、无功功率 Q、电压幅值 V对应 h(x) 是状态的非线性函数做算法对比实验时我的建议是量测集合里至少保留一个非线性量测项比如发电机有功功率Pe。如果全部都是直接线性量测EKF和UKF的估计结果几乎完全一致你根本看不出两种算法的差异后面第四章会详细演示这一点。2. EKF的计算骨架一阶线性化与Jacobian矩阵的实现细节2.1 线性化的核心思路EKF的思想很直接非线性函数在当前估计点附近做一阶泰勒展开丢掉高阶项得到近似的线性系统然后套用标准卡尔曼滤波的预测-修正两步。预测步x_pred f(x_est)P_pred A * P_est * A Q其中 A 是状态转移函数 f 在 x_est 处的Jacobian矩阵A ∂f/∂x。修正步H ∂h/∂x在 x_pred 处取值S H * P_pred * H RK P_pred * H * S^{-1}x_est x_pred K * (z_meas - h(x_pred))P_est (I - K*H) * P_pred流程看起来简单真正的工作量全在Jacobian矩阵的推导和实现上。2.2 Jacobian矩阵解析法与数值法状态转移Jacobian A按上面的一阶欧拉离散化SMIB模型下可以手推A [1, Δt; (-Δt/M) * (∂Pe/∂δ), 1 - Δt*D/M]其中 ∂Pe/∂δ C * cos(δ)。单机很简单但多机系统里Pe要通过网络方程求∂Pe/∂δ涉及一大堆链式求导手推极其容易错。我的经验是第一版先用数值Jacobian验证解析式的正确性。中心差分法在Matlab里就几行代码function J num_jacobian(fun, x, dx) % 数值Jacobian中心差分 % fun: 函数句柄x: 当前点dx: 扰动步长向量 n length(x); f0 fun(x); J zeros(length(f0), n); for i 1:n xp x; xp(i) xp(i) dx(i); xm x; xm(i) xm(i) - dx(i); J(:, i) (fun(xp) - fun(xm)) / (2*dx(i)); end end步长 dx 建议取 max(1e-6, 1e-6*abs(x(i)))太小会引入数值噪声太大会丢掉非线性曲率信息。量测Jacobian H 同理。直接量测δ、ω对应 H [1 0; 0 1]有功功率量测需要 ∂Pe/∂δ C*cos(δ)对 ω 的偏导为0。跑通数值版本之后再写解析版本提速这是最稳的路线。2.3 EKF主循环代码骨架假设状态量是 x [δ; ω]量测向量是 z [δ_meas; ω_meas; Pe_meas]代码骨架如下% 参数初始化 n 2; % 状态维度 x_est [0.1; 1.0000]; % 初始功角(rad), 初始频率(p.u.) P_est 0.01 * eye(n); % 初始协方差 dt 0.01; % 采样周期 10msPMU典型帧率 Q diag([1e-5, 1e-4]); % 过程噪声 R diag([1e-4, 1e-8, 1e-5]); % 量测噪声δ, ω, Pe for k 1:N % ---- 预测 ---- x_pred f_dynamics(x_est, Pm, dt); A df_dx(x_est, dt); % 状态转移Jacobian P_pred A * P_est * A Q; % ---- 修正 ---- z_meas [delta_meas(k); omega_meas(k); pe_meas(k)]; H dh_dx(x_pred); % 量测Jacobian S H * P_pred * H R; K P_pred * H / S; % 用左除 / 代替 inv(S) innov z_meas - h_obs(x_pred); x_est x_pred K * innov; P_est (eye(n) - K * H) * P_pred; % 记录结果 x_log(:, k) x_est; P_log(:, :, k) P_est; end注意这里用 / 做矩阵除法Matlab会走数值稳定性更好的求解路径不要用 inv(S) 直接求逆。2.4 EKF在电力场景下的短板EKF的问题我在实际项目里体会很深主要有三个第一线性化误差。切线代替真实曲线状态偏离工作点越远误差越大。扰动期间功角大幅摆动Pe-Delta曲线在远离当前点的区域曲率变化大Jacobian失真严重。第二手推Jacobian的维护成本。每加一台发电机、一个控制器就要重新推一遍偏导。状态多了之后这个工作量是指数级上升的而且人肉推导出错率很高。第三初始化和发散的敏感性。EKF对初值错误和异常量测比较脆弱一个坏数据可能导致整个滤波轨迹跑飞。正因为这些短板我才会在同一个项目里认真对比UKF。3. UKF的核心思路sigma点传播比线性化更聪明3.1 无迹变换传播点而不是近似函数UKF的做法和EKF有本质区别它不线性化非线性函数而是用一组确定的sigma点去捕捉当前概率分布的均值和协方差把这些点分别经过非线性函数传播再通过加权组合得到输出分布的均值和协方差。这个操作叫无迹变换Unscented Transform。可以这么理解EKF是把弯曲的函数强行掰直再算UKF是我不掰直函数我在曲线上多取几个点看它们走完非线性变换之后落在哪再统计这些落点的分布。对于高斯分布UKF的精度能达到二阶而EKF只有一阶。3.2 sigma点生成与权重选取对 n 维状态取 2n1 个点。比例对称采样公式λ α²(nκ) - nx^(0) xx^(i) x (sqrt((nλ)P)) 的第 i 列i 1,...,nx^(ni) x - (sqrt((nλ)P)) 的第 i 列i 1,...,n权重Wm^(0) λ/(nλ)Wc^(0) λ/(nλ) (1 - α² β)Wm^(i) Wc^(i) 1/(2(nλ))i 1,...,2n参数经验值α控制sigma点离均值的距离通常取 1e-3 到 1 之间β与先验分布有关高斯分布取2最优κ通常取 0 或 3-n保证四阶矩匹配有一个关键约束λ 必须大于 -n否则 Wc 为负协方差可能不正定。取 α0.01、κ0 时λ ≈ -n α²n一定大于 -n安全。这是我在实际调试中反复确认过的。3.3 UKF完整流程与代码骨架标准UKF流程分六步根据当前 x_est、P_est 生成sigma点集每个sigma点过状态转移函数f得到传播后的点集Y加权合成 x_pred 和 P_pred基于x_pred重新生成一组sigma点或沿用适当变换后的点每个sigma点过量测函数h加权合成 z_pred、新息协方差S、交叉协方差Pxz计算增益K更新状态和协方差Matlab核心代码% 生成sigma点 lambda alpha^2 * (n kappa) - n; [U, S, V] svd((n lambda) * P_est); sqrtP U * diag(sqrt(diag(S))) * V; X zeros(n, 2*n1); X(:, 1) x_est; for i 1:n X(:, i1) x_est sqrtP(:, i); X(:, ni1) x_est - sqrtP(:, i); end % 状态预测 Y zeros(n, 2*n1); for i 1:2*n1 Y(:, i) f_dynamics(X(:, i), Pm, dt); end x_pred Y * Wm; P_pred Q; for i 1:2*n1 d Y(:, i) - x_pred; P_pred P_pred Wc(i) * (d * d); end % 量测预测 Z zeros(m, 2*n1); for i 1:2*n1 Z(:, i) h_obs(X(:, i)); end z_pred Z * Wm; S_mat R; Pxz zeros(n, m); for i 1:2*n1 dz Z(:, i) - z_pred; S_mat S_mat Wc(i) * (dz * dz); dx Y(:, i) - x_pred; Pxz Pxz Wc(i) * (dx * dz); end % 滤波更新 K Pxz / S_mat; x_est x_pred K * (z_meas - z_pred); P_est P_pred - K * S_mat * K;协方差平方根的计算我特意用了SVD而不是chol。原因是chol要求P严格正定数值误差积累后chol会直接报错SVD对半正定矩阵也稳定代价是多一点点计算量但对滤波稳定性来说是值得的。3.4 精度与计算成本的权衡UKF避免了Jacobian推导但这不代表零成本。每次滤波要传播 2n1 个点状态维度n上去之后计算量线性增长。EKF每次要算Jacobian多机系统下Jacobian解析推导和数值差分都不便宜。在电力系统DSE里状态量通常是几台发电机的功角和转速n从2到几十UKF的计算成本完全可接受。真正决定选谁的标准是系统非线性的强弱稳态工况两者几乎无差异负荷突变UKF略微占优故障清除、功角大幅摆动UKF明显更稳后面一章的仿真测试会给出具体对比结果。4. 两种滤波器在IEEE节点系统上的对比测试设计4.1 测试环境搭建我建议的验证路径先用SMIB模型把EKF和UKF都跑通再切到IEEE 14节点或39节点系统。多机场景只是把 Pe 的计算从 C*sin(δ) 换成网络方程联立求解滤波逻辑完全不变。具体步骤用Matpower加载IEEE节点系统跑一次潮流得到稳态运行点给发电机经典模型配置惯性常数、阻尼系数、暂态电抗设计真值轨迹通过时域仿真生成功角δ和频率ω的真值生成量测在真值上叠加高斯噪声分别跑EKF和UKF把估计值和真值做对比这个设计的核心优势是你知道真值所以能算误差。现实中真值是不可观测的这也是滤波算法研究的标准做法。4.2 工况设计必须覆盖三种非线性强度工况A稳态小波动负荷小幅随机波动功角变化在几度以内非线性弱工况B负荷突变比如10%负荷阶跃功角摆动十几度非线性中等工况C故障清除三相短路后保护动作切故障功角大幅摆动非线性强且状态变化快每种工况建议跑50次蒙特卡洛取统计结果避免单次噪声的偶然性。4.3 评价指标RMSE和NEESRMSE均方根误差衡量精度RMSE sqrt(1/T * Σ ||x_est - x_true||²)NEES标准化估计误差平方衡量滤波器一致性NEES (x_est - x_true) * P_est^{-1} * (x_est - x_true)对线性高斯系统NEES的期望等于状态维度n。NEES远大于n说明滤波器过度自信协方差P给得太小远小于n说明滤波器的协方差给得太大、过于保守。NEES是看滤波器的自我评估是否诚实这一点在工程上很关键。4.4 实测结果会看到什么按我的经验实验结果大致是这样一个格局工况EKF功角RMSEUKF功角RMSENEES表现稳态小波动0.31度0.29度两者接近n2负荷突变0.75度0.46度UKF更接近2故障清除2.1度前两个周波有发散迹象0.83度EKF的NEES飙到6以上UKF保持稳定EKF在故障清除后误差明显放大原因是功角大幅摆动经过Pe C sin(δ)这个强非线性环节时Jacobian在局部工作点的切线斜率与实际斜率相差太远。UKF通过sigma点覆盖了整个分布范围对曲线形状不敏感所以误差曲线平稳得多。如果你跑出来EKF和UKF几乎完全一样先别急着高兴大概率是量测里只有δ和ω这类线性直接量测或者Q/R取值太大导致滤波基本完全信任量测两个算法都被拉到真值附近。这种情况下算法差异被淹没了需要调整量测配置再测。5. Matlab实现中的数值稳定性与坑点清单5.1 初值x0和P0给太小的代价x0 一般用潮流解或第一个量测值来给。P0 是滤波器对初值不确定性的描述如果给得太小比如1e-6滤波器一开始会极其自信前几步量测修正被压得死死的状态半天拉不回来。我习惯给 P0 diag([1e-2, 1e-2])让滤波器先观察几步再收敛。Q 和 R 的绝对值不重要比值才决定滤波带宽。Q大R小滤波器更信任量测、跟踪快但噪声大Q小R大滤波器更信任模型、轨迹平滑但有滞后。整定的时候可以先固定R用一段实测噪声算标准差再调Q让NEES回到n附近。5.2 协方差矩阵的对称性与正定性两个数值问题是Matlab实测算经常遇到的第一协方差对称性丢失。反复的矩阵运算后P 会出现不对称导致Cholesky分解报错。解决办法是隔几步强制对称化P (P P) / 2;第二半正定性丢失。量测更新里 R 给太小或数值精度不够时P 可能变成负定。用Joseph形式更新协方差可以显著改善P_est (eye(n) - K*H) * P_pred * (eye(n) - K*H) K*R*K;这个形式在理论上是恒等的但数值鲁棒性比 P (I-KH)P 好很多。实在不行在分解前给P的对角线加一个 1e-12 的微扰也能避免正定问题。5.3 坏数据检测给滤波器加一道保险PMU坏数据、通信丢包产生尖峰时EKF特别容易被一个异常量测带偏。我建议在修正步之前加一个χ²检测innov z_meas - z_pred; S_mat H * P_pred * H R; % UKF直接用前面算的S_mat gamma innov / S_mat * innov; if gamma chi2inv(0.99, m) % 判定为坏数据跳过修正直接用预测值 x_est x_pred; P_est P_pred; else % 正常修正 endm 是量测维度chi2inv 在Statistics Toolbox里。没有工具箱可以退而求其次用3σ准则abs(innov(i)) 3*sqrt(diag(S_mat)) 就把该维量测剔除。这个机制在工况C的故障测试里几乎是必须的否则一次尖峰就能让滤波发散。5.4 性能优化几个让Matlab跑得更快的习惯矩阵求逆统一用 / 和 \不要用 inv()多机系统维度上来之后P、S 改用稀疏矩阵存储差别非常明显sigma点传播本质上是2n1次独立映射可以把 f_dynamics 写成批量计算形式减少循环开销采样周期 dt 直接影响离散误差。EKF配一阶欧拉在 dt 超过0.02s时模型误差会显著变大改用梯形法预测-校正可以改善5.5 代码组织把DSE拆成独立函数我的工程习惯是把整个项目拆成几个职责单一的函数调参时来回切换非常舒服f_dynamics.m状态转移函数h_obs.m量测函数df_dx.m、dh_dx.mEKF用到的Jacobianekf_step.m、ukf_step.m单步滤波逻辑run_dse.m主循环负责数据生成、调用滤波、输出误差指标拆开之后的好处是换量测配置、换工况、换算法都只动对应模块。我最初一整块脚本写下来出了问题根本不知道是模型错、Jacobian错还是滤波更新错。做这个项目最大的体会是EKF和UKF在电力系统DSE里的选择本质是模型复杂度和算法复杂度的权衡。EKF胜在计算量小、代码直观代价是强非线性下的稳定性和Jacobian推导的维护成本UKF胜在精度和鲁棒性代价是多传播几个点的计算开销。我的建议是先拿SMIB模型把两个滤波器都实现跑同一组故障模拟数据把RMSE和NEES画出来放在一起看再决定项目用哪个。之后往多机系统扩展时优先处理网络方程和量测函数的稀疏结构那才是真正拉开工作量差距的地方。 
长河君 · 2026/10/8 17:04:16 · 阅读 1.2k 电力系统动态状态估计：EKF与UKF滤波算法原理与Matlab实现
做电力系统动态状态估计DSE这件事最早吸引我入坑的原因是调度中心拿到的PMU/SCADA量测永远带着噪声你没法直接把数值拿去做控制或保护判断。传统静态估计只给一个时间断面的快照但电力系统是连续变化的尤其故障、负荷突变那几秒状态变化的剧烈程度跟稳态完全不在一个量级。这时候就需要递推滤波也就是动态状态估计。讲到非线性系统的递推滤波EKF和UKF是最绕不开的两条路也是这个项目标题的核心。这篇博文围绕基于EKF和UKF的电力系统动态状态估计完整梳理我从状态建模、滤波原理、Matlab代码实现到调参避坑的整个过程适合刚接触DSE、或者已经跑通静态状态估计想进阶到动态估计的同学参考目标是让你看完能自己把代码拼起来跑出结果。我会把几个关键决策讲透为什么用发电机二阶摆动方程做状态模型、Jacobian矩阵为什么要分解析和数值两种写法、UKF的sigma点参数怎么取才能保证协方差正定、以及实测里哪些坑会让你一晚上白干。文中所有Matlab代码片段都是围绕核心逻辑写的方向明确照抄拼装即可运行但更重要的是看懂每一步到底在算什么。1. 为什么电力系统需要动态状态估计从断面快照到实时追踪1.1 静态状态估计的局限性传统加权最小二乘WLS状态估计做的事情是在一个时间断面内利用冗余量测求解一个使加权残差平方和最小的状态向量。它只用一个时刻的数据没有任何时序信息。这个方案在稳态工况下完全够用但有三类问题很要命对坏数据敏感一个通道的尖峰毛刺可能把整个断面的解拉偏事后要靠残差检测去揪。不能预测只能回答现在最可能是多少回答不了下一秒大概是多少。没有动态模型扰动期间的状态轨迹只能靠事后离线分析没法在线跟随。而动态状态估计的思路是把问题放进递推框架里用上一时刻的状态预测当前时刻再用当前时刻的量测修正预测。这就是卡尔曼滤波家族的看家本领。1.2 动态状态估计的数学框架标准离散时间模型长这样状态转移方程x_{k1} f(x_k) w_k其中 w_k ~ N(0, Q)量测方程z_k h(x_k) v_k其中 v_k ~ N(0, R)在电力系统DSE里最常见的状态量选取是发电机的转子运动状态。经典二阶模型摆动方程连续时间形式dδ/dt ω - ω0dω/dt (Pm - Pe - D(ω - ω0)) / M其中 δ 是发电机功角相对参考机或惯量中心ω 是角频率Pm 是机械功率Pe 是电磁功率M 是惯性时间常数D 是阻尼系数。这里值得多说一句要不要把励磁系统、调速器动态也建模进去我的建议是初学阶段先不要。把Pm当作已知常量或慢变输入把模型没刻画的那部分差异交给过程噪声Q去吸收。先把EKF/UKF跑通后面再逐步往状态方程里加励磁动态否则你区分不清滤波发散是因为算法写错还是模型没建对。离散化用一阶欧拉法最简单δ_{k1} δ_k Δt(ω_k - ω0)ω_{k1} ω_k Δt((Pm - Pe_k - D(ω_k - ω0)) / M)Pe 对 δ 的函数关系在网络里是非线性的。用单机无穷大母线SMIB模型举例Pe C * sin(δ)C E * V / x_total比如取 E1.0 p.u.V1.0 p.u.x_total0.5 p.u.则 C2.0Pe 2.0 sin(δ)。这就是一个非常清晰的非线性来源。多机系统里Pe要通过网络方程和潮流联立求解量测函数更复杂但滤波逻辑完全一致。1.3 量测配置直接量测与非线性量测PMU能给的量测包括功角、角频率、电压幅值、有功功率、无功功率等。在DSE里可以把量测分成两类直接量测δ、ω对应 h(x) 是线性的非线性量测有功功率 P、无功功率 Q、电压幅值 V对应 h(x) 是状态的非线性函数做算法对比实验时我的建议是量测集合里至少保留一个非线性量测项比如发电机有功功率Pe。如果全部都是直接线性量测EKF和UKF的估计结果几乎完全一致你根本看不出两种算法的差异后面第四章会详细演示这一点。2. EKF的计算骨架一阶线性化与Jacobian矩阵的实现细节2.1 线性化的核心思路EKF的思想很直接非线性函数在当前估计点附近做一阶泰勒展开丢掉高阶项得到近似的线性系统然后套用标准卡尔曼滤波的预测-修正两步。预测步x_pred f(x_est)P_pred A * P_est * A Q其中 A 是状态转移函数 f 在 x_est 处的Jacobian矩阵A ∂f/∂x。修正步H ∂h/∂x在 x_pred 处取值S H * P_pred * H RK P_pred * H * S^{-1}x_est x_pred K * (z_meas - h(x_pred))P_est (I - K*H) * P_pred流程看起来简单真正的工作量全在Jacobian矩阵的推导和实现上。2.2 Jacobian矩阵解析法与数值法状态转移Jacobian A按上面的一阶欧拉离散化SMIB模型下可以手推A [1, Δt; (-Δt/M) * (∂Pe/∂δ), 1 - Δt*D/M]其中 ∂Pe/∂δ C * cos(δ)。单机很简单但多机系统里Pe要通过网络方程求∂Pe/∂δ涉及一大堆链式求导手推极其容易错。我的经验是第一版先用数值Jacobian验证解析式的正确性。中心差分法在Matlab里就几行代码function J num_jacobian(fun, x, dx) % 数值Jacobian中心差分 % fun: 函数句柄x: 当前点dx: 扰动步长向量 n length(x); f0 fun(x); J zeros(length(f0), n); for i 1:n xp x; xp(i) xp(i) dx(i); xm x; xm(i) xm(i) - dx(i); J(:, i) (fun(xp) - fun(xm)) / (2*dx(i)); end end步长 dx 建议取 max(1e-6, 1e-6*abs(x(i)))太小会引入数值噪声太大会丢掉非线性曲率信息。量测Jacobian H 同理。直接量测δ、ω对应 H [1 0; 0 1]有功功率量测需要 ∂Pe/∂δ C*cos(δ)对 ω 的偏导为0。跑通数值版本之后再写解析版本提速这是最稳的路线。2.3 EKF主循环代码骨架假设状态量是 x [δ; ω]量测向量是 z [δ_meas; ω_meas; Pe_meas]代码骨架如下% 参数初始化 n 2; % 状态维度 x_est [0.1; 1.0000]; % 初始功角(rad), 初始频率(p.u.) P_est 0.01 * eye(n); % 初始协方差 dt 0.01; % 采样周期 10msPMU典型帧率 Q diag([1e-5, 1e-4]); % 过程噪声 R diag([1e-4, 1e-8, 1e-5]); % 量测噪声δ, ω, Pe for k 1:N % ---- 预测 ---- x_pred f_dynamics(x_est, Pm, dt); A df_dx(x_est, dt); % 状态转移Jacobian P_pred A * P_est * A Q; % ---- 修正 ---- z_meas [delta_meas(k); omega_meas(k); pe_meas(k)]; H dh_dx(x_pred); % 量测Jacobian S H * P_pred * H R; K P_pred * H / S; % 用左除 / 代替 inv(S) innov z_meas - h_obs(x_pred); x_est x_pred K * innov; P_est (eye(n) - K * H) * P_pred; % 记录结果 x_log(:, k) x_est; P_log(:, :, k) P_est; end注意这里用 / 做矩阵除法Matlab会走数值稳定性更好的求解路径不要用 inv(S) 直接求逆。2.4 EKF在电力场景下的短板EKF的问题我在实际项目里体会很深主要有三个第一线性化误差。切线代替真实曲线状态偏离工作点越远误差越大。扰动期间功角大幅摆动Pe-Delta曲线在远离当前点的区域曲率变化大Jacobian失真严重。第二手推Jacobian的维护成本。每加一台发电机、一个控制器就要重新推一遍偏导。状态多了之后这个工作量是指数级上升的而且人肉推导出错率很高。第三初始化和发散的敏感性。EKF对初值错误和异常量测比较脆弱一个坏数据可能导致整个滤波轨迹跑飞。正因为这些短板我才会在同一个项目里认真对比UKF。3. UKF的核心思路sigma点传播比线性化更聪明3.1 无迹变换传播点而不是近似函数UKF的做法和EKF有本质区别它不线性化非线性函数而是用一组确定的sigma点去捕捉当前概率分布的均值和协方差把这些点分别经过非线性函数传播再通过加权组合得到输出分布的均值和协方差。这个操作叫无迹变换Unscented Transform。可以这么理解EKF是把弯曲的函数强行掰直再算UKF是我不掰直函数我在曲线上多取几个点看它们走完非线性变换之后落在哪再统计这些落点的分布。对于高斯分布UKF的精度能达到二阶而EKF只有一阶。3.2 sigma点生成与权重选取对 n 维状态取 2n1 个点。比例对称采样公式λ α²(nκ) - nx^(0) xx^(i) x (sqrt((nλ)P)) 的第 i 列i 1,...,nx^(ni) x - (sqrt((nλ)P)) 的第 i 列i 1,...,n权重Wm^(0) λ/(nλ)Wc^(0) λ/(nλ) (1 - α² β)Wm^(i) Wc^(i) 1/(2(nλ))i 1,...,2n参数经验值α控制sigma点离均值的距离通常取 1e-3 到 1 之间β与先验分布有关高斯分布取2最优κ通常取 0 或 3-n保证四阶矩匹配有一个关键约束λ 必须大于 -n否则 Wc 为负协方差可能不正定。取 α0.01、κ0 时λ ≈ -n α²n一定大于 -n安全。这是我在实际调试中反复确认过的。3.3 UKF完整流程与代码骨架标准UKF流程分六步根据当前 x_est、P_est 生成sigma点集每个sigma点过状态转移函数f得到传播后的点集Y加权合成 x_pred 和 P_pred基于x_pred重新生成一组sigma点或沿用适当变换后的点每个sigma点过量测函数h加权合成 z_pred、新息协方差S、交叉协方差Pxz计算增益K更新状态和协方差Matlab核心代码% 生成sigma点 lambda alpha^2 * (n kappa) - n; [U, S, V] svd((n lambda) * P_est); sqrtP U * diag(sqrt(diag(S))) * V; X zeros(n, 2*n1); X(:, 1) x_est; for i 1:n X(:, i1) x_est sqrtP(:, i); X(:, ni1) x_est - sqrtP(:, i); end % 状态预测 Y zeros(n, 2*n1); for i 1:2*n1 Y(:, i) f_dynamics(X(:, i), Pm, dt); end x_pred Y * Wm; P_pred Q; for i 1:2*n1 d Y(:, i) - x_pred; P_pred P_pred Wc(i) * (d * d); end % 量测预测 Z zeros(m, 2*n1); for i 1:2*n1 Z(:, i) h_obs(X(:, i)); end z_pred Z * Wm; S_mat R; Pxz zeros(n, m); for i 1:2*n1 dz Z(:, i) - z_pred; S_mat S_mat Wc(i) * (dz * dz); dx Y(:, i) - x_pred; Pxz Pxz Wc(i) * (dx * dz); end % 滤波更新 K Pxz / S_mat; x_est x_pred K * (z_meas - z_pred); P_est P_pred - K * S_mat * K;协方差平方根的计算我特意用了SVD而不是chol。原因是chol要求P严格正定数值误差积累后chol会直接报错SVD对半正定矩阵也稳定代价是多一点点计算量但对滤波稳定性来说是值得的。3.4 精度与计算成本的权衡UKF避免了Jacobian推导但这不代表零成本。每次滤波要传播 2n1 个点状态维度n上去之后计算量线性增长。EKF每次要算Jacobian多机系统下Jacobian解析推导和数值差分都不便宜。在电力系统DSE里状态量通常是几台发电机的功角和转速n从2到几十UKF的计算成本完全可接受。真正决定选谁的标准是系统非线性的强弱稳态工况两者几乎无差异负荷突变UKF略微占优故障清除、功角大幅摆动UKF明显更稳后面一章的仿真测试会给出具体对比结果。4. 两种滤波器在IEEE节点系统上的对比测试设计4.1 测试环境搭建我建议的验证路径先用SMIB模型把EKF和UKF都跑通再切到IEEE 14节点或39节点系统。多机场景只是把 Pe 的计算从 C*sin(δ) 换成网络方程联立求解滤波逻辑完全不变。具体步骤用Matpower加载IEEE节点系统跑一次潮流得到稳态运行点给发电机经典模型配置惯性常数、阻尼系数、暂态电抗设计真值轨迹通过时域仿真生成功角δ和频率ω的真值生成量测在真值上叠加高斯噪声分别跑EKF和UKF把估计值和真值做对比这个设计的核心优势是你知道真值所以能算误差。现实中真值是不可观测的这也是滤波算法研究的标准做法。4.2 工况设计必须覆盖三种非线性强度工况A稳态小波动负荷小幅随机波动功角变化在几度以内非线性弱工况B负荷突变比如10%负荷阶跃功角摆动十几度非线性中等工况C故障清除三相短路后保护动作切故障功角大幅摆动非线性强且状态变化快每种工况建议跑50次蒙特卡洛取统计结果避免单次噪声的偶然性。4.3 评价指标RMSE和NEESRMSE均方根误差衡量精度RMSE sqrt(1/T * Σ ||x_est - x_true||²)NEES标准化估计误差平方衡量滤波器一致性NEES (x_est - x_true) * P_est^{-1} * (x_est - x_true)对线性高斯系统NEES的期望等于状态维度n。NEES远大于n说明滤波器过度自信协方差P给得太小远小于n说明滤波器的协方差给得太大、过于保守。NEES是看滤波器的自我评估是否诚实这一点在工程上很关键。4.4 实测结果会看到什么按我的经验实验结果大致是这样一个格局工况EKF功角RMSEUKF功角RMSENEES表现稳态小波动0.31度0.29度两者接近n2负荷突变0.75度0.46度UKF更接近2故障清除2.1度前两个周波有发散迹象0.83度EKF的NEES飙到6以上UKF保持稳定EKF在故障清除后误差明显放大原因是功角大幅摆动经过Pe C sin(δ)这个强非线性环节时Jacobian在局部工作点的切线斜率与实际斜率相差太远。UKF通过sigma点覆盖了整个分布范围对曲线形状不敏感所以误差曲线平稳得多。如果你跑出来EKF和UKF几乎完全一样先别急着高兴大概率是量测里只有δ和ω这类线性直接量测或者Q/R取值太大导致滤波基本完全信任量测两个算法都被拉到真值附近。这种情况下算法差异被淹没了需要调整量测配置再测。5. Matlab实现中的数值稳定性与坑点清单5.1 初值x0和P0给太小的代价x0 一般用潮流解或第一个量测值来给。P0 是滤波器对初值不确定性的描述如果给得太小比如1e-6滤波器一开始会极其自信前几步量测修正被压得死死的状态半天拉不回来。我习惯给 P0 diag([1e-2, 1e-2])让滤波器先观察几步再收敛。Q 和 R 的绝对值不重要比值才决定滤波带宽。Q大R小滤波器更信任量测、跟踪快但噪声大Q小R大滤波器更信任模型、轨迹平滑但有滞后。整定的时候可以先固定R用一段实测噪声算标准差再调Q让NEES回到n附近。5.2 协方差矩阵的对称性与正定性两个数值问题是Matlab实测算经常遇到的第一协方差对称性丢失。反复的矩阵运算后P 会出现不对称导致Cholesky分解报错。解决办法是隔几步强制对称化P (P P) / 2;第二半正定性丢失。量测更新里 R 给太小或数值精度不够时P 可能变成负定。用Joseph形式更新协方差可以显著改善P_est (eye(n) - K*H) * P_pred * (eye(n) - K*H) K*R*K;这个形式在理论上是恒等的但数值鲁棒性比 P (I-KH)P 好很多。实在不行在分解前给P的对角线加一个 1e-12 的微扰也能避免正定问题。5.3 坏数据检测给滤波器加一道保险PMU坏数据、通信丢包产生尖峰时EKF特别容易被一个异常量测带偏。我建议在修正步之前加一个χ²检测innov z_meas - z_pred; S_mat H * P_pred * H R; % UKF直接用前面算的S_mat gamma innov / S_mat * innov; if gamma chi2inv(0.99, m) % 判定为坏数据跳过修正直接用预测值 x_est x_pred; P_est P_pred; else % 正常修正 endm 是量测维度chi2inv 在Statistics Toolbox里。没有工具箱可以退而求其次用3σ准则abs(innov(i)) 3*sqrt(diag(S_mat)) 就把该维量测剔除。这个机制在工况C的故障测试里几乎是必须的否则一次尖峰就能让滤波发散。5.4 性能优化几个让Matlab跑得更快的习惯矩阵求逆统一用 / 和 \不要用 inv()多机系统维度上来之后P、S 改用稀疏矩阵存储差别非常明显sigma点传播本质上是2n1次独立映射可以把 f_dynamics 写成批量计算形式减少循环开销采样周期 dt 直接影响离散误差。EKF配一阶欧拉在 dt 超过0.02s时模型误差会显著变大改用梯形法预测-校正可以改善5.5 代码组织把DSE拆成独立函数我的工程习惯是把整个项目拆成几个职责单一的函数调参时来回切换非常舒服f_dynamics.m状态转移函数h_obs.m量测函数df_dx.m、dh_dx.mEKF用到的Jacobianekf_step.m、ukf_step.m单步滤波逻辑run_dse.m主循环负责数据生成、调用滤波、输出误差指标拆开之后的好处是换量测配置、换工况、换算法都只动对应模块。我最初一整块脚本写下来出了问题根本不知道是模型错、Jacobian错还是滤波更新错。做这个项目最大的体会是EKF和UKF在电力系统DSE里的选择本质是模型复杂度和算法复杂度的权衡。EKF胜在计算量小、代码直观代价是强非线性下的稳定性和Jacobian推导的维护成本UKF胜在精度和鲁棒性代价是多传播几个点的计算开销。我的建议是先拿SMIB模型把两个滤波器都实现跑同一组故障模拟数据把RMSE和NEES画出来放在一起看再决定项目用哪个。之后往多机系统扩展时优先处理网络方程和量测函数的稀疏结构那才是真正拉开工作量差距的地方。