电力系统动态状态估计：EKF与UKF的Matlab实现全解析

电力系统动态状态估计这个方向，我在Matlab里从静态加权最小二乘一路折腾到EKF、UKF，最大的感受是：很多资料要么只讲公式，要么直接丢代码，真正把“为什么这么做”和“怎么调通”讲透的不多。这篇我打算把基于EKF和UKF的电力系统动态状态估计从建模、参数设置、代码实现到问题排查完整梳理一遍，尤其是那些仿真时最容易踩的坑，都会用实际经验说明白。想复现这套Matlab代码、做毕业设计或者入门动态状态估计的朋友，按照这趟流程走一遍，基本能把整个逻辑串起来。

在进入正题之前，先明确一下我们到底在做什么。电力系统动态状态估计，简单说就是利用带有噪声的量测数据（比如PMU的电压幅值、相角、功率等），递归推算出系统当前真实的动态状态（比如发电机功角、转速、暂态电动势）。它和传统静态状态估计最大的区别在于引入了时间维度的系统动态模型，而EKF和UKF正是处理这种非线性递推问题最直观的两种滤波器。

1. 电力系统动态状态估计：从静态到动态，问题本身发生了什么变化

1.1 静态状态估计（WLS）的局限：为什么我们非要用卡尔曼滤波

很多教材会把电力系统状态估计讲成“加权最小二乘”的死忠，以至于不少人以为状态估计就是求解一个优化问题。但实际接触过电力系统动态过程的都知道，系统是一个大规模非线性动态系统，不是一张静态断面图。WLS做的是在某一时刻，利用当前时刻的SCADA/PMU量测，找一个最符合量测方程的静态工作点。

问题就在这里：WLS把各个时刻的估计结果看作是相互独立的，没有利用系统自身的运动规律。对于稳态工况，它确实够用，但一旦系统发生扰动（比如切负荷、发电机跳闸或者励磁系统动作），功角、转速、内电动势的变化是连续过程，仅仅靠一个静态快照去拟合，需要等到量测全部更新、迭代收敛才能给出结果，响应速度和滤波平滑效果都不理想。

动态状态估计的思路则是：上一时刻系统的状态通过动态方程演化到当前时刻，再用量测数据做修正。这个过程天然适合卡尔曼滤波框架——预测加更新。而电力系统的动态方程和量测方程又都是非线性的，所以才演化出了EKF（对非线性做局部线性化）和UKF（用采样点近似概率分布）这两条主流路线。在PMU数据逐渐普及的背景下，量测刷新速度从秒级提升到毫秒级，动态估计的价值也随之凸显。

1.2 动态估计的数学模型：状态方程加量测方程，一个都不能少

动态状态估计的出发点是一个离散化的非线性随机系统，一般写成：

[ x_k = f(x_{k-1}) + w_k ] [ z_k = h(x_k) + v_k ]

其中 (x_k) 是第 (k) 个采样时刻的系统状态向量，(z_k) 是量测向量，(w_k) 是过程噪声，(v_k) 是量测噪声。这里的 (f) 是动态转移函数，通常由发电机转子运动方程和励磁绕组动态方程离散化得到；(h) 是量测函数，把状态量映射到量测值上。

我在实际建模时习惯把 (f) 写成连续微分方程，再用数值积分方法离散化。以最常用的三阶单机模型为例：

[ \dot{\delta} = \omega - \omega_s ] [ \dot{\omega} = \frac{\omega_s}{2H}(P_m - P_e - D(\omega - \omega_s)) ] [ \dot{E'q} = \frac{1}{T'{d0}}(E_{fd} - E'_q - (x_d - x'_d)I_d) ]

这里 (\delta) 是功角，(\omega) 是角速度，(E'q) 是暂态电动势，(H) 是惯性常数，(D) 是阻尼系数，(P_m) 是机械功率，(P_e) 是电磁功率，(T'{d0}) 是励磁绕组时间常数。离散化步长通常取采样间隔 (\Delta t)，在PMU场景下一般是 (0.01s \sim 0.02s)。

