从EKF到UKF：Matlab实现电力系统动态状态估计全流程解析

做电力系统动态状态估计这些年，EKF和UKF是我最常用的两个非线性滤波工具。今天这篇就来记录一下我用Matlab从零实现这两种滤波器，并在IEEE标准节点系统上跑通全过程的思路、代码和踩坑记录。内容偏实操，我会把模型怎么建、雅可比怎么求、Sigma点怎么采、参数怎么调一次性讲清楚，适合刚接触动态状态估计的研究生，也想给正在做PMU数据融合的工程师一点参考。

先说结论：EKF胜在计算量小、实现直观，但强非线性场景容易因线性化误差发散；UKF不用求导，精度和鲁棒性都要好一截，代价是多算2n+1个Sigma点。真要在Matlab里跑，两者代码量其实差不多，真正花时间的是模型搭建和参数整定。

1. 为什么要做电力系统动态状态估计

1.1 静态估计的局限性，和动态估计的必然性

传统调度中心用的状态估计，本质上是加权最小二乘(WLS)对一段时间的断面做拟合，依赖SCADA系统慢速周期的数据，通常几分钟才出一个结果。这个时间尺度应对系统稳定运行可以，但遇到负荷快速波动、新能源出力随机变化或者故障后的暂态过程，几分钟前的断面代表不了当前状态。更关键的是，紧急控制决策需要毫秒级到秒级的状态感知，SCADA数据根本喂不过来。

所以有了基于PMU的动态状态估计。PMU以20到100帧每秒的速率同步采集电压电流相量，配合发电机动态方程，就能在线递推地跟踪功角、转速这些状态量。这个“动态”指的不是简单的时变参数，而是把状态估计嵌入到状态方程里，用上一时刻的估计值预测下一时刻，再用实时量测修正，形成闭环递推。

这种思路其实就是卡尔曼滤波的经典框架。电力系统模型是非线性的，所以标准卡尔曼滤波不能用，得用非线性版本。EKF和UKF就是我在实际项目里用到最多的两个非线性卡尔曼变体。搞清楚这两个，就能解决电力系统中绝大多数动态状态估计问题，从单机无穷大系统到多机系统，算法框架是通用的。

1.2 EKF和UKF：非线性卡尔曼滤波的两个主力

EKF的思路很直接：非线性函数在估计点附近做一阶泰勒展开，用雅可比矩阵代替原本的线性状态矩阵。好处是只要会求导，就能用，计算效率高。坏处也明显：线性化误差大，尤其当系统强非线性或者初始误差大的时候，协方差传播不准，容易出现滤波发散。另外求雅可比矩阵容易出错，尤其是多机电力系统里相互耦合的功率方程。

UKF走的是另一条路：用一组确定的Sigma点去近似状态分布，这些点经过非线性函数传播后，再按权重合成新的均值和协方差。整个过程完全不涉及求导，非线性函数就是原样用的。所以UKF的处理精度在二阶以上，比EKF高一阶，对强非线性系统的适应能力也强。

我在做IEEE 9节点系统时，两个都实现了。要说哪个更好用，我的体会是UKF更省心——不需要手推雅可比，参数设置就alpha、beta、kappa三个数，整定好就不太容易发散。但EKF也有存在价值：状态维数比较高的时候，UKF的Sigma点数量会膨胀，计算量增长快，而EKF计算量相对可控。所以实际选型往往是在精度和实时性之间做权衡，不是无脑选UKF。

2. 从原理到代码：EKF和UKF的核心推导

2.1 状态空间模型与离散化

在做动态状态估计之前，必须先确定被控对象的状态方程和量测方程。我用的模型是发电机二阶经典模型，虽然简单，但足以说清楚整个算法流程，适合做教学和基准对比。状态量取功角δ和角速度ω（标幺值形式），输入是机械功率P_m，输出是电磁功率P_e。

连续时间模型如下：

dδ/dt = ω - ω_s dω/dt = (P_m - P_e - D*(ω - ω_s)) / (2H)

其中D是阻尼系数，H是惯性时间常数，ω_s是同步转速。P_e是节点注入功率，和系统其他节点电压、相角非线性相关。对于单机系统可以简化成P_e = E'V/X*sin(δ)，E'是暂态电动势，X是电抗。这个式子是非线性的，正好用来验证EKF和UKF。

离散化我用的是一阶欧拉法，步长dt取0.01秒，对应PMU采样频率100Hz。离散后的状态方程：

delta(k+1) = delta(k) + (omega(k) - omega_s)*dt omega(k+1) = omega(k) + (P_m - P_e(k) - D*(omega(k)-omega_s))/(2H)*dt

