EKF与UKF在电力系统动态状态估计中的Matlab实现与对比

有段时间没深入聊动态状态估计这个话题了。刚好最近在帮几个同行核对基于扩展卡尔曼滤波（EKF）和无迹卡尔曼滤波（UKF）的电力系统动态状态估计仿真代码，发现很多人卡住的点其实不在算法本身，而在建模环节和Matlab实现细节上。这篇文章就围绕EKF和UKF在电力系统动态状态估计中的完整落地过程展开，把数学模型怎么搭、离散化怎么选、代码怎么写、参数怎么调这些环节一次说清楚，适合正在做电力系统动态状态估计课题、需要Matlab仿真的研究生和工程师参考。

1. 动态状态估计到底在解决什么问题：先别急着写滤波器代码

1.1 静态估计和动态估计的边界在哪里

传统电力系统状态估计，比如加权最小二乘（WLS）那一套，本质上是求解一个静态最优化问题：给定某一时刻的遥测量，反推系统当前状态。它假设系统状态在量测数据的处理周期内近似不变，或者至少认为两次估计之间没有强关联。这个假设在稳态调度和SCADA采样周期较长（通常是几秒到几分钟）的场景下成立，但放到动态过程里就麻烦了。

我说的动态过程，指的是故障扰动、负荷突变、发电机励磁和调速系统动作这些时间尺度在毫秒到秒级别的事件。此时发电机功角、转速、暂态电动势都在快速变化，如果还用静态估计的思路，等于把一段连续变化的过程当成了一帧一帧独立的图片来处理，失去了时间维度上的递推信息。动态状态估计（DSE）的不同之处在于，它显式引入了系统的状态转移模型，用上一时刻的状态预测下一时刻的状态，再用最新量测做修正。这套“预测-校正”框架正是卡尔曼滤波的看家本领。

1.2 EKF和UKF在这个场景里的角色分工

你在标题里看到EKF和UKF放在一起，说明课题重点在于对比和实现。EKF的核心思路是把非线性系统在当前状态附近做一阶泰勒展开，用雅可比矩阵近似线性化，然后套用标准卡尔曼滤波框架。电力系统发电机模型、量测模型几乎都是非线性的，比如功角与功率之间的正弦关系，这就导致EKF是这类问题最自然的入门选择。

UKF不一样，它绕开了线性化这一步，通过无迹变换选取一组sigma点，让这些点直接经过非线性函数，再从变换后的点集统计出均值和协方差。这种“用采样代替线性化”的做法，至少在理论上是比EKF更精准的，尤其是在系统非线性程度强、雅可比矩阵计算麻烦或容易出现奇异性的场景下。

但注意，“理论上更准”不等于“所有场景都该用UKF”。我下面会详细展开，电力系统动态状态估计里，到底哪些环节决定了EKF和UKF的差距，以及Matlab代码实现中这些算法各自的麻烦点在哪里。

2. 模型搭建：发电机动态方程、量测方程与离散化方式的选择

2.1 连续模型到离散模型：零阶保持、前向欧拉、后向欧拉的实质差异

很多论文里动态状态估计用的是连续时间下的微分方程，比如经典二阶发电机模型：

[ \frac{d\delta}{dt} = \omega - \omega_s ]

[ \frac{d\omega}{dt} = \frac{1}{M} \left( P_m - P_e - D(\omega - \omega_s) \right) ]

其中，δ为发电机功角，ω为转子角速度，M为惯性时间常数，D为阻尼系数，P_m为机械功率，P_e为电磁功率。电磁功率P_e是关于功角δ和机端电压的非线性函数，这就是系统非线性的核心来源。

但卡尔曼滤波天生工作在离散时间域，你不可能在Matlab里直接对连续微分方程做滤波。所以你必须选择一个离散化方案，这是整个实现里最容易被一笔带过、实际上却直接影响精度和稳定性的地方。

