电力系统动态状态估计：EKF与UKF算法的Matlab实现与对比

在电力系统的在线监测与运行控制里，动态状态估计一直是个绕不开的核心话题。用扩展卡尔曼滤波（EKF）和无迹卡尔曼滤波（UKF）去跟踪发电机功角、转速这些动态状态，是目前学术研究和工程尝试里最主流的做法。这套方案落地成Matlab代码以后，既能跑通仿真验证算法效果，又能为后续接入PMU实测数据打底。这篇博文围绕这套代码的实现思路、算法选型、建模细节和排坑经验展开，凡是正在做电力系统动态状态估计、或者准备在Matlab里实现卡尔曼系列算法的朋友，都能从中找到可以直接参考的路径。

1. 动态状态估计为何是电力系统的"实时眼"

1.1 静态状态估计的短板：数据与模型的滞后

传统电力系统状态估计以加权最小二乘为代表，它处理的是某一断面下的静态快照，利用SCADA系统采集的遥测、遥信数据，解算出母线电压幅值和相角。这套体系运行了几十年，在稳态工况下表现尚可，但面对现在的电网形态，短板越来越明显。

SCADA的数据刷新率通常在几秒到几十秒一次，量测到达时间不严格同步。当系统发生负荷快速波动、新能源出力突变、短路故障后的暂态过程，SCADA拿到的数据往往已经"过时"了。再加上传统状态估计模型以代数方程为主，不包含发电机转子运动这类动态方程，它给出的结果本质上是一个准稳态解，无法反映功角、转速在动态过程中的演化轨迹。换句话说，它在时间尺度上就达不到"动态追踪"的要求。

有人可能会说：既然有PMU（同步相量测量单元）了，不是可以直接测功角、转速吗？问题是PMU只能测电气量（电压、电流相量），功角和转速属于机械量，PMU并不能直接完全测量。虽然可以间接推导，但实际中量测噪声、线路模型误差、角度基准漂移都会让直接推导的结果不可靠。所以需要动态状态估计，把系统模型和PMU高精度量测融合起来，通过滤波算法实时估计内部动态状态。

1.2 动态状态估计的核心任务：在线追踪发电机的"运动轨迹"

动态状态估计的目标很明确：基于发电机转子运动方程、励磁绕组动态方程等状态方程，结合PMU提供的快速量测，在每一个时间步递归估计出发电机的动态状态。最具代表性的状态量包括转子功角δ、转速偏差ω、暂态电动势的d轴和q轴分量（对应发电机四阶模型）或更详细的状态。

这样做的工程价值很直接。调度和稳控系统需要知道当前发电机是否处于稳定运行范围，功角摆开有多大，转速偏差是否收敛，当系统受扰动后整个过程是否向稳定方向演化。如果这些状态能被实时、平滑、去噪地估计出来，后续的紧急控制、低频减载、失步预测就有了一组可靠的输入信号。也正因为如此，动态状态估计被普遍看作电力系统在线动态安全分析的数据前置环节。

要完成这个任务，滤波器选型就非常关键。卡尔曼滤波框架天然适合"状态方程预测+量测更新"的递归模式，而系统中的状态方程和量测方程都是非线性的，这就引出EKF和UKF两个经典选项。

2. EKF和UKF的原理与选型逻辑

2.1 EKF：对非线性做一阶线性化，代价是精度和稳定性

扩展卡尔曼滤波的思路非常直接：在每一个时间步，把非线性状态方程和量测方程在当前估计点附近做一阶泰勒展开，用雅可比矩阵代替线性卡尔曼滤波中的状态转移矩阵和量测矩阵，然后套用标准的预测-更新流程。

预测阶段的核心计算是：

• 状态预测：x_pred = f(x_est)，其中f是离散化的机电暂态状态方程；
• 协方差预测：P_pred = A * P_est * A' + Q，A是状态方程对状态的雅可比矩阵；
• 卡尔曼增益：K = P_pred * H' * (H * P_pred * H' + R)^(-1)，H是量测方程对状态的雅可比矩阵；
• 状态更新：x_est = x_pred + K * (z - h(x_pred))；
• 协方差更新：P_est = (I - K * H) * P_pred。

