尧图网页设计 
最近几年做电力系统动态状态估计的人越来越多了PMU逐步普及之后大家发现传统的静态状态估计加权最小二乘那一套已经喂不饱实时监控的需求。动态状态估计的核心就是利用系统的动态模型发电机转子运动方程这类和量测数据实时递推得到系统状态。而在这个场景里EKF和UKF是绕不开的两个名字——一个靠泰勒展开线性化一个靠sigma点采样思路完全不同但目的都是把非线性模型塞进卡尔曼滤波的框架里。这篇博文就把这件事彻底讲透。我会从为什么要做动态状态估计讲起把EKF和UKF的原理差异拆开对比然后给出完整的电力系统状态空间模型最后用MATLAB把两者实现出来对比仿真结果再把我踩过的坑和调试经验一并分享出来。代码结构我会写成可以直接改参数复用的形式方便你在自己的研究或者课程设计里快速上手。1. 为什么需要动态状态估计静态那套已经不够用了1.1 静态状态估计的局限性传统电力系统状态估计不管是加权最小二乘WLS还是快速分解法本质都是在某个时间断面上根据SCADA量测电压幅值、注入功率、支路潮流去求解一组非线性代数方程得到系统的电压幅值和相角。这套方法在调度中心跑了三四十年稳定性没有问题但它的前提是系统处于稳态或准稳态量测刷新周期是秒级甚至分钟级。问题来了。当系统发生扰动、振荡或者故障后的恢复过程时状态量功角、转速是在毫秒到秒级别快速变化的。SCADA的采样率根本跟不上而且静态估计把时间上相邻的断面当成独立问题处理完全扔掉了系统动态方程里蕴含的信息。这就像你开车只看仪表盘照片还是每隔几秒才拍一张根本没法判断下一秒会不会失控。动态状态估计的思路完全不一样——它把发电机转子运动方程等动态模型作为状态转移方程把PMU量测同步相量、频率变化率作为观测方程用卡尔曼滤波家族的算法逐拍递推实时估计功角、转速这类动态状态量。PMU的采样率通常是每秒30~60帧这为动态估计提供了充足的量测基础。1.2 为什么偏偏是卡尔曼滤波家族状态估计问题本质上就是“带噪声的最优滤波问题”卡尔曼滤波提供了线性高斯条件下的最优解。但电力系统的量测方程和动态方程几乎都是非线性的功角的正弦项在输出方程里天然存在发电机的电功率公式里也有非线性耦合。面对非线性系统工程上基本就是两条路近似解析EKF和近似采样UKF。如果你用粒子滤波虽然理论上能处理任意非线性非高斯但计算量在实时场景下很难扛住。所以EKF和UKF成了电力系统动态状态估计里最主流的两个选择。前者实现简单、计算量小后者精度更高、不需要求雅可比矩阵两者之间的取舍其实很有意思后面专门展开聊。2. EKF和UKF的核心原理解析线性化与统计线性化的分岔路2.1 EKF的思路泰勒展开打到一阶EKF的思想很直白既然系统是非线性的我就在当前状态估计点附近把状态转移函数和量测函数做一阶泰勒展开取雅可比矩阵作为线性化后的状态转移矩阵和量测矩阵然后直接套用标准卡尔曼滤波的预测—更新框架。以一发电机系统为例设状态向量为$x [\delta, \omega]^T$动态方程为$$x_{k1} f(x_k, u_k) w_k$$量测方程为$$z_k h(x_k) v_k$$其中$w_k \sim N(0, Q)$是过程噪声$v_k \sim N(0, R)$是量测噪声。EKF的递推分为两步。预测步$$\hat{x}{k|k-1} f(\hat{x}{k-1|k-1}, u_{k-1})$$$$P_{k|k-1} F_k P_{k-1|k-1} F_k^T Q$$其中$F_k \left. \frac{\partial f}{\partial x} \right|{\hat{x}{k-1|k-1}}$。更新步$$K_k P_{k|k-1} H_k^T \left( H_k P_{k|k-1} H_k^T R \right)^{-1}$$$$\hat{x}{k|k} \hat{x}{k|k-1} K_k \left( z_k - h(\hat{x}_{k|k-1}) \right)$$$$P_{k|k} \left( I - K_k H_k \right) P_{k|k-1}$$这个流程看着简单但实际用起来有两个痛点。第一雅可比矩阵$F_k$和$H_k$的推导极其容易出错尤其是多机系统里状态变量多、方程耦合强的时候一个符号错了整套滤波直接发散第二一阶线性化在系统非线性强度高的时候误差很大特别是功角摆动幅度大的暂态过程中线性化误差会直接导致估计偏差甚至发散。2.2 UKF的思路sigma点告诉你什么是“无迹”UKF的出发点完全不同。它不再去做解析线性化而是利用“对非线性函数传播概率分布比直接对非线性函数做线性近似更容易”这个思想。具体做法是在当前状态均值附近按照一定的规则选取一组sigma点让这些点的均值和协方差与原分布一致然后把这些点逐个通过非线性函数再用传递后的点加权重建均值和协方差。sigma点选取规则最常用的是对称策略。设状态维度为$n$取$2n1$个点$$\chi^{(0)} \bar{x}, \quad W_m^{(0)} \frac{\lambda}{n\lambda}$$$$\chi^{(i)} \bar{x} \left( \sqrt{(n\lambda)P} \right)_i, \quad W_m^{(i)} \frac{1}{2(n\lambda)}$$$$\chi^{(in)} \bar{x} - \left( \sqrt{(n\lambda)P} \right)_i, \quad W_c^{(i)} \frac{1}{2(n\lambda)}$$其中$\lambda \alpha^2(n\kappa) - n$$\alpha$控制sigma点的散布范围通常取$10^{-3}$到$1$$\kappa$是次级缩放参数通常取$0$或$3-n$。预测时把每个sigma点代入非线性状态转移函数$$\chi^{(i)}{k|k-1} f\left( \chi^{(i)}{k-1|k-1}, u_{k-1} \right)$$加权得到预测均值和协方差$$\hat{x}{k|k-1} \sum{i0}^{2n} W_m^{(i)} \chi^{(i)}_{k|k-1}$$$$P_{k|k-1} \sum_{i0}^{2n} W_c^{(i)} \left( \chi^{(i)}{k|k-1} - \hat{x}{k|k-1} \right) \left( \cdots \right)^T Q$$更新步类似把sigma点通过量测函数$h$得到量测预测均值和协方差以及状态与量测的互协方差再算出卡尔曼增益。整个过程不需要计算任何雅可比矩阵所以UKF也被称为“免求导”滤波方法。2.3 两者对比你真的需要UKF吗对比项EKFUKF核心手段一阶泰勒展开求雅可比sigma点采样统计近似实现难度需要手推雅可比易错不需要求导实现统一非线性程度弱非线性时精度可观强非线性下精度远高于EKF计算量每步一次雅可比一次协方差传播每步$2n1$次非线性函数传播数值稳定性强非线性时易发散相对稳健但需保证$P$正定适用范围弱非线性系统首选强非线性或雅可比难推导的系统从这张表可以看出一个关键结论如果你的系统非线性不强、雅可比好推导EKF就是够用且高效的选择但电力系统的PMU量测方程里充满正弦余弦项暂态过程中功角摆开幅度大我实测下来UKF的估计精度普遍比EKF高一个量级左右尤其是对转速$\omega$的估计。代价是多算了$2n$次非线性函数传播但对单机系统$n2$来说多算4次函数几乎可以忽略。提示在多机系统中状态维度$n$会随发电机台数线性增长UKF每拍要传播$2n1$个sigma点计算量优势会被稀释。实际工程中常把UKF用于单机或小规模系统验证大规模系统可以考虑事后估计精度和算力的平衡。3. 状态空间建模把发电机摆动方程写成滤波器能用的形式3.1 状态量的选择功角和转速是底线电力系统动态状态估计中最经典的状态向量是发电机的功角$\delta$和转速偏差$\omega - \omega_s$。功角决定了发电机转子之间的相对位置转速偏差决定了系统的频率动态。更精细的模型还会把暂态电动势$E_q$加入状态向量形成三阶或四阶模型但基础的两阶摇摆方程足以讲清楚EKF和UKF的实现逻辑也便于验证算法效果。我采用最常用的单机无穷大母线SMIB系统作为仿真对象。这样做的好处是物理概念清晰、解析模型简单而且结果容易和直接数值仿真对照验证。扩展到多机系统时只需要把每台发电机的状态拼接成大状态向量动态方程变成每台机各写一组摇摆方程量测方程叠加PMU量测即可算法框架完全不变。3.2 动态方程状态转移函数SMIB系统中发电机的转子运动方程摇摆方程为$$\frac{d\delta}{dt} \omega - \omega_s$$$$\frac{d\omega}{dt} \frac{1}{2H} \left( P_m - P_e - D(\omega - \omega_s) \right)$$其中$H$是发电机的惯性时间常数秒$D$是阻尼系数$P_m$是机械功率在仿真中通常假设恒定$P_e$是电磁功率。而电磁功率本身就是功角的非线性函数$$P_e \frac{E V_\infty}{X_\Sigma} \sin\delta$$这里$E$是发电机暂态电动势$V_\infty$是无穷大母线电压幅值$X_\Sigma$是发电机暂态电抗与线路电抗之和。把这个连续方程离散化用最简单的向前欧拉法$$\delta_{k1} \delta_k \Delta t \cdot (\omega_k - \omega_s)$$$$\omega_{k1} \omega_k \frac{\Delta t}{2H} \left( P_m - \frac{E V_\infty}{X_\Sigma} \sin\delta_k - D(\omega_k - \omega_s) \right)$$注意这里$\Delta t$的选择很关键。PMU量测时间戳通常是0.02s或0.033s对应50Hz/60Hz系统的一个或两个周波$\Delta t$直接取量测间隔即可。欧拉法在$\Delta t0.02s$时精度足够但如果你把$\Delta t$取得过大比如0.1s以上积分误差会明显影响滤波效果这时候建议换成四阶龙格-库塔法离散。3.3 量测方程PMU给了什么就用什么PMU可以直接提供节点的电压相量幅值和相角也可以通过相量推算线路电流和功率。这里采用最直接的量测发电机母线电压幅值$V_t$和电压相角$\theta_t$以及母线注入功率$P_e$和$Q_e$。以电压相角为例它和功角之间有明确的机电关系$$\theta_t \delta - \arctan\left( \frac{X_d P_e}{V_t^2 X_d Q_e} \right)$$公式里$X_d$是发电机暂态电抗。为了演示简洁我直接采用更简单的量测模型$$z h(x) \begin{bmatrix} P_e \ \theta_t \ V_t \end{bmatrix} \begin{bmatrix} \frac{E V_\infty}{X_\Sigma} \sin\delta \ \delta \ V_t(\delta) \end{bmatrix}$$其中电压幅值$V_t$也可以用功角的非线性表达式给出。实际上量测方程取什么形式取决于你能拿到哪些PMU量测。核心原则是量测方程必须能写成状态量的可微或可直接采样传播的函数。3.4 噪声矩阵Q和R的设置滤波器的命门过程噪声协方差矩阵$Q$和量测噪声协方差矩阵$R$的设置直接决定滤波器的性能。$Q$描述的是模型误差和未建模动态的统计特性$R$描述的是PMU量测的误差水平。我的经验值PMU幅值量测噪声标准差大约为0.5%~1%标幺值计相角量测噪声标准差约为0.2~0.5度约0.0035~0.0087弧度。如果量测是$z[P_e, \theta_t]^T$那么$$R \begin{bmatrix} (0.01)^2 0 \ 0 (0.005)^2 \end{bmatrix}$$$Q$的取法更需要细心。$Q$太小滤波器过于信任模型量测稍有偏差就更新不充分$Q$太大滤波器过度信任量测噪声抑制能力变差。常见的做法是把$Q$设成对角阵对角线元素取状态量物理范围的小百分比平方。例如功角$\delta$的量程是弧度级别取$\sigma_\delta 0.01$ rad则$Q_{11} 0.0001$转速偏差$\omega-\omega_s$的量程是标幺值取$\sigma_\omega 0.005$则$Q_{22} 2.5 \times 10^{-5}$。提示在仿真调试阶段我习惯先设较大的$Q$比如$10^{-3}$级别观察滤波器是否发散再逐步减小。如果滤波器出现“看似收敛但估计值有明显偏置”的情况通常是$Q$设得过小模型误差没有被充分补偿。4. MATLAB代码实现从零搭起EKF和UKF4.1 代码结构总览模块化才是王道我用MATLAB实现这套代码时刻意把它拆成了几个独立模块方便单独调试和复用文件/函数作用system_params.m定义系统参数结构体H、D、Xd、Vinf等f_dynamics.m状态转移函数离散化后的摇摆方程h_measurement.m量测函数ekf_filter.mEKF完整滤波函数ukf_filter.mUKF完整滤波函数generate_true_trajectory.m生成真实状态轨迹和量测数据main_demo.m主脚本串联整个仿真流程这么拆的好处是你要换系统参数、换量测配置、替换滤波器算法都只需要改动对应的一个模块而不是在一坨代码里翻找。4.2 系统参数与真实轨迹生成首先定义系统参数% system_params.m function params system_params() params.H 5.0; % 惯性时间常数 (s) params.D 2.0; % 阻尼系数 (p.u.) params.Xd 0.3; % 发电机暂态电抗线路电抗之和 (p.u.) params.E_prime 1.05; % 暂态电动势 (p.u.) params.V_inf 1.0; % 无穷大母线电压 (p.u.) params.omega_s 1.0; % 同步转速 (p.u.) params.Pm 0.8; % 机械功率 (p.u.) params.dt 0.02; % 采样间隔 (s) end真实轨迹生成用了仿真加噪声的方式。用一个较大的机械功率扰动作为激励让系统发生功角摆动然后用四阶龙格-库塔精确积分生成“真实”状态轨迹再叠加上高斯白噪声模拟量测% generate_true_trajectory.m function [x_true, z_meas, t] generate_true_trajectory(params, T_sim) dt params.dt; N round(T_sim / dt); t (0:N-1) * dt; % 初始状态功角在平衡点附近转速等于同步转速 delta0 asin(params.Pm * params.Xd / (params.E_prime * params.V_inf)); x_true zeros(2, N); x_true(:,1) [delta0; params.omega_s]; % 在第100采样时刻施加一个机械功率阶跃扰动 for k 1:N-1 Pm_k params.Pm; if k 100 Pm_k params.Pm * 1.2; % 20%机械功率阶跃 end x_true(:, k1) rk4_step(x_true(:,k), Pm_k, params); end % 生成量测电磁功率电压相角 z_true zeros(2, N); for k 1:N z_true(:, k) h_measurement(x_true(:,k), params); end % 叠加上量测噪声 R_diag [0.01^2; 0.005^2]; z_meas z_true sqrt(R_diag) .* randn(2, N); end上面用到的$rk4_step$函数就是标准的四阶龙格-库塔单步积分代码不复杂关键是让真实轨迹足够精确、让后续的状态估计有个可信的对照基准。4.3 EKF的两个雅可比手推还是数值差分EKF实现的精髓在雅可比矩阵。我建议在原理验证阶段用手推解析雅可比这样代码运行效率高、结果可信但如果你的系统复杂到雅可比推导困难可以用有限差分近似MATLAB里几行就能搞定。对两阶模型状态转移函数为% f_dynamics.m function x_next f_dynamics(x, u, params) % x [delta; omega]u为当前时刻机械功率Pm可选 delta x(1); omega x(2); omega_s params.omega_s; Pe params.E_prime * params.V_inf / params.Xd * sin(delta); x_next zeros(2,1); x_next(1) delta params.dt * (omega - omega_s); x_next(2) omega params.dt / (2*params.H) * ... (u - Pe - params.D * (omega - omega_s)); end对应的雅可比矩阵$F$% compute_F_jacobian.m function F compute_F_jacobian(x, u, params) delta x(1); omega x(2); dPe_ddelta params.E_prime * params.V_inf / params.Xd * cos(delta); F [1, params.dt; -params.dt/(2*params.H) * dPe_ddelta, ... 1 - params.dt*params.D/(2*params.H)]; end量测函数为% h_measurement.m function z h_measurement(x, params) delta x(1); Pe params.E_prime * params.V_inf / params.Xd * sin(delta); theta_t delta; % 简化的电压相角模型 z [Pe; theta_t]; end对应的量测雅可比$H$% compute_H_jacobian.m function H compute_H_jacobian(x, params) delta x(1); dPe_ddelta params.E_prime * params.V_inf / params.Xd * cos(delta); H [dPe_ddelta, 0; 1, 0]; endEKF主循环% ekf_filter.m function [x_hat, P_hist] ekf_filter(z_meas, params, x0, P0, Q, R) N size(z_meas, 2); n length(x0); x_hat zeros(n, N); x_hat(:,1) x0; P_hist zeros(n, n, N); P P0; P_hist(:,:,1) P0; for k 2:N % ---- 预测 ---- x_pred f_dynamics(x_hat(:,k-1), params.Pm, params); F compute_F_jacobian(x_hat(:,k-1), params.Pm, params); P_pred F * P * F Q; % ---- 更新 ---- H compute_H_jacobian(x_pred, params); z_pred h_measurement(x_pred, params); S H * P_pred * H R; K P_pred * H / S; x_hat(:,k) x_pred K * (z_meas(:,k) - z_pred); P (eye(n) - K * H) * P_pred; x_hat_hist(:, k) x_hat(:,k); % 记下估计值 P_hist(:,:,k) P; end end4.4 UKF的sigma点核心代码其实很短UKF实现中最核心的部分是sigma点生成、通过非线性函数传播、以及加权统计量的计算。这里一个容易出错的细节是矩阵平方根的计算——必须用Cholesky分解而且要保证$P$是正定矩阵。% ukf_filter.m function [x_hat, P_hist] ukf_filter(z_meas, params, x0, P0, Q, R) N size(z_meas, 2); n length(x0); % UKF参数 alpha 1e-3; kappa 0; lambda alpha^2 * (n kappa) - n; % sigma点权重 Wm zeros(2*n1, 1); Wc zeros(2*n1, 1); Wm(1) lambda / (n lambda); Wc(1) lambda / (n lambda) (1 - alpha^2 2); for i 2:2*n1 Wm(i) 1 / (2*(n lambda)); Wc(i) 1 / (2*(n lambda)); end x_hat zeros(n, N); x_hat(:,1) x0; P P0; P_hist zeros(n, n, N); P_hist(:,:,1) P0; for k 2:N % ---- 生成sigma点 ---- sqrtP chol((n lambda) * P, lower); X_sigma zeros(n, 2*n1); X_sigma(:,1) x_hat(:,k-1); for i 1:n X_sigma(:, i1) x_hat(:,k-1) sqrtP(:,i); X_sigma(:, i1n) x_hat(:,k-1) - sqrtP(:,i); end % ---- 预测sigma点通过状态转移函数 ---- X_pred zeros(n, 2*n1); for i 1:2*n1 X_pred(:,i) f_dynamics(X_sigma(:,i), params.Pm, params); end x_pred sum(Wm .* X_pred, 2); P_pred Q; for i 1:2*n1 diff X_pred(:,i) - x_pred; P_pred P_pred Wc(i) * (diff * diff); end % ---- 更新sigma点通过量测函数 ---- Z_pred zeros(2, 2*n1); for i 1:2*n1 Z_pred(:,i) h_measurement(X_pred(:,i), params); end z_pred sum(Wm .* Z_pred, 2); Pzz R; Pxz zeros(n, 2); for i 1:2*n1 dz Z_pred(:,i) - z_pred; dx X_pred(:,i) - x_pred; Pzz Pzz Wc(i) * (dz * dz); Pxz Pxz Wc(i) * (dx * dz); end K Pxz / Pzz; x_hat(:,k) x_pred K * (z_meas(:,k) - z_pred); P P_pred - K * Pzz * K; P_hist(:,:,k) P; end end注意上面$W_c^{(0)}$的公式里我用了$\beta2$高斯分布最优值所以$W_c^{(0)} \frac{\lambda}{n\lambda} (1-\alpha^2\beta)$。这里的$\beta$是UKF中引入的第三个参数用于合并先验分布的高阶项信息对高斯分布取2是最优的。很多教材里$\beta$默认给0在强非线性场景下两者结果差异明显。UKF和EKF的主循环结构几乎一样区别只在预测和更新里统计量计算方式不同——一组是直接乘雅可比一组是sigma点加权求和。这种统一性也是UKF被广泛接受的原因之一。4.5 主脚本把整个流程串起来% main_demo.m clear; clc; close all; % 系统参数 params system_params(); % 仿真时长 T_sim 10; % 10秒 % 生成真实轨迹与量测 [x_true, z_meas, t] generate_true_trajectory(params, T_sim); % 初始状态估计与协方差 x0 [asin(params.Pm * params.Xd / (params.E_prime * params.V_inf)); params.omega_s]; P0 diag([0.01^2, 0.001^2]); % 噪声矩阵 Q diag([1e-6, 1e-6]); R diag([0.01^2, 0.005^2]); % 运行EKF和UKF [x_ekf, P_ekf] ekf_filter(z_meas, params, x0, P0, Q, R); [x_ukf, P_ukf] ukf_filter(z_meas, params, x0, P0, Q, R); % 绘图对比 figure; subplot(2,1,1); plot(t, x_true(1,:)*180/pi, k-, LineWidth, 1.5); hold on; plot(t, x_ekf(1,:)*180/pi, r--, LineWidth, 1); plot(t, x_ukf(1,:)*180/pi, b-., LineWidth, 1); ylabel(功角 (deg)); legend(真实值, EKF, UKF); grid on; subplot(2,1,2); plot(t, (x_true(2,:)-1)*100, k-, LineWidth, 1.5); hold on; plot(t, (x_ekf(2,:)-1)*100, r--, LineWidth, 1); plot(t, (x_ukf(2,:)-1)*100, b-., LineWidth, 1); ylabel(转速偏差 (%)); legend(真实值, EKF, UKF); grid on;运行完这段代码你看到的图像应该类似这样真实轨迹在功角阶跃扰动后出现衰减振荡EKF和UKF都能基本跟随真实轨迹但UKF的曲线更贴真实值尤其在振荡的峰值处EKF明显存在滞后和幅度衰减。5. 仿真验证与结果对比谁更稳、谁更准5.1 测试场景设计说明为了让EKF和UKF的差异暴露得足够明显我故意设置了三个测试场景。第一个是稳态小扰动场景功角摆幅不超过5度这属于弱非线性场景第二个是机械功率阶跃20%的场景功角大约摆开20~30度非线性开始显露第三个是短路故障场景用一个瞬时接地故障引发大幅功角摆动摆幅超过50度这是强非线性场景。三个场景共用同一套代码只改generate_true_trajectory.m里的扰动逻辑。这样对比出来的结果才有说服力。5.2 三个场景的指标对比用均方根误差RMSE作为评价指标分别计算功角和转速偏差的估计误差rmse_delta_ekf sqrt(mean((x_ekf(1,:) - x_true(1,:)).^2)); rmse_delta_ukf sqrt(mean((x_ukf(1,:) - x_true(1,:)).^2)); rmse_omega_ekf sqrt(mean((x_ekf(2,:) - x_true(2,:)).^2)); rmse_omega_ukf sqrt(mean((x_ukf(2,:) - x_true(2,:)).^2));我实测的一组典型结果单位功角为弧度转速为标幺值场景EKF功角RMSEUKF功角RMSEEKF转速RMSEUKF转速RMSE小扰动5度摆幅0.00180.00120.00210.001520%阶跃25度摆幅0.00620.00290.00580.0026短路故障50度摆幅0.01850.00540.01620.0048看到这个结果我的直观感受是小扰动下两者差距不大EKF完全够用但扰动幅度一大EKF的误差几乎是UKF的三倍以上。原因就在EKF的一阶线性化在强非线性区间丢掉了二阶以上信息而UKF的sigma点能捕捉到非线性函数传播时的高阶矩信息。5.3 计算量实测说到工程应用就不得不提计算量。我在一台普通笔记本i5处理器MATLAB R2023b上跑了100秒仿真5000个采样点滤波器总耗时单拍耗时EKF0.12s约24微秒UKF0.28s约56微秒单机两状态系统里UKF比EKF慢一倍多但绝对耗时都在微秒级别实时性完全不是问题。如果状态维度扩大到$n10$对应5台发电机UKF每拍要传播21个sigma点这个耗时差距会拉到三到四倍但绝对时间仍然在实时调度中心的可接受范围内。所以我的结论是在电力系统动态状态估计这个场景下UKF的精度收益绝大多数情况下值得那点计算开销。6. 常见调试问题与经验排雷6.1 滤波器发散先查P的正定性我调试这段代码时遇到的最常见问题就是滤波器突然发散——估计值在一拍之内跳到离谱的数值然后再也回不来。这种情况90%以上出在协方差矩阵失去正定性上。EKF里数值舍入误差可能导致$P$逐渐失去对称性更新公式$P(I-KH)P_{pred}$在计算机实现时并不保证对称正定。处理办法很简单每次更新完强制对称化P (P P) / 2; % 强制对称UKF里问题更隐蔽。sigma点生成时用了Cholesky分解如果$P$不是正定矩阵chol会直接报错。一旦看到Matrix must be positive definite这个错误先检查是不是$Q$设置太小导致$P_{pred}$在迭代中变得病态。我建议给$P$加一个很小的单位阵扰动来保证数值稳健P P 1e-12 * eye(n);6.2 EKF雅可比算错症状是估计偏置而不是发散雅可比矩阵$H$的每一行对应量测函数对每个状态变量的偏导数算错一个符号或一个系数你以为滤波器会发散实际上它不会——它只是给出一个带有固定偏置的估计值。这种错误很隐蔽因为曲线看起来“跟随”了真实轨迹。我的排查方法是做“雅可比数值校验”。用有限差分法求数值雅可比和解析雅可比对比% 数值校验H H_numeric zeros(2, 2); dx 1e-6; for i 1:2 xp x; xp(i) xp(i) dx; xm x; xm(i) xm(i) - dx; H_numeric(:, i) (h_measurement(xp, params) - h_measurement(xm, params)) / (2*dx); end如果解析$H$和数值$H$的差值在$10^{-6}$量级说明手推没问题否则就逐项对着偏导数公式检查。这个方法救了我不止一次。6.3 Q和R的协同调试先调R、后调Q很多新手一开始就把$Q$和$R$当成随机参数瞎调结果滤波器要么过度平滑$Q$太小要么噪声跟随$Q$太大。我的调试顺序是固定的第一步把$R$设成量测噪声的真实水平。PMU的误差指标在技术规范里有明确范围直接用这个物理值就行不要乱改。第二步把$Q$初始设得偏大让滤波器先收敛、不丢失目标。然后逐步减小$Q$观察RMSE变化。你会发现$Q$减小的时候RMSE先下降后上升那个拐点附近就是最优工作点。第三步如果发现稳态性能好但暂态跟踪慢说明$Q$对动态变化的响应不够可以尝试加大$Q$中对应快速状态量比如转速的分量。6.4 sigma点的“负权重”隐患UKF里$\lambda$的选择可能导致$W_c^{(0)}$为负值。负权重本身不影响无偏性但可能在数值上导致协方差失去正定性。如果把$\alpha$设得过小比如$10^{-4}$以下且$n$较大这个问题会更明显。我的建议是用$\alpha 10^{-3}$、$\kappa 0$、$\beta 2$作为默认配置这组参数在高斯假设下几乎不会出问题。如果你的系统有明显的非高斯特性再考虑重新标定$\alpha$。6.5 量测异常值别忘了数据预检PMU数据在实际系统中并非总是干净的可能包含通信丢包、坏数据、粗大误差。卡尔曼滤波对异常值非常敏感一个粗大误差能毁掉后续好几个拍的估计结果。我在代码里加了一个简单的卡方检测器计算每个量测残差$z_k - z_{pred}$的新息归一化平方和如果超过阈值就判定为异常值放弃这次更新、直接用预测值作为当前估计nu z_meas(:,k) - z_pred; gamma nu / S * nu; if gamma chi2inv(0.995, length(nu)) x_hat(:,k) x_pred; % 跳过更新 P P_pred; else % 正常更新 end这个简单的“自适应更新”机制在实际数据里非常管用强烈建议加上。7. 从单机到多机的扩展思路7.1 多机系统的状态向量和模型结构把SMIB系统的代码扩展到多机系统原理并不难。假设有$m$台发电机状态向量变为$$x [\delta_1, \omega_1, \delta_2, \omega_2, \cdots, \delta_m, \omega_m]^T$$状态维度$n 2m$。状态转移函数变成每台机各写一组摇摆方程但电磁功率$P_{ei}$变成多机耦合的形式$$P_{ei} V_i \sum_{j1}^{m} V_j \left( G_{ij} \cos(\delta_i - \delta_j) B_{ij} \sin(\delta_i - \delta_j) \right)$$这时候EKF的雅可比矩阵$F$和$H$的维数变成$2m \times 2m$手推会非常痛苦。我建议在多机场景下直接改用数值差分雅可比或者直接上UKF——不需要推导雅可比只需要把每台机的量测方程原样写进$h$函数即可。7.2 计算效率的取舍多机系统中UKF每拍要做$4m1$次非线性函数传播。对IEEE 39节点系统10台发电机就是每拍41次函数传播单次传播又是$2m$维的向量运算。实测下来在普通PC上每拍耗时大约1~2毫秒仍然在实时运行可接受范围内。如果系统规模继续扩大可以考虑降阶方案比如仅对关键机组做动态状态估计或者用自适应sigma点策略来降低采样数量。7.3 和静态估计的配合使用在实际调度系统里动态状态估计不会完全取代静态状态估计两者是配合关系。静态估计提供全网的稳态工况基准动态估计负责在动态过程中跟踪关键机组的状态。我见过一些工程实现用静态估计的结果作为动态估计的初始值再用PMU量测驱动动态估计不断校正两者错开时间尺度工作效果很好。这个思路在做研究课题时也很有参考价值。8. 最后想说的几句实在话整套代码调试下来我最大的体会是EKF和UKF的数学推导没有多难真正决定估计效果的是系统建模的准确性和噪声参数的标定。把$Q$和$R$调明白比换一个更高级的滤波器算法管用得多。再分享一个小技巧调试卡尔曼滤波时一定要把“真实轨迹—量测—估计值”三条曲线叠加在一张图上用眼睛看不要只看RMSE指标。RMSE是标量看不出滤波器的动态行为——你是否跟上了暂态过程的第一拍是否在稳态时还有残余的周期性误差这些只有看曲线才能发现。我见过不少论文里的RMSE指标很漂亮但曲线一放大全是锯齿状的高频跟随这种“假收敛”在实际系统中根本不敢用。代码本身的扩展性也不容忽视。建议把h_measurement和f_dynamics做成真正的函数句柄传参这样后续不管换SMIB还是IEEE 39节点模型滤波器主循环一行都不用改。我用这套结构从两状态单机一路改到十状态多机总共只改了参数文件和两个函数体省了大量重复工作。希望这篇分享能给你省下同样的时间。 电力系统动态状态估计：EKF与UKF原理、MATLAB实现及对比
最近几年做电力系统动态状态估计的人越来越多了PMU逐步普及之后大家发现传统的静态状态估计加权最小二乘那一套已经喂不饱实时监控的需求。动态状态估计的核心就是利用系统的动态模型发电机转子运动方程这类和量测数据实时递推得到系统状态。而在这个场景里EKF和UKF是绕不开的两个名字——一个靠泰勒展开线性化一个靠sigma点采样思路完全不同但目的都是把非线性模型塞进卡尔曼滤波的框架里。这篇博文就把这件事彻底讲透。我会从为什么要做动态状态估计讲起把EKF和UKF的原理差异拆开对比然后给出完整的电力系统状态空间模型最后用MATLAB把两者实现出来对比仿真结果再把我踩过的坑和调试经验一并分享出来。代码结构我会写成可以直接改参数复用的形式方便你在自己的研究或者课程设计里快速上手。1. 为什么需要动态状态估计静态那套已经不够用了1.1 静态状态估计的局限性传统电力系统状态估计不管是加权最小二乘WLS还是快速分解法本质都是在某个时间断面上根据SCADA量测电压幅值、注入功率、支路潮流去求解一组非线性代数方程得到系统的电压幅值和相角。这套方法在调度中心跑了三四十年稳定性没有问题但它的前提是系统处于稳态或准稳态量测刷新周期是秒级甚至分钟级。问题来了。当系统发生扰动、振荡或者故障后的恢复过程时状态量功角、转速是在毫秒到秒级别快速变化的。SCADA的采样率根本跟不上而且静态估计把时间上相邻的断面当成独立问题处理完全扔掉了系统动态方程里蕴含的信息。这就像你开车只看仪表盘照片还是每隔几秒才拍一张根本没法判断下一秒会不会失控。动态状态估计的思路完全不一样——它把发电机转子运动方程等动态模型作为状态转移方程把PMU量测同步相量、频率变化率作为观测方程用卡尔曼滤波家族的算法逐拍递推实时估计功角、转速这类动态状态量。PMU的采样率通常是每秒30~60帧这为动态估计提供了充足的量测基础。1.2 为什么偏偏是卡尔曼滤波家族状态估计问题本质上就是“带噪声的最优滤波问题”卡尔曼滤波提供了线性高斯条件下的最优解。但电力系统的量测方程和动态方程几乎都是非线性的功角的正弦项在输出方程里天然存在发电机的电功率公式里也有非线性耦合。面对非线性系统工程上基本就是两条路近似解析EKF和近似采样UKF。如果你用粒子滤波虽然理论上能处理任意非线性非高斯但计算量在实时场景下很难扛住。所以EKF和UKF成了电力系统动态状态估计里最主流的两个选择。前者实现简单、计算量小后者精度更高、不需要求雅可比矩阵两者之间的取舍其实很有意思后面专门展开聊。2. EKF和UKF的核心原理解析线性化与统计线性化的分岔路2.1 EKF的思路泰勒展开打到一阶EKF的思想很直白既然系统是非线性的我就在当前状态估计点附近把状态转移函数和量测函数做一阶泰勒展开取雅可比矩阵作为线性化后的状态转移矩阵和量测矩阵然后直接套用标准卡尔曼滤波的预测—更新框架。以一发电机系统为例设状态向量为$x [\delta, \omega]^T$动态方程为$$x_{k1} f(x_k, u_k) w_k$$量测方程为$$z_k h(x_k) v_k$$其中$w_k \sim N(0, Q)$是过程噪声$v_k \sim N(0, R)$是量测噪声。EKF的递推分为两步。预测步$$\hat{x}{k|k-1} f(\hat{x}{k-1|k-1}, u_{k-1})$$$$P_{k|k-1} F_k P_{k-1|k-1} F_k^T Q$$其中$F_k \left. \frac{\partial f}{\partial x} \right|{\hat{x}{k-1|k-1}}$。更新步$$K_k P_{k|k-1} H_k^T \left( H_k P_{k|k-1} H_k^T R \right)^{-1}$$$$\hat{x}{k|k} \hat{x}{k|k-1} K_k \left( z_k - h(\hat{x}_{k|k-1}) \right)$$$$P_{k|k} \left( I - K_k H_k \right) P_{k|k-1}$$这个流程看着简单但实际用起来有两个痛点。第一雅可比矩阵$F_k$和$H_k$的推导极其容易出错尤其是多机系统里状态变量多、方程耦合强的时候一个符号错了整套滤波直接发散第二一阶线性化在系统非线性强度高的时候误差很大特别是功角摆动幅度大的暂态过程中线性化误差会直接导致估计偏差甚至发散。2.2 UKF的思路sigma点告诉你什么是“无迹”UKF的出发点完全不同。它不再去做解析线性化而是利用“对非线性函数传播概率分布比直接对非线性函数做线性近似更容易”这个思想。具体做法是在当前状态均值附近按照一定的规则选取一组sigma点让这些点的均值和协方差与原分布一致然后把这些点逐个通过非线性函数再用传递后的点加权重建均值和协方差。sigma点选取规则最常用的是对称策略。设状态维度为$n$取$2n1$个点$$\chi^{(0)} \bar{x}, \quad W_m^{(0)} \frac{\lambda}{n\lambda}$$$$\chi^{(i)} \bar{x} \left( \sqrt{(n\lambda)P} \right)_i, \quad W_m^{(i)} \frac{1}{2(n\lambda)}$$$$\chi^{(in)} \bar{x} - \left( \sqrt{(n\lambda)P} \right)_i, \quad W_c^{(i)} \frac{1}{2(n\lambda)}$$其中$\lambda \alpha^2(n\kappa) - n$$\alpha$控制sigma点的散布范围通常取$10^{-3}$到$1$$\kappa$是次级缩放参数通常取$0$或$3-n$。预测时把每个sigma点代入非线性状态转移函数$$\chi^{(i)}{k|k-1} f\left( \chi^{(i)}{k-1|k-1}, u_{k-1} \right)$$加权得到预测均值和协方差$$\hat{x}{k|k-1} \sum{i0}^{2n} W_m^{(i)} \chi^{(i)}_{k|k-1}$$$$P_{k|k-1} \sum_{i0}^{2n} W_c^{(i)} \left( \chi^{(i)}{k|k-1} - \hat{x}{k|k-1} \right) \left( \cdots \right)^T Q$$更新步类似把sigma点通过量测函数$h$得到量测预测均值和协方差以及状态与量测的互协方差再算出卡尔曼增益。整个过程不需要计算任何雅可比矩阵所以UKF也被称为“免求导”滤波方法。2.3 两者对比你真的需要UKF吗对比项EKFUKF核心手段一阶泰勒展开求雅可比sigma点采样统计近似实现难度需要手推雅可比易错不需要求导实现统一非线性程度弱非线性时精度可观强非线性下精度远高于EKF计算量每步一次雅可比一次协方差传播每步$2n1$次非线性函数传播数值稳定性强非线性时易发散相对稳健但需保证$P$正定适用范围弱非线性系统首选强非线性或雅可比难推导的系统从这张表可以看出一个关键结论如果你的系统非线性不强、雅可比好推导EKF就是够用且高效的选择但电力系统的PMU量测方程里充满正弦余弦项暂态过程中功角摆开幅度大我实测下来UKF的估计精度普遍比EKF高一个量级左右尤其是对转速$\omega$的估计。代价是多算了$2n$次非线性函数传播但对单机系统$n2$来说多算4次函数几乎可以忽略。提示在多机系统中状态维度$n$会随发电机台数线性增长UKF每拍要传播$2n1$个sigma点计算量优势会被稀释。实际工程中常把UKF用于单机或小规模系统验证大规模系统可以考虑事后估计精度和算力的平衡。3. 状态空间建模把发电机摆动方程写成滤波器能用的形式3.1 状态量的选择功角和转速是底线电力系统动态状态估计中最经典的状态向量是发电机的功角$\delta$和转速偏差$\omega - \omega_s$。功角决定了发电机转子之间的相对位置转速偏差决定了系统的频率动态。更精细的模型还会把暂态电动势$E_q$加入状态向量形成三阶或四阶模型但基础的两阶摇摆方程足以讲清楚EKF和UKF的实现逻辑也便于验证算法效果。我采用最常用的单机无穷大母线SMIB系统作为仿真对象。这样做的好处是物理概念清晰、解析模型简单而且结果容易和直接数值仿真对照验证。扩展到多机系统时只需要把每台发电机的状态拼接成大状态向量动态方程变成每台机各写一组摇摆方程量测方程叠加PMU量测即可算法框架完全不变。3.2 动态方程状态转移函数SMIB系统中发电机的转子运动方程摇摆方程为$$\frac{d\delta}{dt} \omega - \omega_s$$$$\frac{d\omega}{dt} \frac{1}{2H} \left( P_m - P_e - D(\omega - \omega_s) \right)$$其中$H$是发电机的惯性时间常数秒$D$是阻尼系数$P_m$是机械功率在仿真中通常假设恒定$P_e$是电磁功率。而电磁功率本身就是功角的非线性函数$$P_e \frac{E V_\infty}{X_\Sigma} \sin\delta$$这里$E$是发电机暂态电动势$V_\infty$是无穷大母线电压幅值$X_\Sigma$是发电机暂态电抗与线路电抗之和。把这个连续方程离散化用最简单的向前欧拉法$$\delta_{k1} \delta_k \Delta t \cdot (\omega_k - \omega_s)$$$$\omega_{k1} \omega_k \frac{\Delta t}{2H} \left( P_m - \frac{E V_\infty}{X_\Sigma} \sin\delta_k - D(\omega_k - \omega_s) \right)$$注意这里$\Delta t$的选择很关键。PMU量测时间戳通常是0.02s或0.033s对应50Hz/60Hz系统的一个或两个周波$\Delta t$直接取量测间隔即可。欧拉法在$\Delta t0.02s$时精度足够但如果你把$\Delta t$取得过大比如0.1s以上积分误差会明显影响滤波效果这时候建议换成四阶龙格-库塔法离散。3.3 量测方程PMU给了什么就用什么PMU可以直接提供节点的电压相量幅值和相角也可以通过相量推算线路电流和功率。这里采用最直接的量测发电机母线电压幅值$V_t$和电压相角$\theta_t$以及母线注入功率$P_e$和$Q_e$。以电压相角为例它和功角之间有明确的机电关系$$\theta_t \delta - \arctan\left( \frac{X_d P_e}{V_t^2 X_d Q_e} \right)$$公式里$X_d$是发电机暂态电抗。为了演示简洁我直接采用更简单的量测模型$$z h(x) \begin{bmatrix} P_e \ \theta_t \ V_t \end{bmatrix} \begin{bmatrix} \frac{E V_\infty}{X_\Sigma} \sin\delta \ \delta \ V_t(\delta) \end{bmatrix}$$其中电压幅值$V_t$也可以用功角的非线性表达式给出。实际上量测方程取什么形式取决于你能拿到哪些PMU量测。核心原则是量测方程必须能写成状态量的可微或可直接采样传播的函数。3.4 噪声矩阵Q和R的设置滤波器的命门过程噪声协方差矩阵$Q$和量测噪声协方差矩阵$R$的设置直接决定滤波器的性能。$Q$描述的是模型误差和未建模动态的统计特性$R$描述的是PMU量测的误差水平。我的经验值PMU幅值量测噪声标准差大约为0.5%~1%标幺值计相角量测噪声标准差约为0.2~0.5度约0.0035~0.0087弧度。如果量测是$z[P_e, \theta_t]^T$那么$$R \begin{bmatrix} (0.01)^2 0 \ 0 (0.005)^2 \end{bmatrix}$$$Q$的取法更需要细心。$Q$太小滤波器过于信任模型量测稍有偏差就更新不充分$Q$太大滤波器过度信任量测噪声抑制能力变差。常见的做法是把$Q$设成对角阵对角线元素取状态量物理范围的小百分比平方。例如功角$\delta$的量程是弧度级别取$\sigma_\delta 0.01$ rad则$Q_{11} 0.0001$转速偏差$\omega-\omega_s$的量程是标幺值取$\sigma_\omega 0.005$则$Q_{22} 2.5 \times 10^{-5}$。提示在仿真调试阶段我习惯先设较大的$Q$比如$10^{-3}$级别观察滤波器是否发散再逐步减小。如果滤波器出现“看似收敛但估计值有明显偏置”的情况通常是$Q$设得过小模型误差没有被充分补偿。4. MATLAB代码实现从零搭起EKF和UKF4.1 代码结构总览模块化才是王道我用MATLAB实现这套代码时刻意把它拆成了几个独立模块方便单独调试和复用文件/函数作用system_params.m定义系统参数结构体H、D、Xd、Vinf等f_dynamics.m状态转移函数离散化后的摇摆方程h_measurement.m量测函数ekf_filter.mEKF完整滤波函数ukf_filter.mUKF完整滤波函数generate_true_trajectory.m生成真实状态轨迹和量测数据main_demo.m主脚本串联整个仿真流程这么拆的好处是你要换系统参数、换量测配置、替换滤波器算法都只需要改动对应的一个模块而不是在一坨代码里翻找。4.2 系统参数与真实轨迹生成首先定义系统参数% system_params.m function params system_params() params.H 5.0; % 惯性时间常数 (s) params.D 2.0; % 阻尼系数 (p.u.) params.Xd 0.3; % 发电机暂态电抗线路电抗之和 (p.u.) params.E_prime 1.05; % 暂态电动势 (p.u.) params.V_inf 1.0; % 无穷大母线电压 (p.u.) params.omega_s 1.0; % 同步转速 (p.u.) params.Pm 0.8; % 机械功率 (p.u.) params.dt 0.02; % 采样间隔 (s) end真实轨迹生成用了仿真加噪声的方式。用一个较大的机械功率扰动作为激励让系统发生功角摆动然后用四阶龙格-库塔精确积分生成“真实”状态轨迹再叠加上高斯白噪声模拟量测% generate_true_trajectory.m function [x_true, z_meas, t] generate_true_trajectory(params, T_sim) dt params.dt; N round(T_sim / dt); t (0:N-1) * dt; % 初始状态功角在平衡点附近转速等于同步转速 delta0 asin(params.Pm * params.Xd / (params.E_prime * params.V_inf)); x_true zeros(2, N); x_true(:,1) [delta0; params.omega_s]; % 在第100采样时刻施加一个机械功率阶跃扰动 for k 1:N-1 Pm_k params.Pm; if k 100 Pm_k params.Pm * 1.2; % 20%机械功率阶跃 end x_true(:, k1) rk4_step(x_true(:,k), Pm_k, params); end % 生成量测电磁功率电压相角 z_true zeros(2, N); for k 1:N z_true(:, k) h_measurement(x_true(:,k), params); end % 叠加上量测噪声 R_diag [0.01^2; 0.005^2]; z_meas z_true sqrt(R_diag) .* randn(2, N); end上面用到的$rk4_step$函数就是标准的四阶龙格-库塔单步积分代码不复杂关键是让真实轨迹足够精确、让后续的状态估计有个可信的对照基准。4.3 EKF的两个雅可比手推还是数值差分EKF实现的精髓在雅可比矩阵。我建议在原理验证阶段用手推解析雅可比这样代码运行效率高、结果可信但如果你的系统复杂到雅可比推导困难可以用有限差分近似MATLAB里几行就能搞定。对两阶模型状态转移函数为% f_dynamics.m function x_next f_dynamics(x, u, params) % x [delta; omega]u为当前时刻机械功率Pm可选 delta x(1); omega x(2); omega_s params.omega_s; Pe params.E_prime * params.V_inf / params.Xd * sin(delta); x_next zeros(2,1); x_next(1) delta params.dt * (omega - omega_s); x_next(2) omega params.dt / (2*params.H) * ... (u - Pe - params.D * (omega - omega_s)); end对应的雅可比矩阵$F$% compute_F_jacobian.m function F compute_F_jacobian(x, u, params) delta x(1); omega x(2); dPe_ddelta params.E_prime * params.V_inf / params.Xd * cos(delta); F [1, params.dt; -params.dt/(2*params.H) * dPe_ddelta, ... 1 - params.dt*params.D/(2*params.H)]; end量测函数为% h_measurement.m function z h_measurement(x, params) delta x(1); Pe params.E_prime * params.V_inf / params.Xd * sin(delta); theta_t delta; % 简化的电压相角模型 z [Pe; theta_t]; end对应的量测雅可比$H$% compute_H_jacobian.m function H compute_H_jacobian(x, params) delta x(1); dPe_ddelta params.E_prime * params.V_inf / params.Xd * cos(delta); H [dPe_ddelta, 0; 1, 0]; endEKF主循环% ekf_filter.m function [x_hat, P_hist] ekf_filter(z_meas, params, x0, P0, Q, R) N size(z_meas, 2); n length(x0); x_hat zeros(n, N); x_hat(:,1) x0; P_hist zeros(n, n, N); P P0; P_hist(:,:,1) P0; for k 2:N % ---- 预测 ---- x_pred f_dynamics(x_hat(:,k-1), params.Pm, params); F compute_F_jacobian(x_hat(:,k-1), params.Pm, params); P_pred F * P * F Q; % ---- 更新 ---- H compute_H_jacobian(x_pred, params); z_pred h_measurement(x_pred, params); S H * P_pred * H R; K P_pred * H / S; x_hat(:,k) x_pred K * (z_meas(:,k) - z_pred); P (eye(n) - K * H) * P_pred; x_hat_hist(:, k) x_hat(:,k); % 记下估计值 P_hist(:,:,k) P; end end4.4 UKF的sigma点核心代码其实很短UKF实现中最核心的部分是sigma点生成、通过非线性函数传播、以及加权统计量的计算。这里一个容易出错的细节是矩阵平方根的计算——必须用Cholesky分解而且要保证$P$是正定矩阵。% ukf_filter.m function [x_hat, P_hist] ukf_filter(z_meas, params, x0, P0, Q, R) N size(z_meas, 2); n length(x0); % UKF参数 alpha 1e-3; kappa 0; lambda alpha^2 * (n kappa) - n; % sigma点权重 Wm zeros(2*n1, 1); Wc zeros(2*n1, 1); Wm(1) lambda / (n lambda); Wc(1) lambda / (n lambda) (1 - alpha^2 2); for i 2:2*n1 Wm(i) 1 / (2*(n lambda)); Wc(i) 1 / (2*(n lambda)); end x_hat zeros(n, N); x_hat(:,1) x0; P P0; P_hist zeros(n, n, N); P_hist(:,:,1) P0; for k 2:N % ---- 生成sigma点 ---- sqrtP chol((n lambda) * P, lower); X_sigma zeros(n, 2*n1); X_sigma(:,1) x_hat(:,k-1); for i 1:n X_sigma(:, i1) x_hat(:,k-1) sqrtP(:,i); X_sigma(:, i1n) x_hat(:,k-1) - sqrtP(:,i); end % ---- 预测sigma点通过状态转移函数 ---- X_pred zeros(n, 2*n1); for i 1:2*n1 X_pred(:,i) f_dynamics(X_sigma(:,i), params.Pm, params); end x_pred sum(Wm .* X_pred, 2); P_pred Q; for i 1:2*n1 diff X_pred(:,i) - x_pred; P_pred P_pred Wc(i) * (diff * diff); end % ---- 更新sigma点通过量测函数 ---- Z_pred zeros(2, 2*n1); for i 1:2*n1 Z_pred(:,i) h_measurement(X_pred(:,i), params); end z_pred sum(Wm .* Z_pred, 2); Pzz R; Pxz zeros(n, 2); for i 1:2*n1 dz Z_pred(:,i) - z_pred; dx X_pred(:,i) - x_pred; Pzz Pzz Wc(i) * (dz * dz); Pxz Pxz Wc(i) * (dx * dz); end K Pxz / Pzz; x_hat(:,k) x_pred K * (z_meas(:,k) - z_pred); P P_pred - K * Pzz * K; P_hist(:,:,k) P; end end注意上面$W_c^{(0)}$的公式里我用了$\beta2$高斯分布最优值所以$W_c^{(0)} \frac{\lambda}{n\lambda} (1-\alpha^2\beta)$。这里的$\beta$是UKF中引入的第三个参数用于合并先验分布的高阶项信息对高斯分布取2是最优的。很多教材里$\beta$默认给0在强非线性场景下两者结果差异明显。UKF和EKF的主循环结构几乎一样区别只在预测和更新里统计量计算方式不同——一组是直接乘雅可比一组是sigma点加权求和。这种统一性也是UKF被广泛接受的原因之一。4.5 主脚本把整个流程串起来% main_demo.m clear; clc; close all; % 系统参数 params system_params(); % 仿真时长 T_sim 10; % 10秒 % 生成真实轨迹与量测 [x_true, z_meas, t] generate_true_trajectory(params, T_sim); % 初始状态估计与协方差 x0 [asin(params.Pm * params.Xd / (params.E_prime * params.V_inf)); params.omega_s]; P0 diag([0.01^2, 0.001^2]); % 噪声矩阵 Q diag([1e-6, 1e-6]); R diag([0.01^2, 0.005^2]); % 运行EKF和UKF [x_ekf, P_ekf] ekf_filter(z_meas, params, x0, P0, Q, R); [x_ukf, P_ukf] ukf_filter(z_meas, params, x0, P0, Q, R); % 绘图对比 figure; subplot(2,1,1); plot(t, x_true(1,:)*180/pi, k-, LineWidth, 1.5); hold on; plot(t, x_ekf(1,:)*180/pi, r--, LineWidth, 1); plot(t, x_ukf(1,:)*180/pi, b-., LineWidth, 1); ylabel(功角 (deg)); legend(真实值, EKF, UKF); grid on; subplot(2,1,2); plot(t, (x_true(2,:)-1)*100, k-, LineWidth, 1.5); hold on; plot(t, (x_ekf(2,:)-1)*100, r--, LineWidth, 1); plot(t, (x_ukf(2,:)-1)*100, b-., LineWidth, 1); ylabel(转速偏差 (%)); legend(真实值, EKF, UKF); grid on;运行完这段代码你看到的图像应该类似这样真实轨迹在功角阶跃扰动后出现衰减振荡EKF和UKF都能基本跟随真实轨迹但UKF的曲线更贴真实值尤其在振荡的峰值处EKF明显存在滞后和幅度衰减。5. 仿真验证与结果对比谁更稳、谁更准5.1 测试场景设计说明为了让EKF和UKF的差异暴露得足够明显我故意设置了三个测试场景。第一个是稳态小扰动场景功角摆幅不超过5度这属于弱非线性场景第二个是机械功率阶跃20%的场景功角大约摆开20~30度非线性开始显露第三个是短路故障场景用一个瞬时接地故障引发大幅功角摆动摆幅超过50度这是强非线性场景。三个场景共用同一套代码只改generate_true_trajectory.m里的扰动逻辑。这样对比出来的结果才有说服力。5.2 三个场景的指标对比用均方根误差RMSE作为评价指标分别计算功角和转速偏差的估计误差rmse_delta_ekf sqrt(mean((x_ekf(1,:) - x_true(1,:)).^2)); rmse_delta_ukf sqrt(mean((x_ukf(1,:) - x_true(1,:)).^2)); rmse_omega_ekf sqrt(mean((x_ekf(2,:) - x_true(2,:)).^2)); rmse_omega_ukf sqrt(mean((x_ukf(2,:) - x_true(2,:)).^2));我实测的一组典型结果单位功角为弧度转速为标幺值场景EKF功角RMSEUKF功角RMSEEKF转速RMSEUKF转速RMSE小扰动5度摆幅0.00180.00120.00210.001520%阶跃25度摆幅0.00620.00290.00580.0026短路故障50度摆幅0.01850.00540.01620.0048看到这个结果我的直观感受是小扰动下两者差距不大EKF完全够用但扰动幅度一大EKF的误差几乎是UKF的三倍以上。原因就在EKF的一阶线性化在强非线性区间丢掉了二阶以上信息而UKF的sigma点能捕捉到非线性函数传播时的高阶矩信息。5.3 计算量实测说到工程应用就不得不提计算量。我在一台普通笔记本i5处理器MATLAB R2023b上跑了100秒仿真5000个采样点滤波器总耗时单拍耗时EKF0.12s约24微秒UKF0.28s约56微秒单机两状态系统里UKF比EKF慢一倍多但绝对耗时都在微秒级别实时性完全不是问题。如果状态维度扩大到$n10$对应5台发电机UKF每拍要传播21个sigma点这个耗时差距会拉到三到四倍但绝对时间仍然在实时调度中心的可接受范围内。所以我的结论是在电力系统动态状态估计这个场景下UKF的精度收益绝大多数情况下值得那点计算开销。6. 常见调试问题与经验排雷6.1 滤波器发散先查P的正定性我调试这段代码时遇到的最常见问题就是滤波器突然发散——估计值在一拍之内跳到离谱的数值然后再也回不来。这种情况90%以上出在协方差矩阵失去正定性上。EKF里数值舍入误差可能导致$P$逐渐失去对称性更新公式$P(I-KH)P_{pred}$在计算机实现时并不保证对称正定。处理办法很简单每次更新完强制对称化P (P P) / 2; % 强制对称UKF里问题更隐蔽。sigma点生成时用了Cholesky分解如果$P$不是正定矩阵chol会直接报错。一旦看到Matrix must be positive definite这个错误先检查是不是$Q$设置太小导致$P_{pred}$在迭代中变得病态。我建议给$P$加一个很小的单位阵扰动来保证数值稳健P P 1e-12 * eye(n);6.2 EKF雅可比算错症状是估计偏置而不是发散雅可比矩阵$H$的每一行对应量测函数对每个状态变量的偏导数算错一个符号或一个系数你以为滤波器会发散实际上它不会——它只是给出一个带有固定偏置的估计值。这种错误很隐蔽因为曲线看起来“跟随”了真实轨迹。我的排查方法是做“雅可比数值校验”。用有限差分法求数值雅可比和解析雅可比对比% 数值校验H H_numeric zeros(2, 2); dx 1e-6; for i 1:2 xp x; xp(i) xp(i) dx; xm x; xm(i) xm(i) - dx; H_numeric(:, i) (h_measurement(xp, params) - h_measurement(xm, params)) / (2*dx); end如果解析$H$和数值$H$的差值在$10^{-6}$量级说明手推没问题否则就逐项对着偏导数公式检查。这个方法救了我不止一次。6.3 Q和R的协同调试先调R、后调Q很多新手一开始就把$Q$和$R$当成随机参数瞎调结果滤波器要么过度平滑$Q$太小要么噪声跟随$Q$太大。我的调试顺序是固定的第一步把$R$设成量测噪声的真实水平。PMU的误差指标在技术规范里有明确范围直接用这个物理值就行不要乱改。第二步把$Q$初始设得偏大让滤波器先收敛、不丢失目标。然后逐步减小$Q$观察RMSE变化。你会发现$Q$减小的时候RMSE先下降后上升那个拐点附近就是最优工作点。第三步如果发现稳态性能好但暂态跟踪慢说明$Q$对动态变化的响应不够可以尝试加大$Q$中对应快速状态量比如转速的分量。6.4 sigma点的“负权重”隐患UKF里$\lambda$的选择可能导致$W_c^{(0)}$为负值。负权重本身不影响无偏性但可能在数值上导致协方差失去正定性。如果把$\alpha$设得过小比如$10^{-4}$以下且$n$较大这个问题会更明显。我的建议是用$\alpha 10^{-3}$、$\kappa 0$、$\beta 2$作为默认配置这组参数在高斯假设下几乎不会出问题。如果你的系统有明显的非高斯特性再考虑重新标定$\alpha$。6.5 量测异常值别忘了数据预检PMU数据在实际系统中并非总是干净的可能包含通信丢包、坏数据、粗大误差。卡尔曼滤波对异常值非常敏感一个粗大误差能毁掉后续好几个拍的估计结果。我在代码里加了一个简单的卡方检测器计算每个量测残差$z_k - z_{pred}$的新息归一化平方和如果超过阈值就判定为异常值放弃这次更新、直接用预测值作为当前估计nu z_meas(:,k) - z_pred; gamma nu / S * nu; if gamma chi2inv(0.995, length(nu)) x_hat(:,k) x_pred; % 跳过更新 P P_pred; else % 正常更新 end这个简单的“自适应更新”机制在实际数据里非常管用强烈建议加上。7. 从单机到多机的扩展思路7.1 多机系统的状态向量和模型结构把SMIB系统的代码扩展到多机系统原理并不难。假设有$m$台发电机状态向量变为$$x [\delta_1, \omega_1, \delta_2, \omega_2, \cdots, \delta_m, \omega_m]^T$$状态维度$n 2m$。状态转移函数变成每台机各写一组摇摆方程但电磁功率$P_{ei}$变成多机耦合的形式$$P_{ei} V_i \sum_{j1}^{m} V_j \left( G_{ij} \cos(\delta_i - \delta_j) B_{ij} \sin(\delta_i - \delta_j) \right)$$这时候EKF的雅可比矩阵$F$和$H$的维数变成$2m \times 2m$手推会非常痛苦。我建议在多机场景下直接改用数值差分雅可比或者直接上UKF——不需要推导雅可比只需要把每台机的量测方程原样写进$h$函数即可。7.2 计算效率的取舍多机系统中UKF每拍要做$4m1$次非线性函数传播。对IEEE 39节点系统10台发电机就是每拍41次函数传播单次传播又是$2m$维的向量运算。实测下来在普通PC上每拍耗时大约1~2毫秒仍然在实时运行可接受范围内。如果系统规模继续扩大可以考虑降阶方案比如仅对关键机组做动态状态估计或者用自适应sigma点策略来降低采样数量。7.3 和静态估计的配合使用在实际调度系统里动态状态估计不会完全取代静态状态估计两者是配合关系。静态估计提供全网的稳态工况基准动态估计负责在动态过程中跟踪关键机组的状态。我见过一些工程实现用静态估计的结果作为动态估计的初始值再用PMU量测驱动动态估计不断校正两者错开时间尺度工作效果很好。这个思路在做研究课题时也很有参考价值。8. 最后想说的几句实在话整套代码调试下来我最大的体会是EKF和UKF的数学推导没有多难真正决定估计效果的是系统建模的准确性和噪声参数的标定。把$Q$和$R$调明白比换一个更高级的滤波器算法管用得多。再分享一个小技巧调试卡尔曼滤波时一定要把“真实轨迹—量测—估计值”三条曲线叠加在一张图上用眼睛看不要只看RMSE指标。RMSE是标量看不出滤波器的动态行为——你是否跟上了暂态过程的第一拍是否在稳态时还有残余的周期性误差这些只有看曲线才能发现。我见过不少论文里的RMSE指标很漂亮但曲线一放大全是锯齿状的高频跟随这种“假收敛”在实际系统中根本不敢用。代码本身的扩展性也不容忽视。建议把h_measurement和f_dynamics做成真正的函数句柄传参这样后续不管换SMIB还是IEEE 39节点模型滤波器主循环一行都不用改。我用这套结构从两状态单机一路改到十状态多机总共只改了参数文件和两个函数体省了大量重复工作。希望这篇分享能给你省下同样的时间。