因为dt较小，欧拉法精度够用。如果仿真步长更大，可能要改用梯形法或Runge-Kutta。量测方程就直接取P_e(k)和V(k)（如果可观测电压幅值）。这里要特别提醒：状态量和量测量的单位要统一，用标幺值就别混入有名值，不然卡尔曼增益会算得一团糟。

2.2 EKF的线性化：雅可比矩阵从哪来

EKF的核心在于两个雅可比矩阵：状态转移矩阵F和量测矩阵H。对于上面的离散方程，F是状态方程对状态向量的偏导数。因为状态方程里omega(k+1)含有P_e(k)，而P_e是δ的函数，所以F不是简单的常数矩阵，需要逐项求偏导。

以单机系统为例，设状态x=[delta; omega]，量测z=P_e，雅可比矩阵H就是：

H = [dP_e/ddelta, 0]

其中dP_e/ddelta = E'V/X*cos(delta)。如果量测还包含电压幅值V，那V对delta、omega也有偏导，这里不展开。

实际写Matlab代码时，我建议先用符号运算做验证，再手写解析式。因为手写偏导容易漏项，尤其多机系统里功率对相角的偏导直接用B矩阵和G矩阵表示时更容易写错。我踩过的坑就是在推导H矩阵时把B和G的位置搞反，导致滤波一会就发散。

EKF主循环的逻辑很简单：先根据当前状态和协方差做预测，再用量测更新。预测部分用状态方程递推，协方差用F、Q传播；更新部分用H算卡尔曼增益，然后修正状态和协方差。代码结构后面会给出。

2.3 UKF的Sigma点采样：不求导也能传播不确定性

UKF的思想是无迹变换。假设状态x的均值为x_mean，协方差为P，那么构造2n+1个Sigma点，其中n是状态维度，常见生成方式：

X_0 = x_mean X_i = x_mean + (sqrt((n+lambda)*P))_i , i=1..n X_i = x_mean - (sqrt((n+lambda)*P))_i , i=n+1..2n

这里的lambda = alpha^2*(n+kappa) - n，alpha控制Sigma点的散布程度，kappa是次级缩放参数。sqrt((n+lambda)*P)表示对矩阵做Cholesky分解得到的下三角矩阵的每一列。这组Sigma点通过非线性函数传播后，按权重计算新的均值和协方差。

权重定义：

W_m_0 = lambda/(n+lambda) W_c_0 = lambda/(n+lambda) + (1 - alpha^2 + beta) W_m_i = W_c_i = 1/(2*(n+lambda)), i=1..2n

beta是用来包含状态分布先验信息的参数，对高斯分布取2是最优。UKF代码比EKF长一些，但胜在不需要求雅可比，直接把你写好的状态方程和量测方程当黑盒用就行。这也意味着，如果以后换了更复杂的发电机模型，只要替换函数句柄，滤波框架完全不用动。

3. Matlab实现的关键环节与实操步骤

3.1 仿真数据生成：用MATPOWER搭IEEE 9节点系统

真试验证算法，我推荐用MATPOWER导入标准节点系统。IEEE 9节点系统是最经典的三机九节点，数据量不大，但足够体现多机交互的动态。MATPOWER里提供case9.m文件，直接加载就能得到节点、线路、发电机的稳态参数。先做一次潮流计算，得到系统稳态工作点，作为动态模型的初值。

接下来生成“真值轨迹”：我在稳态工作点基础上设置一个扰动，比如在某个节点加一个持续0.1秒的三相短路故障，然后用数值积分（ode15s或自己写的R-K）跑一段时间，得到功角、转速、节点功率的演化轨迹。这个轨迹就是真实状态，用来生成量测数据和控制滤波效果评估。

生成量测时，我分别在功率量测和功角量测上叠加高斯白噪声。PMU的幅值测量误差可以按0.01到0.02标幺值算，相角误差按0.01弧度算，转换成协方差就能得到量测噪声矩阵R。不能用一个凭空的R，最好结合PMU的出厂指标和实际测试结果来定。

3.2 EKF代码实现与Q/R参数整定

EKF代码不难写，但有几个细节决定成败。先说初始化：状态初值在真实初值附近加一个随机偏差，模拟实际情况下只知道大概状态。协方差P0设为一个对角阵，对角线可以取状态误差的平方。比如初值偏差0.05弧度，P0对角元素就取0.0025。Q矩阵表示过程噪声协方差，代表模型误差和扰动的不确定性，一般也设成对角阵。

我常用的Q整定经验是：先用小量级（比如1e-6到1e-4）跑，看滤波轨迹是否平滑；如果估计结果过于依赖量测，波动大，就适当增大Q；如果估计结果跟不上真实轨迹变化，说明Q太小，模型预测过于自信。这是一种“黑盒整定”的笨办法，但在没有精确噪声统计的时候很实用。R则根据传感器精度来，如果PMU标称误差是1%，R对角元素就取(0.01)^2=1e-4。