EKF看似简单，实战中却有几处暗坑。电力系统状态量的数量级跨层很大，功角以弧度计通常在零点几到几之间，转速偏差以标幺值或rad/s计，范围也在零点几量级，但不同状态对应的雅可比矩阵元素可能相差几十倍到上百倍。这会造成协方差矩阵的病态，影响数值稳定性。此外，EKF只是用一阶泰勒近似逼近非线性函数，系统非线性强（比如故障后的剧烈摆动阶段）时，线性化误差会被放大，滤波结果可能明显偏离真值甚至发散。

2.2 UKF：用一组sigma点绕开求导，逼近非线性传播

无迹卡尔曼滤波的核心思想是无迹变换。它不再对非线性函数做泰勒展开，而是选取一组带权值的sigma点，让这些点经过非线性函数传播后，再用加权统计方法重建均值与协方差。

具体说，如果状态维度是n，那么选取2n+1个sigma点，用尺度参数调整点在均值附近的散布程度，其中常见的参数配置包括α（决定sigma点的散布，通常取1e-3~1之间）、β（用于融入先验分布信息，高斯分布取2）、κ（通常取0或3-n）。这些sigma点通过状态方程和量测方程传播后，用它们计算预测均值和预测协方差，再相同框架下计算卡尔曼增益并更新状态。

UKF最直接的好处是不需要手推雅可比矩阵。对于电力系统这种状态方程和量测方程形式复杂、含三角函数和代数变量消除过程的场景，避免求导是节省人力和降低出错率的巨大优势。在非线性程度较高的暂态过程中，UKF一般能比EKF获得更高的估计精度，因为它对非线性函数的传播精度能达到三阶矩水平（对高斯分布而言），而不是像EKF那样只保留一阶。

代价也很明确，计算量大约是EKF的2n+1倍，因为每个时间步需要多次求值非线性函数。对n=4~9的发电机动态状态估计来说，这个计算量在Matlab仿真和现代硬件上完全可接受，所以UKF的性价比相当高。

2.3 两个放一起比较：什么时候选谁

下表是我在实际项目里总结的选型对照：

我的建议是：如果你刚开始做动态状态估计，先用电台模型把UKF实现跑通，因为代码链路短、不容易在求导环节卡住；如果论文或工程方案中同时需要对比方法，再把EKF补上，两者对比也能更直观地凸显算法差异。这套Matlab代码把两种滤波器都实现了，正好可以同台对比。

3. 系统建模：发电机动态模型与量测模型

3.1 用几阶发电机模型？如何离散化

动态状态估计的模型选择，精度和复杂度要平衡。经典二阶模型（摇摆方程）只描述转子运动，状态变量是功角δ和转速偏差ω，结构最简单，适合快速验证滤波器代码是否正确收敛。更贴合实际的是四阶模型，在δ、ω基础上增加d轴暂态电动势Ed'和q轴暂态电动势Eq'，描述励磁绕组和阻尼绕组的动态，能反映暂态过程中的电压变化。

以四阶模型为例，连续时间状态方程可写成：

dδ/dt = ω * ω_b - ω_s（不同文献量纲写法有差异，另一种常用写法是dδ/dt = ω - ω_s）

dω/dt = (Pm - Pe - D*(ω - ω_s)) / (2H)