如果从椅子上站起来做类比，状态方程就像“人的运动趋势”——你知道一个人上一秒在什么位置、速度是多少，就能粗略推断下一秒大概在哪；量测方程则像“摄像头拍到的画面”——不一定完全准确，但能修正你的推断。卡尔曼滤波干的就是把这俩拧到一块，既相信模型趋势，又相信量测修正，最终得到一个比两者单独使用都更准的状态估计。

1.3 采样步长和动态模型的选择，直接决定滤波能不能跟上真实轨迹

步长选择是个很实际的问题。如果PMU是 30帧/秒、60帧/秒，采样间隔只有几十毫秒，那么系统状态在这么短时间内的变化量通常不大，用一阶欧拉离散化问题不大。但如果你处理的是仿真数据或者SCADA数据，采样间隔到了几百毫秒甚至秒级，一阶欧拉误差就会明显偏大，我一般会换成四阶Runge-Kutta（RK4）来对 (f) 积分。

模型选择上，经典二阶模型（只含 (\delta) 和 (\omega)）最简单，适合做算法验证；三阶模型加了 (E'_q) 动态，能反映励磁系统影响，是目前研究和工程文献中最常用的折中方案；更详细的四阶、六阶模型虽然贴近实际，但对EKF来说雅可比矩阵推导会非常痛苦，而且对参数精度要求更高，新手不建议一上来就上高阶模型。

2. EKF和UKF的核心原理拆解：线性化和无损变换，路径不同但目标一致

2.1 EKF：对非线性系统进行局部线性化，用泰勒展开一阶近似

EKF的思路很直接：既然卡尔曼滤波只适用于线性系统，那我就在当前估计点附近把非线性函数展开、丢掉高阶项，得到一个近似的线性系统，然后套标准卡尔曼滤波公式。

具体来说，在每个采样时刻：

1. 计算状态转移矩阵 (F_{k-1} = \frac{\partial f}{\partial x}\big|{x=\hat{x}{k-1}})，以及量测矩阵 (H_k = \frac{\partial h}{\partial x}\big|{x=\hat{x}{k|k-1}})。
2. 状态预测：(\hat{x}{k|k-1} = f(\hat{x}{k-1}))。
3. 协方差预测：(P_{k|k-1} = F_{k-1}P_{k-1}F_{k-1}^T + Q)。
4. 增益计算：(K_k = P_{k|k-1}H_k^T(H_kP_{k|k-1}H_k^T + R)^{-1})。
5. 状态更新：(\hat{x}k = \hat{x}{k|k-1} + K_k(z_k - h(\hat{x}_{k|k-1})))。
6. 协方差更新：(P_k = (I - K_kH_k)P_{k|k-1})。

用起来最大的难点就在 (F) 和 (H) 的求导上。对于一个十几维的三阶多机状态空间，手推雅可比矩阵很容易出错。而且EKF的本质是用一个线性近似去代替真实非线性映射，当系统在强非线性区域（比如故障后功角摆开幅度很大）工作时，一阶近似误差会放大，滤波精度会显著下降，甚至发散。

我之前做一个双机三阶模型时，手推 (H_k) 矩阵整整花了半天，结果仿真中还是出现了估计偏差逐渐增大的情况，后来一核对，是量测方程对 (E'_q) 的偏导漏了一项。这种错误在EKF实现里极其常见。

2.2 UKF：用无迹变换（UT）代替线性化，用采样点近似分布

UKF绕开了“求雅可比”这个痛点，核心是无迹变换（Unscented Transform）。它的思想是：对一个随机向量，与其直接线性化非线性函数，不如在原分布中选取一组精心构造的Sigma点，把这些点通过非线性函数传播，再从传播后的点集统计出均值和协方差。

对 (n) 维状态向量，通常选取 (2n+1) 个Sigma点：

[ \chi_0 = \bar{x} ] [ \chi_i = \bar{x} + \left(\sqrt{(n+\lambda)P_x}\right)i, \quad i=1,\dots,n ] [ \chi_i = \bar{x} - \left(\sqrt{(n+\lambda)P_x}\right){i-n}, \quad i=n+1,\dots,2n ]

对应的权重为：

[ W_0^{(m)} = \frac{\lambda}{n+\lambda} ] [ W_0^{(c)} = \frac{\lambda}{n+\lambda} + (1 - \alpha^2 + \beta) ] [ W_i^{(m)} = W_i^{(c)} = \frac{1}{2(n+\lambda)}, \quad i=1,\dots,2n ]

其中 (\lambda = \alpha^2(n+\kappa) - n)。(\alpha) 控制Sigma点离均值的距离，通常取 (1e-3 \sim 1)；(\kappa) 取值通常为 (0)；(\beta) 用来融入先验分布信息，高斯分布下取 (2)。

这些Sigma点经过 (f) 和 (h) 的非线性传播之后，加权平均得到预测均值和预测协方差。整个过程不需要任何求导操作，而且对非线性函数的近似精度至少达到二阶，明显优于EKF的一阶线性化。

用生活化的比喻来说，EKF是“把一条弯曲的路每段都看成直线，找最短方向”，UKF则是“在路上安排一批观察员，分别走到不同位置，再汇总他们看到的情况”。对于电力系统这种强非线性、非线性的量测方程，UKF通常更稳。

2.3 EKF和UKF怎么选：算力、雅可比、非线性强度的综合权衡

我在项目里同时实现了两种算法，实际使用下来可以给出一个很直观的对比：

如果状态维度在10维以下，UKF的计算开销完全可以接受。我实测在普通笔记本上跑IEEE 14节点系统的三阶发电机模型，UKF单步运行时间大约比EKF慢1.5到2倍，但换来的是更好的数值稳定性和更高的估计精度。如果状态维度超过20维，UKF的Sigma点数量会明显上升，此时EKF可能更适合实时性要求高的场景。

我个人给初学者推荐的路径是：先用EKF跑通整条流程，感受一下预测-更新循环的每一步在干什么；再切换到UKF，你就能立刻体会到“不用推雅可比”的快乐。两种算法共用同一套状态方程和量测方程，切换起来并不困难。

2.4 从零复现的顺序建议：先单机后多机，先简化后精细

状态估计的仿真最容易犯的错就是把初始系统搞得太复杂。我见过不少同学一上来就选IEEE 118节点、六阶详细模型，结果连滤波器发散都找不到原因。建议复现顺序是：

1. 单机无穷大系统，三阶发电机模型，状态维度3，量测取有功、无功和机端电压幅值。先把EKF调通，观察收敛情况。
2. 单机系统换UKF，对比两种滤波器的估计误差。
3. 扩展到两机或三机系统，状态维度变成6或9，加入线路功率量测，测试多机间状态互相影响下的滤波性能。
4. 再上IEEE 14节点等标准算例，这时候状态维度一般在12到30左右，可以系统验证算法在大网下的表现。

我在第二步对比时就发现，UKF对初值和噪声参数的鲁棒性明显更强，EKF在初值偏差大时容易在一开始出现大误差尖峰，需要两三个采样周期才能拉回来。这个差异在写论文时是很自然的对比素材。

3. 核心细节与参数设置：状态选什么、量测怎么配、Q和R怎么给

3.1 状态量选取：三阶发电机模型是最常见的选择，维度适中

多数文献和Matlab实现项目采用的动态状态估计状态向量是基于发电机转子运动方程和励磁动态的三阶模型，即每个发电机取：

[ x_i = [\delta_i, \omega_i, E'_{q,i}]^T ]

如果有 (n_g) 台发电机，系统状态维度就是 (3n_g)。比如IEEE 14节点系统中通常有5台同步发电机，状态维度就是15。我还见过只取 (\delta_i, \omega_i) 的二阶模型，状态维度降到 (2n_g)，算法跑起来更轻快，但无法刻画励磁系统对内电动势动态的影响，估计精度有限。

选择状态量时要特别注意一个工程细节：功角 (\delta) 是相对于参考机或者无穷大母线的相对量。单机无穷大系统里直接以无穷大母线为参考，没有问题；多机系统里必须明确参考机，否则状态向量不唯一，滤波会跑偏。

3.2 量测配置：电压幅值、注入功率和线路功率的搭配

量测向量(z_k)设计直接决定了可观测性和估计精度。工程上常见的量测包括：

• 发电机机端电压幅值 (V_i)。
• 发电机有功注入 (P_i) 和无功注入 (Q_i)。
• 关键母线电压幅值和相角（PMU量测）。
• 关键线路的潮流 (P_{ij}, Q_{ij})。

量测方程 (h(x)) 本质上是潮流方程的逆映射。以发电机节点为例，有功注入可以写成：

[ P_i = V_i \sum_{j \in N_i} Y_{ij} V_j \cos(\theta_i - \theta_j - \alpha_{ij}) ]

其中 (Y_{ij}) 是导纳矩阵元素，(\alpha_{ij}) 是导纳相角。由于状态量中包含 (E'_q)，机端电压幅值和相角需要通过与电抗和节点电压的关系进一步关联，这一块在Matlab实现里占代码比重最大。

我的经验是：量测冗余度要有，比如每个发电机的 (P, Q, V) 至少全部量测，才能保证滤波器不会退化；但量测数量也不宜过多到让计算量爆炸。先用“每个发电机节点配置有功、无功、电压幅值”这一套组合，基本够用。

3.3 协方差矩阵Q和R：整篇代码里最需要耐心调整的部分

过程噪声协方差 (Q) 反映模型的不确定性。模型越粗略、负荷波动越大，(Q) 应该给得越大。我常用对角形式：

[ Q = \text{diag}(q_{\delta}, q_{\omega}, q_{E'_q}, \dots) ]

对于三阶发电机模型，一个能工作的初始值是 (q_{\delta}=1e-6)，(q_{\omega}=1e-4)，(q_{E'_q}=1e-4)。量测噪声协方差 (R) 则可以依据传感器精度来定，PMU的幅值量测误差通常在 (0.1% \sim 1%)，相角误差在 (0.1^\circ \sim 1^\circ)。取电压幅值方差 (0.01^2)，有功无功功率方差 (0.02^2)，是个比较合理的起点。

Q和R都小、滤波会过度信任模型，导致量测修正很弱，输出曲线“过于平滑”但可能偏离真值；Q和R都大、滤波会剧烈抖动甚至发散。调参方法说白了就是：

1. 先把R设成真实噪声水平，只调Q。
2. 从极小值逐步增大Q，观察估计曲线是否在量测噪声中仍能平滑跟随真实轨迹。
3. 如果出现振荡型发散，优先减小Q中的转速项；如果出现滞后性误差，优先增大Q。

这个环节没有理论一步到位的方法，只能靠仿真观察。我当初在调EKF时，因为Q给得过小，结果滤波功能完全失效，状态轨迹几乎就是开环动态模型的预测，量测噪声根本没被压下去。

4. 实操全过程：从初始化到结果评估，一步步跑通Matlab代码

4.1 滤波器初始化：状态初值、协方差初值、噪声矩阵的确定

状态初值最好从一个潮流解获得。先运行MATPOWER或者自己写个牛顿-拉夫逊潮流，得到稳态的功角、机端电压和注入功率，然后反推 (E'_q) 的初值。如果图省事，也可以用“平坦启动”：

[ \hat{x}_0 = [0, 2\pi f_s, 1.0, 0, 2\pi f_s, 1.0, \dots]^T ]

其中功角初值给0，转速给同步转速（弧度每秒），(E'_q) 给1.0。平坦启动在单机系统中通常能正常收敛，多机系统里最好还是做个潮流初始化。

协方差初值 (P_0) 表示对初始状态的信任程度。给太大会导致前几步增益很高、状态跳变剧烈；给太小会导致滤波收敛极慢。我一般取：

[ P_0 = \text{diag}(0.1, 0.01, 0.1, \dots) ]

这个数量级下，功角和内电动势的不确定性稍大一些，转速不确定性小一些，比较符合实际物理直觉。

4.2 EKF主循环实现：离散化、雅可比和预测-更新全流程

EKF的核心循环可以用下面的Matlab代码结构来表示。以单机三阶模型为例，先定义状态方程函数：

function dx = gen_dyn(x, u, params) % x = [delta; omega; Eqp] % u = [Pm; Efd] delta = x(1); omega = x(2); Eqp = x(3); wb = params.wb; H = params.H; D = params.D; Td0 = params.Td0p; xd = params.xd; xdp = params.xdp; Id = Eqp / xdp; % 简化假设，忽略端电压影响 Pe = Eqp * 1.0 / xdp; % 实际应写全潮流表达式 dx = [ omega - wb; wb / (2*H) * (u(1) - Pe - D*(omega - wb)); (u(2) - Eqp - (xd - xdp)*Id) / Td0; ]; end

实际工程中 (P_e) 要写成与网络方程耦合的完整形式，这里只展示结构。离散化用RK4：

function x_next = rk4_step(fun, x, u, params, dt) k1 = fun(x, u, params); k2 = fun(x + 0.5*dt*k1, u, params); k3 = fun(x + 0.5*dt*k2, u, params); k4 = fun(x + dt*k3, u, params); x_next = x + (dt/6)*(k1 + 2*k2 + 2*k3 + k4); end

EKF的雅可比矩阵 (F) 和 (H)，如果在多机系统下手推太麻烦，可以用数值差分代替：

F = zeros(n, n); dx = 1e-6; for i = 1:n xp = x; xp(i) = x(i) + dx; xm = x; xm(i) = x(i) - dx; F(:, i) = (f_func(xp, u, params) - f_func(xm, u, params)) / (2*dx); end

数值差分求导在步长合适时精度完全够用，而且比我当时手推公式出错概率低得多。求完 (F) 和 (H) 后，套标准EKF公式执行预测、增益、更新即可。

4.3 UKF实现要点：Sigma点生成、非线性传播、协方差修复

UKF主循环比EKF更统一，因为它不需要区分 (F) 和 (H) 的求法，全部通过Sigma点传播。Matlab中可以这样写生成Sigma点：

function [X, Wm, Wc] = ut_sigma(x, P, alpha, beta, kappa) n = numel(x); lambda = alpha^2 * (n + kappa) - n; cov_sqrt = chol((n + lambda) * P, 'lower'); % 要求P正定 X = zeros(n, 2*n+1); X(:, 1) = x; for i = 1:n X(:, i+1) = x + cov_sqrt(:, i); X(:, i+n+1) = x - cov_sqrt(:, 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); Wm(2:end) = 1 / (2*(n + lambda)); Wc(2:end) = 1 / (2*(n + lambda)); end

注意chol要求协方差矩阵正定。滤波过程中由于数值误差，(P) 可能变得半正定，导致Cholesky分解失败。我通常在每次调用ut_sigma前加一个检查：

[~, flag] = chol(P, 'lower'); if flag > 0 P = P + 1e-10 * eye(n); % 加一个小对角阵恢复正定性 end

这个技巧虽然简单，但能救回无数次仿真崩溃。

4.4 算例设计与结果解读：怎么判断滤波器到底“好”还是“坏”

我建议用一个简单的单机无穷大系统做数值算例，设定机械功率在 (t=2s) 时阶跃增加5%，模拟负荷扰动。量测由真值叠加高斯白噪声生成：电压幅值噪声标准差0.01，注入功率噪声标准差0.02。滤波步长取0.01s，仿真时长10s，共1000步。

跑完之后看两个指标。一是状态变量的跟踪曲线，功角、转速、(E'_q) 的估计曲线应该紧跟真值，而不是紧追噪声；二是计算均方根误差（RMSE）：

[ \text{RMSE}x = \sqrt{\frac{1}{N}\sum{k=1}^{N}(x_k - \hat{x}_k)^2} ]

我在同样条件下跑EKF和UKF，UKF的功角RMSE通常比EKF小20%到30%，转速RMSE也类似。如果某个量测通道突然出现粗差，EKF的估计会被拉偏一大截，UKF由于Sigma点传播方式的平滑作用，抗粗差能力会略强一些，但依然需要在算法外层加坏数据检测。

另一个常被忽略的评估点是滤波的“一致性”，可以看归一化误差平方（NIS）。如果NIS长期远大于理论自由度，说明 (Q) 或者 (R) 设置不合理，滤波器对自身不确定性估计不准。这个指标写论文时很加分，调参时也很有用。

5. 常见问题与排查技巧实录：这些坑我真的都踩过

5.1 初学者高频问题速查表

这个表里的每一项背后都是一次惨痛调试。我记得有一次在UKF里连续发散，最后发现是初始 (P_0) 里转速项的方差给到了 100，相当于完全不确定初始转速，增益被放得巨大，系统直接失去了稳定性。

5.2 调试技巧：把滤波器拆开来一步一步验证

遇到“滤波发散”先别急着调参数，把预测和更新拆开看。

第一步，先看开环预测。把量测更新关掉，让状态只走动态方程，如果预测轨迹与真值偏差在几个步长内尚可接受，说明模型和离散化没问题；如果预测本身已经飞了，那不是滤波器的问题，是微分方程或RK4实现的问题。

第二步，单独检查量测方程。用已知状态算一遍 (h(x))，和对应量测的物理值比较。很多时候量测方程里的导纳矩阵转置问题、相角单位问题，都会导致更新量算错方向。

第三步，检查残差序列 (z_k - h(\hat{x}_{k|k-1}))。正常工作的滤波器残差应该近似零均值白噪声。如果残差有系统性的偏置，说明模型有偏差或者量测方程有错误；如果残差方差过大，说明R给小了或者Q给大了。

这个“从开环到闭环”的分步调试方法，比我当年一上来就盯着卡尔曼增益、协方差矩阵瞎猜高效得多。

5.3 工程级避坑经验：坏数据、相角参考和时间对齐

后面做实际PMU数据或者更接近工程的算例时，有几个问题是算法本身之外的，更容易让整个项目翻车。

坏数据问题。动态估计对量测中的粗差很敏感。我通常在外层加一个残差检测：当某通道的新息绝对值超过 (3\sigma) 时，临时把对应的 (R) 对角元调大（相当于降低这路量测的信任度），而不是直接丢弃通道，这样能避免滤波器状态跳变。这个方法虽然简单，但在工程算例里非常有效。

相角参考问题。PMU相角是相对于GPS秒脉冲的绝对相角，而发电机功角是相对于内部转子位置的物理量，两者之间存在一个旋转坐标系的对齐问题。不做坐标变换，直接把PMU相角当成量测值喂给滤波器，功角估计会出现持续漂移。

时间对齐问题。不同量测装置的上送延迟不一致就会导致滤波精度下降，甚至振荡。我在仿真数据中一般假设完全对齐，但用真实数据时必须做插值对齐。在Matlab里可以用interp1把不同通道的量测统一到同一时间网格上。

5.4 关于性能评估和论文写作的补充建议

如果这篇博文是给毕业设计或者小论文用的，建议把实验部分设计成三组对比：EKF、UKF、以及一个基准方案（比如稳定的静态WLS）。三组算法用同一套测试数据和同一个工况，统计RMSE和平均迭代时间。

这里有个实用的细节：在WLS基准方案中，PMU高频数据可以每隔几个PMU帧才解算一次断面；而EKF和UKF是逐帧递推，不需要迭代求解非线性优化问题。所以动态估计的真正优势不只是估计精度，还包括计算延迟的可控性。把这一点写清楚，论文的“动机”部分就站得住了。

另外，可以画一个“估计误差箱线图”或者“多工况RMSE对比表”，比单纯画三条时间曲线更有说服力。我还习惯加一个“参数敏感性分析”，把Q增大10倍、减小10倍前后的RMSE变化列成表，这表明算法对参数扰动有一定的鲁棒性，也是评审人很喜欢看到的内容。

写在最后，几句掏心窝的经验

这套代码我从头到尾实现、调试、改写了不下三遍，最大的体会是：电力系统动态状态估计的难点从来不在卡尔曼滤波公式本身，而在前面的系统建模和后面的参数调测。EKF和UKF的公式，一页纸就能写完；但把发电机动态方程写对、把量测方程映射对、把Q和R调到物理上说得通，才是真正花时间的地方。

最后分享一个个人习惯：我会在仿真脚本里把“预测状态曲线”、“只靠量测修正后的估计曲线”、“真值曲线”三条线画在同一张图上。这样一眼能看出滤波器是在“信任模型”还是“信任量测”，调参时非常有方向感。很多同学一看到发散就慌，其实只要把预测轨迹和量测轨迹分开看，问题基本都能定位到具体环节上。你先从单机模型把这个流程跑通，再往多机系统扩展，这套方法论是通用的。