电力系统动态状态估计：EKF与UKF的Matlab实现与调参指南

最近在帮几个师弟调电力系统动态状态估计的仿真代码，发现大家卡住的位置都差不多——模型都能写出来，卡尔曼增益公式也背得熟，但真到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，ω0=2π×50。
• 发电机用经典二阶模型，E'、X'd、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回避了雅可比，它的核心是构造2n+1个sigma点，n是状态维数。设状态均值x̄、协方差P，使用比例对称采样：

λ = α²(n + κ) - n

S = chol((n + λ)P)，取Cholesky分解的下三角矩阵

sigma点集合为

χ₀ = x̄，χ_i = x̄ + S_i，χ_{i+n} = 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 - ẑ)(...)ᵀ + R

P_{xz} = Σ W_i^c (χ_i - x̂)(ζ_i - ẑ)ᵀ

然后增益K = P_{xz} P_{zz}^{-1}，状态和协方差更新。

从实现角度看，UKF比EKF多的工作主要集中在sigma点生成和两组权重数组管理上，核心计算量约是EKF的三倍左右。但换来的是不需要任何导数信息，在量测函数很复杂的时候优势极其明显。

2.3 两者在电力系统DSE场景下的对比

我在这次仿真里的实测结论是，单机系统正常工况下两者精度差异不大，EKF甚至因为数值噪声更小显得更稳。但把扰动加大、功角摆动幅度拉高后，UKF的跟踪优势就开始体现出来。如果以后扩展到大电网多机系统，我会优先考虑UKF。

3. 电力系统动态模型构建与离散化

3.1 发电机二阶状态方程的建立

单机无穷大系统的经典二阶模型写成连续形式：

dδ/dt = ω - ω0

dω/dt = (Pm - Pe - D(ω - ω0)) / M

其中M = 2H/ω0，H为惯性时间常数，D为阻尼系数，Pm为机械功率，Pe为电磁功率。在经典模型中，Pe由暂态电势E'、机端电压Vt和直轴暂态电抗X'd决定：

Pe = (E' Vt / X'd) 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 / X'd) sin δ + v₃

把第三项加进去，量测方程就变成非线性的了，EKF在更新步就必须要算h对x的雅可比，这个设计能更加真实地考验UKF的值逼近能力。

量测噪声标准差我设成pa.delta=0.005rad、pa.omega=0.005pu、pa.Pe=0.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(:, k+1) = 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(:,k+1) = x_pred + K_k * (z_meas(:,k+1) - z_pred); P_ekf = (eye(2) - K_k * H_k) * P_pred; end

这里我用的是K = P_pred * H' / S_k，Matlab里“/”等价于右除乘以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]; % 2n+1列，每列是一个sigma点

权重数组提前算好：

n = 2; alpha = 1e-2; beta = 2; kappa = 0; lambda = alpha^2 * (n + kappa) - n; Wm = ones(2*n+1, 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*n+1); for j = 1:2*n+1 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*n+1); for j = 1:2*n+1 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(:,k+1) = x_pred + K * (z_meas(:,k+1) - 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]

从q=1e-4开始，依次试到q=1e-8，观察滤波曲线的平滑程度与延迟。q大了滤波器会过度信任量测，跟着噪声走，曲线毛刺多；q小了滤波器过度信任模型，真实轨迹发生突变时跟踪滞后明显。最后我在单机系统里取q=1e-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.2s，EKF的跟踪曲线出现了一个较为明显的振荡，UKF则更平滑地跟上了过渡过程。

这个结果完全符合理论预期——阶跃扰动让功角轨迹的局部非线性增强，一阶线性化的EKF在这个时刻出现精度损失，而sigma点传播的UKF对非线性分布拟合得更好。角速度估计上两种滤波器的差距相对小一些，因为ω的动态方程在扰动后很快回到线性区域。

这张表可以直观看到两种算法在相同参数下的精度数据（模拟结果）：

单步耗时是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点传播循环尽量向量化。我的经验是，对n=2的简单模型，向量化能快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，你会发现后面这个无非就是“换了一种算概率分布的方式”，架构上并没有本质障碍。