dEq'/dt = (-Eq' + Ef + (Xd - Xd')*Id) / Tdo'

dEd'/dt = (-Ed' - (Xq - Xq')*Iq) / Tqo'

需要注意，这里电气量（Id、Iq、Pe）并非状态变量的直接函数，还要通过代数网络方程消去。要将这些连续方程变成离散状态空间形式，最简单的做法是采用一阶欧拉离散，它直观易实现，但步长小时精度才够。实际中我建议采用四阶Runge-Kutta离散，这样在PMU量测步长（通常20ms到100ms）下仍能保持足够的精度。

3.2 量测方程怎么构造

量测方程建模是这套代码里最容易出问题的地方。以机端PMU量测为例，假设量测量为机端电压幅值Vt、机端电压相角θ，或者还可以加发电机输出有功功率Pe和无功功率Qe。

量测方程的形式并不唯一，取决于你选择哪个坐标系和母线方程。如果机端电压相角θ与功角δ之间的角度差作为变量Y = θ - δ，那么量测方程就可以写成关于状态量和端电流的非线性方程：

Vt = sqrt(Vd^2 + Vq^2)

θ = δ + atan(Vq / Vd)（具体符号需根据参考轴定义调整）

这样量测方程天然是非线性的，而且量测矩阵H（EKF需要）推导起来比状态方程更繁琐，因为要经过潮流网络方程间接求导。这也是很多人在EKF实现中被卡住的环节，后面排查章节我会给出验证手段。

3.3 量测噪声、过程噪声和初值怎么定

卡尔曼滤波里Q和R矩阵的取值直接影响滤波性能。R矩阵可以从PMU的技术参数中估计，例如电压幅值误差0.1%~0.2%，相角误差0.01~0.02弧度，把它对角化设置为量测协方差矩阵即可。Q矩阵则复杂得多，它代表过程噪声，本质是对模型误差的容忍度。模型省略了励磁调节器、调速器动态，这些未知输入都会以过程噪声的形式体现。

初值的设置基于稳态潮流解：先跑一次潮流计算，得到发电机端电压和功角的初始稳态值，并以该稳态点计算其他状态初值。协方差P_est初始值通常设置为一个对角矩阵，数值不宜太小，若设置过小会让滤波器过于"自信"，导致初期的量测更新跟不上真值变化；太大又会让前几步估计波动明显。一般取各状态初值平方的一定比例即可。

4. Matlab代码实现：从框架到核心函数

4.1 工程化的模块划分

写Matlab代码不要一个脚本塞到底，我建议按模块划分：主脚本用于设置参数、生成仿真数据、调用滤波器并绘图；系统模型函数负责状态方程和量测方程；滤波器函数独立封装EKF和UKF。这样不同算例之间切换只需要修改参数，代码复用性高。

主脚本中需要定义开关量控制滤波器类型，例如用method='ekf'或method='ukf'选择调用哪个函数，最终统一输出估计状态序列，方便对比。

4.2 状态方程和量测方程的Matlab封装

状态方程函数建议写成统一形式：[x_next] = dst_func(x, u, dt, para)，其中u表示输入（如机械功率Pm、励磁电压Ef），para结构体放发电机和时间常数等参数。量测方程函数则写成[z] = meas_func(x, para)。两个函数需要保持输入输出接口固定，因为UKF的sigma点传播和EKF的预测阶段都要反复调用它们。

四阶Runge-Kutta离散实现的关键代码大致是：

function xn = f_rk4(x, u, dt, para) k1 = f_cont(x, u, para); k2 = f_cont(x + 0.5*dt*k1, u, para); k3 = f_cont(x + 0.5*dt*k2, u, para); k4 = f_cont(x + dt*k3, u, para); xn = x + (dt/6)*(k1 + 2*k2 + 2*k3 + k4); end

其中f_cont是连续时间状态方程，负责根据当前状态和输入计算导数。

4.3 EKF核心代码逻辑

EKF的关键点在于雅可比矩阵的计算。这里给出两种方式：一是解析法，用符号工具箱求出A和H；二是数值差分法，用有限差分近似。后者实现简单，但要注意步长选择，太大会引入截断误差，太小会淹没在浮点误差中。

核心更新代码如下：

function [x_est, P_est] = ekf_update(x_est, P_est, z, u, dt, para) % 预测 x_pred = f_rk4(x_est, u, dt, para); A = jacobian_state(x_est, u, dt, para); P_pred = A * P_est * A' + para.Q; % 量测预测 z_pred = meas_func(x_pred, para); H = jacobian_meas(x_pred, para); S = H * P_pred * H' + para.R; K = P_pred * H' / S; % 更新 x_est = x_pred + K * (z - z_pred); P_est = (eye(length(x_est)) - K * H) * P_pred; end

这里我特意加了数值差分辅助函数，当你对解析导数没有十足把握时，先用有限差分结果对比验证：

function A = jacobian_state(x, u, dt, para) n = length(x); A = zeros(n, n); f0 = f_rk4(x, u, dt, para); h = 1e-6; for i = 1:n x_pert = x; x_pert(i) = x_pert(i) + h; f_pert = f_rk4(x_pert, u, dt, para); A(:, i) = (f_pert - f0) / h; end end

4.4 UKF核心代码逻辑：sigma点生成与权重计算

UKF的实现关键在无迹变换和权重计算。按常用的比例修正形式，权重计算如下：

function [chi, Wm, Wc, c] = ut_transform(x, P, alpha, beta, kappa) n = numel(x); lambda = alpha^2 * (n + kappa) - n; c = n + lambda; % 计算协方差平方根 sqrtP = chol((n + lambda) * P, 'lower'); chi = zeros(n, 2*n + 1); chi(:, 1) = x; for i = 1:n chi(:, i+1) = x + sqrtP(:, i); chi(:, 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) = lambda / (n + lambda) + (1 - alpha^2 + beta); for i = 2:2*n + 1 Wm(i) = 1 / (2*(n + lambda)); Wc(i) = Wm(i); end end

注意这里用了chol函数求平方根，后续如果遇到"矩阵非正定导致chol失败"的问题，参见第5章排查办法。UKF预测和更新阶段的代码逻辑是把每个sigma点分别通过状态方程和量测方程传播，然后加权统计均值和协方差。

4.5 仿真数据生成：让"真值"有处可查

做滤波器最怕没有参照物。这套代码里，我把"真实系统"用精细的仿真模型生成：状态方程中不加过程噪声、量测方程输出后叠加高斯白噪声模拟PMU量测。这样既能得到noisy量测序列z，又能保留干净的状态真值x_true用于计算估计误差。

生成方式类似：

for k = 1:N x_true(:, k+1) = f_rk4(x_true(:, k), u(:, k), dt, para); z(:, k) = meas_func(x_true(:, k+1), para) + mvnrnd(zeros(1, m), R)'; end

然后以z作为滤波器的输入，比较x_est与x_true。常用的性能指标包括RMSE（均方根误差）和最大绝对误差，还可以记录单次滤波的总耗时。

4.6 参数整定与初始化示例

下面是我在单机无穷大算例中的一组典型参数，可供参考：

• 发电机惯性常数H = 3.5s，阻尼系数D = 2.0；
• 同步电抗Xd = 1.8，暂态电抗Xd' = 0.3，Xq = 1.7，Xq' = 0.45；
• 励磁绕组时间常数Tdo' = 7.5s，阻尼绕组时间常数Tqo' = 0.45s；
• 量测噪声：电压幅值标准差0.005，相角标准差0.01弧度；
• Q矩阵对角元取1e-4到1e-3量级，R矩阵按上述噪声标准差平方设置；
• 仿真时长5s，步长0.01s，PMU量测每0.02s采一次（50Hz）。

这套参数跑下来，我见过UKF的功角RMSE在0.001~0.01弧度量级，转速RMSE在1e-4量级，具体数值取决于故障场景强度。

5. 实测中遇到的问题与排查技巧实录

5.1 滤波器发散：最常见的翻车现场

动态状态估计里最让人头疼的就是滤波发散：前几步看起来正常，忽然某一步估计值突变，之后误差越来越大，直接把曲线跑飞。根据我调试的经历，原因通常归为几类。

一是Q矩阵设置不合理。Q太小意味着过程模型被过度信任，一旦模型与实际系统存在偏差，滤波器会不断把误差归结到量测上，输出出现震荡；Q太大会导致滤波过于依赖量测，失去平滑去噪能力。排查办法是先固定R，用一组对比仿真扫描Q的数量级，观察RMSE曲线变化。

二是初值严重偏离真值。在仿真中可以把初值从真值偏移10%~20%测试滤波器随时间收敛的能力。若发散，优先检查状态方程离散化是否稳定，尤其二阶模型里ω的量纲很容易写错。攻角对时间的导数到底是ω还是ω-1（标幺值基准下的转速偏差）必须要一致，否则状态演化方向都是错的。

三是协方差阵数值病态。电力系统状态量量纲跨度过大时，P矩阵条件数可能达到10^8甚至更高，EKF中尤甚。解决办法是归一化处理：功角用弧度，转速直接用rad/s的量纲，不要混用标幺值和其他量纲系统；或者对P做定期对称化和特征值裁剪，把极小负特征值修剪到0。

5.2 Jacobian矩阵算错：EKF的隐形炸弹

EKF对雅可比矩阵的准确性极其敏感。H矩阵算错一个符号，结果可能不是数值偏差，而是直接发散。最典型的错误集中在量测方程对功角的导数上，因为量测方程不仅含有状态变量的三角函数项，还要通过电流Id/Iq间接依赖状态，这层复合求导极易漏项。

我建议的验证方法是有限差分对照法：先写解析H矩阵，再用数值差分生成参考矩阵，两者对比误差在1e-4以内基本可以认定正确。这个验证可以做成一次性辅助脚本，不放进主循环，避免影响性能。还有个小技巧：当EKF和UKF在同样的条件下UKF正常而EKF发散，优先怀疑雅可比矩阵而不是滤波器结构，因为UKF不依赖雅可比。

5.3 UKF的协方差非正定问题

UKF虽然避开了雅可比求导，但引入了另一个典型问题：协方差矩阵P在递推中可能失去正定性，导致chol分解失败，报错信息通常是"Matrix must be positive definite"。

这类问题多出现在初始P矩阵设置不当、量测噪声过小导致S矩阵接近奇异、或者步骤中浮点误差累积。应对方案：

• 在每次更新后强制对称化：P = (P + P') / 2；
• 诊断矩阵最小特征值，若接近0或为负，加一个很小的正对角阵（比如1e-8 * eye(n)）作为正则项；
• 改用平方根UKF（SR-UKF），它直接递推协方差的Cholesky因子，从根上避免非正定问题，代价是实现更复杂。
• 检查参数kapppa，当n较大时取kappa=3-n会导致lambda为负，传播中可能出现负权值，这也是数值不稳定的来源之一。建议使用GTF（Gaussian-to-First）参数组合，即alpha=1，beta=0，kappa=3-n，或者经典alpha=1e-3，beta=2，kappa=0的组合，并确认lambda满足必要的条件。

5.4 离散化步长和采样周期怎么选

PMU上送速率多为10~50帧每秒，而发电机暂态过程的时间常数从毫秒级（次暂态分量）到秒级（功角摆动）都有。离散化步长要能覆盖快变过程，否则系统微分方程的数值解本身就发散。

我的经验是：仿真/状态预测的积分步长取0.001~0.01s，量测更新步长取0.02~0.1s，两者之间可以不等距——积分步长小时，量测更新之间做多个积分小步，等更新时刻到了，再用当时的量测做校正。这样滤波器结构更贴近实际，也利于对比不同PMU速率对估计精度的影响。感兴趣的话还可以做一个"量测丢失率"测试，模拟PMU丢帧场景，看看滤波器在缺少某些时刻量测时是否能靠模型预测继续维持基本可用。这也是动态状态估计走向实际应用时不得不面对的问题。

6. 算例对比：EKF与UKF的实测表现

6.1 同一场景下的性能对比结果

以某单机无穷大系统为算例，在1s时设置一个三相短路故障、1.15s切除的暂态场景，分别用EKF和UKF做动态状态估计，我测试得到的结果具有明显的代表性。

EKF在功角摆动幅度较大、系统非线性较强的时段出现了明显的估计偏差，最大功角估计误差约0.02弧度，转速误差约2e-3 rad/s。而UKF在整个过程中保持了更平滑的估计轨迹，功角最大误差约0.005弧度，转速误差约5e-4 rad/s，精度高了一个量级左右。计算耗时方面，n=4时UKF每步需要9次函数求值，耗时约为EKF的3~5倍（因为EKF每次迭代还需要计算雅可比，所以倍数低于理论预期的2n+1）。

当然，在系统非线性很弱、接近稳态运行时，EKF与UKF的估计结果几乎重合，这时候UEKF在计算速度上的优势就体现出来了。这就印证了选型阶段的结论：强非线性暂态过程选UKF，轻载稳定场景EKF性价比更高。

6.2 面对量测噪声差异时的鲁棒性

我还做过一组量测噪声敏感度测试，把量测噪声水平从0.005逐步加大到0.05，观察两种滤波器的RMSE变化。结果是UKF的性能下降更缓慢，在量测噪声较大时依然能维持可用的估计质量；EKF则对噪声的容忍度明显更低，噪声大时滤波轨迹出现较明显的跟随滞后。这背后的原理是UKF在更新步骤中用的协方差传播更准确，对量测异常值的反应更平缓。

6.3 从单机走向多机系统时要注意什么

当从单机系统扩展到多机系统（比如IEEE 9节点、39节点），状态维度上升，状态方程和量测方程之间的耦合变复杂，量测方程需要包含互联母线的电压相量信息。此时有两个建议：一是状态可分割为每台机一个局部滤波器（分散式动态状态估计），局部量测覆盖本机及相邻母线，降低全局维度，适合并行计算；二是如果坚持集中式全状态估计，务必关注P矩阵维度增大后的数值稳定性，Q矩阵也要相应调整，因为不同发电机之间过程噪声相关性需要建模，否则多机系统背景下滤波器容易误发散。

7. 扩展方向：这套代码还能往哪走

做滤波算法和仿真代码，最忌讳的就是"跑通一次就封存"。这套EKF和UKF代码框架其实可以平滑迁移到更多场景。

叶片尺度上，可以替换发电机模型，例如加入励磁调节器（AVR）和调速器（Governor）的简单动态模型，状态维度扩展到8~10维，观察动态状态估计在闭环控制影响下的表现。再进一步，把光伏、储能逆变器的动态模型纳入估计范围，研究新型电力系统下多时间尺度的状态估计问题，这在当前碳中和背景下很有研究价值。

算法层面，可以对比更先进的滤波方案：中心差分卡尔曼滤波器（CDKF）、平方根UKF、容积卡尔曼滤波（CKF），还有粒子滤波和H∞滤波，这些方法在处理更强的非线性和非高斯噪声时各有擅长。动态状态估计结合深度学习也是热门，比如用神经网络自动调整Q、R矩阵，或者用RNN学习模型误差残差作为过程噪声补偿，都能在现有Matlab框架上扩展。

把PMU实测数据接入这套代码也不难，只需把量测数据来源从仿真生成换成数据文件读取，注意时间对齐和标幺值换算即可。到这一步，代码就不再只是算法验证工具，而是具备工程原型潜质的应用底座了。

回头再看，整个动态状态估计项目最核心的心得其实只有两条：第一，建模一致性比算法选择更重要，状态方程和量测方程一旦在量纲、参考系、离散化上出错，什么滤波器都救不回来；第二，Q和R的整定不是一次性的，配合不同场景反复调参，才能真正体会每种滤波器的性格。这套Matlab代码作为起点，足够你在EKF和UKF之间来回切换、设置各种故障场景并评估性能了。踩过那些发散的坑、调过参数之后，你对动态状态估计的理解会扎实很多。