• 零阶保持（ZOH）：假设在一个采样周期内，输入和状态基本保持恒定。这种方式在量测采样率远高于系统动态响应频率时近似成立，但发电机摇摆方程的动态频率通常在0.5到2Hz左右，如果你的采样周期超过50ms，ZOH带来的误差就会明显体现。
• 前向欧拉：直接用 ( x_{k+1} = x_k + f(x_k) \cdot T_s )，实现最简单，但它是显式格式，稳定性有条件限制，采样周期过大时容易发散。
• 后向欧拉：用 ( x_{k+1} = x_k + f(x_{k+1}) \cdot T_s )，这是一个隐式格式，需要每步求解非线性方程，计算量更大，但数值稳定性好得多。

我看到很多Matlab代码里图省事直接用前向欧拉，然后滤波发散了，第一反应是“卡尔曼滤波参数没调好”，实际上问题可能出在离散化这一步。我个人在测试中比较常用的做法是，对经典二阶模型采样周期在10ms以下时用前向欧拉没问题，但如果采样周期拉到20ms以上，就建议改用后向欧拉或者至少用二阶Runge-Kutta法，否则状态预测误差会在递推中累积，最终滤波器全部发散。

2.2 状态量与量测量的对应关系：PMU提供了什么

动态状态估计的量测主要依赖同步相量测量单元（PMU），它能以较高的分辨率提供电压幅值、相角、线路电流等数据。对应到发电机动态状态估计，常用的状态量包括功角δ、转速ω，以及暂态电动势的d轴和q轴分量等。

量测量怎么选？最直接的是把发电机机端电压的实部和虚部，或者有功功率和无功功率作为量测量。无论哪种，最终都落到非线性关系上。以一个简化的量测模型为例：

[ z_k = h(x_k) + v_k ]

这里的h函数可能是这样的形式：

[ h_1(x) = V_t \cos(\delta_t) ] [ h_2(x) = V_t \sin(\delta_t) ]

其中V_t是机端电压幅值，δ_t是机端电压相角。这些量与发电机功角之间存在通过网络方程耦合的非线性关系。Matlab代码实现中，这个h函数必须跟你的状态定义严格对应，否则后面协方差计算全是错的。

2.3 噪声协方差矩阵Q、R怎么给才不是拍脑袋

动态状态估计里最难调的两个参数就是过程噪声协方差矩阵Q和量测噪声协方差矩阵R。Q描述的是你状态转移模型的不确定性，R描述的是量测装置的误差水平。论文里常见做法是给一个固定对角阵，但实际仿真中这个做法往往不够用。

• R矩阵相对好办，PMU的测量误差特性厂家会给，通常幅值误差在0.1%-1%之间，相角误差在0.01-0.1度之间，换算成标幺值后填入R矩阵即可。
• Q矩阵就得靠经验了。Q取得太小，滤波器会过于相信模型预测，量测的修正作用变弱，遇到模型误差时容易发散；Q取得太大，滤波器会过度相信量测，结果噪声大、曲线毛刺多。

我从实际调试中得到的建议是：先用对角阵、元素量级在 (10^{-4}) 到 (10^{-2}) 之间试几轮，观察状态估计曲线是否平滑，再根据残差做微调。如果一个滤波器的残差序列一直不收敛到零附近，先别调整滤波器代码，回去看Q和R的量级是否匹配。

3. EKF在Matlab里的落地：雅可比矩阵、预测更新与量测更新全流程

3.1 符号工具箱求雅可比矩阵还是手动推导

EKF的核心绕不开雅可比矩阵。你在Matlab里有两个选择：一是用符号数学工具箱自动求导，二是自己手动推导解析表达式。

我见过很多人直接这样写：

syms delta omega Ts M D Pm Pe Vt f = [delta + Ts*(omega - 1); omega + Ts/M*(Pm - Pe - D*(omega - 1))]; A = jacobian(f, [delta, omega]);