EKF主循环的Matlab伪代码如下：

for k = 2:N % 预测 x_pred = f(x_est(:,k-1), u, dt); F = compute_F(x_pred); P_pred = F * P_est(:,:,k-1) * F' + Q; % 更新 H = compute_H(x_pred); K = P_pred * H' / (H * P_pred * H' + R); x_est(:,k) = x_pred + K * (z(:,k) - h(x_pred)); P_est(:,:,k) = (eye(n) - K * H) * P_pred; end

compute_F和compute_H是手写的雅可比函数，h是量测方程。注意每一步H要在预测点x_pred处计算，不是在上一时刻的估计点算，这是新手容易忽略的。

3.3 UKF代码实现与参数调优

UKF的实现可以抽象成几步：给定当前均值和协方差，生成Sigma点，然后分别通过状态方程和量测方程传播，再计算增益。Matlab里可以用函数句柄传入状态方程，不用像EKF那样维护雅可比，代码复用性更好。

我贴上核心部分的伪代码：

for k = 2:N % 生成Sigma点 [X_sig, Wm, Wc] = generateSigmaPoints(x_est(:,k-1), P_est(:,:,k-1), alpha, beta, kappa); % 状态传播 X_sig_pred = f(X_sig, u, dt); x_pred = X_sig_pred * Wm; P_pred = (X_sig_pred - x_pred) * diag(Wc) * (X_sig_pred - x_pred)' + Q; % 量测传播 Z_sig = h(X_sig_pred); z_pred = Z_sig * Wm; Pzz = (Z_sig - z_pred) * diag(Wc) * (Z_sig - z_pred)' + R; Pxz = (X_sig_pred - x_pred) * diag(Wc) * (Z_sig - z_pred)'; % 更新 K = Pxz / Pzz; x_est(:,k) = x_pred + K * (z(:,k) - z_pred); P_est(:,:,k) = P_pred - K * Pzz * K'; end

generateSigmaPoints里要对P做Cholesky分解，所以P必须保持正定。前面提到，alpha取1e-3时Sigma点离均值很近，对于状态变化剧烈的场景反而容易低估不确定性；我实际跑下来的经验是alpha取0.01到0.1之间比较稳，kappa取0，beta取2。如果Cholesky分解报错，说明P矩阵不是正定，可以在分解前加一个很小的对角阵，比如1e-12*eye(n)，手动保证数值稳定。

4. 结果对比与性能分析

4.1 精度对比：RMSE指标怎么算

评价滤波精度，我一般用均方根误差(RMSE)来对比，分别统计状态量（功角和转速）以及输出功率的估计误差。公式不复杂：RMSE = sqrt(mean((x_est - x_true).^2))。注意要排除前几个点的暂态，我通常让滤波器跑0.5秒后再开始统计，给滤波器一个收敛过程。

用IEEE 9节点系统的三台发电机状态做对比，在相同初值偏差、噪声和故障扰动下，UKF的功角RMSE通常比EKF小30%到50%。原因还是EKF的线性化误差，尤其在系统经历大的扰动、功角摆动幅度比较大的时段，EKF的协方差传播精度不足，修正效果打折扣。

下面是我一次典型仿真中得到的数值（具体数值因扰动设置略有差异，但趋势一致）：

UKF在非线性更强的时段优势会更明显；系统平稳时两者差别不大。

4.2 计算效率和实时性对比

精度之外，计算效率也是选型的重要考量。我在同一台机器、相同步数和数据长度下分别跑了100次，记录单步平均耗时。UKF生成Sigma点、传播每个点，计算量大约是EKF的1.5到2倍。对于9节点系统，状态维度n=6（三台发电机，每台2个状态），UKF需要生成13个Sigma点，还不算太夸张。但如果把模型换成四阶甚至六阶，状态维度上到20以上，UKF的计算负担会明显上升。

如果你要做实时在线估计，且系统规模不大，UKF完全可以跑在实时仿真器上。如果状态维度很大，且模型有比较平滑的非线性特性，EKF的计算优势就体现出来了。实际项目中也可以考虑“混合策略”：系统运行平稳时用EKF省资源，检测到大扰动时切到UKF提高精度，这个思路后续可以扩展。

4.3 噪声协方差Q/R的敏感性分析

Q和R的匹配程度直接影响滤波质量。我把Q设成固定对角阵，发现当Q过小时，滤波器对新量测的信任度很低，估计轨迹会“滞后”于真值；当Q过大时，状态被噪声污染，估计曲线毛刺明显。R的设定错误也类似：R过大，滤波器过度相信模型，无法跟踪突变；R过小，量测噪声会被放大，估计值波动剧烈。

建议在调试阶段先用一组名义Q/R跑出基准，然后固定其中一个参数，把另一个乘以系数0.1到10扫描一遍，看RMSE变化。这个过程可以用脚本批量跑，画出RMSE随参数变化的曲线，帮助理解系统的灵敏方向。实测下来，EKF对Q/R比UKF更敏感，因为EKF的协方差传播不准，如果R设小了，卡尔曼增益会异常偏高，极易发散；UKF在这方面稳定得多，这也是我推荐新手先用UKF起步的原因。

5. 常见问题与排查技巧

5.1 滤波器发散，先查这几个地方

我在调试EKF时遇到过最典型的问题就是滤波发散：估计值飞掉，曲线直接冲出坐标系。排查顺序很有讲究：先看状态方程对不对，再看量测方程对不对，最后才怀疑参数问题。

具体检查点包括：初值是否落在可行域内（δ不可能超出一堆电抗决定的稳定边界）；P0是否反映了真实的不确定性，如果P0设得过大，会让早期增益过大；Q是否过小，导致模型预测误差被低估；R是否过小，导致过度相信量测；雅可比是否计算正确，这个我专门用符号工具箱验证过。

如果一时查不出问题，可以在滤波前几步打印新息序列（z - h(x_pred)）的均值和协方差。新息均值不该持续偏大，新息协方差应该和理论值大致匹配。如果新息序列呈强相关性或偏差很大，多半是模型与量测不匹配。

5.2 雅可比矩阵写错，用符号工具自查

EKF的雅可比矩阵是手写最容易出错的地方。多机系统的状态方程里，P_e对每个发电机功角都有偏导，交叉项特别容易漏。我的做法是先用Matlab Symbolic Toolbox把f和h的解析式写出来，用jacobian函数求偏导，然后对比手写的函数数值结果。比如随机撒10个状态点，比较手写雅可比和符号雅可比的计算结果，看到所有点误差小于1e-8再放心用。

还有一个细节：雅可比是离散化后的状态方程对状态的偏导，不是连续微分方程的偏导。如果你直接用连续模型雅可比套到离散递推里，步长一大会有明显的截断误差。实际调试中我发现，步长0.01s时连续雅可比还能用，但步长大于0.05s时误差就比较明显了，最好自己对离散方程重新推导雅可比，或者直接用数值逼近的雅可比（但那会增加计算量）。

5.3 UKF三个参数怎么设最省心

UKF参数里alpha、beta、kappa，网上有各种说法。我结合自己的测试给出一个推荐起点：alpha=0.01，beta=2，kappa=0。这个组合在大多数电力系统动态估计场景下都是个不错的起点。alpha越小，Sigma点越集中，对于接近线性的区域精度高，但对于强非线性可能低估不确定性；alpha在0.01到0.1之间是一个折中区间。beta=2是所有高斯分布场景下的理论最优，没什么好调的。kappa取0保证半正定（为3-n可能是另一个常用选择，但很多文献实际测试对单峰高斯影响很小）。

如果Cholesky分解报错，我的经验是把alpha调大一点，或者给P加上一个很小的对角扰动，不要轻易调kappa到负值，因为负的kappa可能会导致权重为负，进而协方差失去正定性。

5.4 步长与采样周期不匹配

动态状态估计对时间离散化很敏感。我的仿真里用dt=0.01s做状态递推，量测假设每0.02s来一个（对应50Hz采样），所以滤波器每两步更新一次。这样做的原因是模拟PMU实际采样间隔可能大于积分步长。如果直接把观测更新步长设成和状态递推步长一样，结果看起来没问题，但和真实场景有差距。

另一种情况是PMU采样率很高，量测每0.01s来一次，但状态递推步长也需要0.01s，这样倒简单。最怕的是把状态递推步长设得很大，比如0.1s，那么EKF和UKF的离散化误差都会显著增大，估计性能变差。我的经验是：状态递推步长至少要比系统主要动态时间常数小一个数量级，电力系统机电振荡周期大约在0.5到2s之间，所以0.01s到0.02s的步长是合适的。

如果模型包含更快的电磁暂态，那就要用更小的步长，或者把状态增广，这已经超出常规动态状态估计的范围了。如果读者做的是这部分，可以留言交流，我对电磁暂态下卡尔曼滤波也做过一些尝试，但那是另一个话题了。

最后再分享一个实用习惯。我写Matlab代码时，会把EKF和UKF封装成两个通用函数，输入是状态方程句柄、量测方程句柄、初始状态、协方差以及Q/R，输出是状态估计序列。这样后续想换模型、换节点系统，只需要改函数句柄和初始参数，滤波器主体完全不用动。配合MATPOWER批量生成场景，可以快速做多工况对比。这一套下来，不仅省了重复写代码的时间，也让我更专注于模型本身的问题。