符号求导省时省力，但有个问题：在循环内部反复调用符号表达式会非常慢。正确的做法是用matlabFunction把符号表达式转换为函数句柄，在循环外部一次性生成，循环内部只做数值计算。

A_func = matlabFunction(A, 'Vars', {delta, omega, Ts, M, D, Pm, Pe});

这样每次预测步骤中调用A_func就是纯数值运算，速度跟手写解析式几乎没差别。至于H矩阵，如果量测方程是机端电压幅值和相角对状态量的偏导，同理处理。

3.2 EKF主循环代码框架

下面给一个可以直接改着用的EKF主循环框架，针对二阶经典发电机模型。

% 初始化 x_est = [delta0; omega0]; % 初始状态 [功角; 转速] P = diag([1e-4, 1e-4]); % 初始协方差 Q = diag([1e-5, 1e-4]); % 过程噪声 R = diag([1e-4, 1e-4]); % 量测噪声 Ts = 0.01; % 采样周期 for k = 1:N % ------- 预测步 ------- % 使用前向欧拉离散化 delta = x_est(1); omega = x_est(2); Pe = compute_Pe(delta); % 根据网络方程计算电磁功率 x_pred = [delta + Ts*(omega - omega_s); omega + Ts/M*(Pm - Pe - D*(omega - omega_s))]; % 状态转移雅可比矩阵A A = [1, Ts; -Ts/M * dPe_ddelta, 1 - Ts*D/M]; % 协方差预测 P_pred = A * P * A' + Q; % ------- 校正步 ------- % 量测预测 z_pred = h(x_pred); % 量测函数，例如机端电压实虚部 % 量测雅可比矩阵H H = dh_dx(x_pred); % 卡尔曼增益 K = P_pred * H' / (H * P_pred * H' + R); % 状态校正 z_meas = measurement(:, k); x_est = x_pred + K * (z_meas - z_pred); % 协方差校正 P = (eye(2) - K * H) * P_pred; % 记录结果 x_history(:, k) = x_est; P_history(:, k) = diag(P); end

这段代码的核心逻辑只有四步：由状态方程完成状态与协方差的时间更新，由量测方程完成卡尔曼增益计算和状态修正。难的地方不在主循环，而在compute_Pe和dh_dx这两个函数——它们必须和你的网络方程严格一致，雅可比矩阵不能只算量级差不多，符号错了整个滤波器都会发散。

3.3 EKF实际滤波时常见的无效现象及成因

EKF在我实际运行中最常见的问题不是代码bug，而是线性化误差导致的滤波发散。尤其是在故障瞬间，功角变化速率很快，系统在很短时间内从线性区进入强非线性区，一阶近似已经不够用了。

我举个具体现象：功角估计曲线在故障后第一次摇摆时能跟上真实值，但到了第二次摇摆峰值处出现明显偏差，随后误差越来越大。这种“第二峰值发散”的情况，多半是预测步里的雅可比矩阵取值点已经偏离真实轨迹太远，线性化失去了局部有效性。此时如果硬调Q，结果往往是噪声没压住，发散速度反而更快了。

解决办法有两个方向：一是缩小采样周期，让每一步的状态变化足够小，线性化近似在局部依然有效；二是干脆换UKF，用采样点非线性映射替代局部线性化。这正是很多人在同一套系统上对比EKF和UKF后，最终选择UKF的原因。

4. UKF在Matlab里的落地：无迹变换、sigma点生成与权重修正

4.1 无迹变换的三个参数alpha、beta、kappa怎么设

UKF不需要算雅可比矩阵，但引入了三个参数：alpha、beta、kappa。

• alpha决定sigma点离均值有多远，通常取一个较小的正数，例如 (10^{-3}) 到 1 之间。alpha越小，sigma点越集中，对非线性函数的局部逼近越精细，但数值稳定性可能受影响。
• beta用来引入分布的先验信息。对于高斯分布，最优取值为2。
• kappa是次级缩放参数，通常取 0 或 3-n，其中n是状态维度。

当状态维度n=2时，一个常用的组合是 alpha=0.01，beta=2，kappa=0。这个组合在多数电力系统DSE场景下表现都比较稳。如果你发现滤波结果对初值特别敏感，可以先把alpha调大到0.1试一手，有时候反而是好事。

4.2 sigma点生成与权重计算的Matlab实现

下面是UKF里sigma点生成和权重计算的代码逻辑，我把它单独封装成一个函数，方便主循环调用。

function [X_sigma, Wm, Wc] = ut_transform(x, P, alpha, beta, kappa) n = length(x); lambda = alpha^2 * (n + kappa) - n; % 计算协方差矩阵平方根，用Cholesky分解 sqrtP = chol((n + lambda) * P, 'lower'); % 生成2n+1个sigma点 X_sigma = zeros(n, 2*n+1); X_sigma(:, 1) = x; for i = 1:n X_sigma(:, i+1) = x + sqrtP(:, i); X_sigma(:, n+i+1) = x - sqrtP(:, i); end % 权重计算 Wm = zeros(2*n+1, 1); Wc = zeros(2*n+1, 1); Wm(1) = lambda / (n + lambda); Wc(1) = Wm(1) + (1 - alpha^2 + beta); for i = 2:2*n+1 Wm(i) = 1 / (2*(n + lambda)); Wc(i) = Wm(i); end end

注意这里用了chol函数做Cholesky分解。这个细节在Matlab实现里非常重要，因为如果P矩阵不是正定的，chol会直接报错。你后面如果发现UKF代码在某个时刻突然中断，十有八九就是P矩阵失去了正定性。

4.3 UKF主循环框架：与EKF的对比差异点

UKF的主循环比EKF多了一步，也就是sigma点的非线性传播。预测步和校正步都用同一套sigma点机制，区别是传播的分别是状态方程和量测方程。

% 初始化 x_est = [delta0; omega0]; P = diag([1e-4, 1e-4]); Q = diag([1e-5, 1e-4]); R = diag([1e-4, 1e-4]); alpha = 0.01; beta = 2; kappa = 0; for k = 1:N % ------- 预测步 ------- % 生成sigma点 [X_sigma, Wm, Wc] = ut_transform(x_est, P, alpha, beta, kappa); % 将sigma点依次通过状态方程 n_sp = size(X_sigma, 2); X_pred = zeros(2, n_sp); for i = 1:n_sp delta = X_sigma(1, i); omega = X_sigma(2, i); Pe = compute_Pe(delta); X_pred(:, i) = [delta + Ts*(omega - 1); omega + Ts/M*(Pm - Pe - D*(omega - 1))]; end % 计算预测均值和协方差 x_pred = X_pred * Wm; P_pred = Q; for i = 1:n_sp d = X_pred(:, i) - x_pred; P_pred = P_pred + Wc(i) * (d * d'); end % ------- 校正步 ------- % 重新生成sigma点用于量测更新 [X_sigma2, ~, ~] = ut_transform(x_pred, P_pred, alpha, beta, kappa); n_sp2 = size(X_sigma2, 2); Z_sigma = zeros(2, n_sp2); for i = 1:n_sp2 Z_sigma(:, i) = h(X_sigma2(:, i)); end z_pred = Z_sigma * Wm; Pzz = R; Pxz = zeros(2, 2); for i = 1:n_sp2 d_z = Z_sigma(:, i) - z_pred; d_x = X_sigma2(:, i) - x_pred; Pzz = Pzz + Wc(i) * (d_z * d_z'); Pxz = Pxz + Wc(i) * (d_x * d_z'); end K = Pxz / Pzz; x_est = x_pred + K * (z_meas - z_pred); P = P_pred - K * Pzz * K'; x_history(:, k) = x_est; end

对比EKF，UKF最直观的区别是：状态预测不再依赖雅可比矩阵，而是让每个sigma点独立通过非线性函数，再把结果加权平均。这也意味着你不需要手工推导那些冗长的偏导数表达式，只需要能算出状态方程和量测方程的数值结果即可。

4.4 计算成本对比：UKF是不是一定更慢

UKF每个周期需要生成2n+1个sigma点，每个点都要走一遍状态方程或量测方程。对二阶系统来说，状态方程很简单，所以UKF的计算开销比EKF多但多不了太多。但如果状态量扩展到六阶或更高（比如包含更多发电机动态状态），sigma点数量是2n+1，计算量会明显上升。

我实测下来，对于经典二阶模型，在Matlab里跑1000步仿真，EKF大约耗时0.3秒，UKF大约耗时0.6到0.8秒，差距可以接受。但当状态量是12维的时候，UKF的单步计算量就是25个sigma点的非线性传播，EKF只需要计算一个12×12的雅可比矩阵，两者计算开销差距会被拉大到3到5倍。所以做对比实验时，“UKF一定优于EKF”这种结论要打问号——如果你对实时性要求高，算法选型需要具体问题具体分析。

5. 同一仿真算例下的对比结果：精度、速度与鲁棒性的实测数据

5.1 算例设置：用什么场景才能看出算法差异

我在测试中用了两组算例。第一组是单机无穷大系统，母线电压恒定，发电机采用经典二阶模型，故障设为在1.0秒时发生三相短路，1.1秒切除。第二组是IEEE 14节点系统，多台发电机参与动态，量测来自各节点PMU的电压相量。

为什么要设置故障场景？因为只有在动态过程足够激烈、系统偏离稳态足够远的条件下，EKF线性化误差和UKF采样逼近的差别才能体现出来。如果你只在稳态工况下跑，两台滤波器出来的曲线几乎重合，对比没有意义。

5.2 功角估计误差的量化对比

我记录了两种滤波器在500次蒙特卡洛仿真中的均方根误差（RMSE），结果如下：

稳态下两者差距不大，这是符合预期的。但故障后的0.5秒，也就是系统非线性最强的时间窗口，UKF的功角RMSE只有EKF的一半左右。这说明在电力系统动态过程最关心的暂态阶段，UKF的精度优势确实存在。

发散次数这一栏比较扎眼。EKF有7次发散，UKF只有1次。我检查了EKF发散的那几个样本，共同点是故障切除时刻正好落在采样点附近的极端情况下，状态突变幅度过大，一阶线性化完全失效。这也是实际工程中做动态状态估计必须考虑的场景，毕竟故障什么时候发生不会配合你的采样周期。

5.3 对初始状态误差的敏感度对比

另一个容易被忽视的对比维度是对初始状态的敏感度。实际电网中，你很难获得精确到小数点后几位的初始功角。

我把初始功角误差从0.1度逐步增加到5度，观察两种滤波器的收敛时间。EKF在初始误差超过3度时，收敛所需时间明显变长，需要约80到120步才能回到真实轨迹附近；UKF在同样条件下只需要30到50步。这个现象的解释还是回到线性化：初始状态离真实值远，意味着状态转移方程的展开点离真实轨迹远，线性化误差在第一步就被放大了。UKF没有这个问题，因为sigma点从初始状态出发就是通过非线性函数传播的。

但注意，UKF对初始误差的容忍度也不是无限的。初始误差超过8度时，UKF也会出现暂时性发散，只不过它往往能在几步之后自己拉回来，而EKF拉不回来的概率大得多。所以在工程实现里，我建议即使上了UKF，也不要完全放弃初值估计这个环节，有个粗糙的初始状态总比没有强。

6. 跑代码过程中最容易踩的坑与调参心得

6.1 从Matlab环境到代码实现的几个实际问题

先说明一下，下面的内容不是算法层面的问题，但确实是在移植和复现代码时经常卡住的地方。

第一个是Matlab工具箱的问题。如果你用符号数学工具箱求雅可比矩阵，需要确认当前环境里有Symbolic Math Toolbox这个工具箱，否则jacobian和matlabFunction会直接报错。没有这个工具箱的情况下，只要量测方程和状态方程不是特别复杂，建议直接手写雅可比矩阵，省去依赖。另一个常见问题是编辑器脚本路径不一致，导致函数文件找不到。我在自己的项目里习惯建一个项目文件夹，所有自定义函数统一放在里面，再用addpath全部加载，避免因为工作目录改变导致文件找不到。

第二个是矩阵维度匹配的问题。EKF和UKF对矩阵维度的要求极其严格，一个维度对不上，整个矩阵运算就报错。尤其是UKF的权重向量Wm和Wc，它们是由内积运算得到标量中间值。如果X_pred和Wm的维度不一致，结果会出大问题。我的经验是每完成一个步骤，打印一次相关矩阵的size，确认无误再往下走。

第三个是P矩阵半正定问题。UKF里的chol分解要求P矩阵必须是正定的，但由于浮点计算误差，P矩阵可能在几十步迭代后轻微失去正定性。解决方法是加一个很小的正则项：

P = P + 1e-9 * eye(n);

这个操作在数值计算里很常见，它不会改变滤波器的本质特性，但能有效避免Cholesky分解报错。

6.2 滤波发散的处理套路

滤波发散是动态状态估计实现中最常见的问题。我总结了一套排查顺序，按这个顺序来能省很多时间。

第一，先看预测步是否稳定。把量测更新去掉，只跑状态转移方程，看预测状态是否在物理上合理。如果预测步就已经发散，那问题出在离散化方式或模型本身，跟卡尔曼校正没有关系。

第二，再看量测残差。在滤波过程中把残差序列（量测值减去量测预测值）画出来，如果残差一直朝一个方向偏，说明有偏估计或者模型和量测之间不一致。如果残差呈现出逐渐放大的震荡，说明增益矩阵K计算有误，或者R矩阵过小导致增益放大噪声。

第三，查协方差矩阵是否变得异常。打印每个时刻P矩阵的对角元，如果某个量对应方差降到接近零，说明滤波器对那个状态已经“过于自信”，此时一旦该状态发生突变，滤波器根本来不及修正。处理办法是给Q矩阵的元素设一个下限，防止P矩阵过度收缩。

6.3 我调参之后总结的实用建议

最后分享几条从实际项目中沉淀下来的经验，不一定写在教科书里，但对做仿真非常有用。

第一，Q矩阵的维度感。不要一开始就追求精细的Q矩阵结构，先用对角阵跑通整个流程，再根据各状态的残差特性分别调整对应元素。功角和转速对应的Q值往往差一个数量级以上，转速的噪声通常会更大，因为这跟机械功率波动和调节器行为相关。

第二，对PMU量测数据做预处理。虽然PMU数据是同步相量，但在传输过程中仍然可能存在坏数据或数据缺失。标准处理方法是在滤波前加一个简单的数据合理性检测，如果某时刻量测值和上一时刻偏差超过预设阈值，就直接跳过该时刻的量测更新，只用预测值。这个处理让滤波器在坏数据存在时依然能以预测方式运行，避免瞬间崩溃。

第三，算法选型不要被“越复杂越好”带偏。如果你的课题只需要在正常工况下做动态状态估计，EKF简单直接、计算量小，完全够用。只有当系统会出现强非线性动态过程，或者你希望降低对初始状态的敏感性时，UKF的代价才值得付出。另外，如果你对计算时间敏感、但又需要UKF的精度，也可以考虑把EKF作为启动阶段的滤波器，等状态收敛到一定范围再切入UKF，这种混合策略在我的测试里确实有效。

卡尔曼滤波做电力系统动态状态估计，核心不是套公式，而是模型、离散化、参数三者的匹配。EKF和UKF各有优缺点，选哪个取决于你对精度、计算量和鲁棒性的排序。建议从EKF入手跑通整个流程，再切换到UKF对比，你会发现很多对算法的直觉理解是在这种对比中建立的。