企业官网建设流程全解析

项目标题: 基于扩展(EKF)和无迹卡尔曼滤波(UKF)的电力系统动态状态估计(Matlab代码实现)摘要描述: 暂无关键词: 暂无相关热搜词: EKF, UKF, 电力系统动态状态估计, Matlab

搞电力系统状态估计这块的同行应该都有感触：传统的静态状态估计（加权最小二乘法那套）算一个断面还行，但真要拿它去跟踪功角摇摆、母线电压的连续波动，基本是力不从心。这几年新能源大规模接入，系统动态特性比过去复杂了不止一个量级，我手里的项目里需要实时跟踪动态过程的场景越来越多，所以动态状态估计（Dynamic State Estimation, DSE）几乎成了绕不开的工具。

而在DSE的实现方案里，扩展卡尔曼滤波（EKF）和无迹卡尔曼滤波（UKF）又是最主流的两条技术路线。说白了，EKF的思路就是把非线性模型在当前工作点做一阶泰勒展开，硬凑成线性系统再用标准卡尔曼滤波的框架；UKF则是利用无迹变换，用一组精心挑选的Sigma点去逼近状态分布，不需要求导，也不会把高阶非线性信息直接扔掉。

这篇文章我不打算给你堆公式堆到头晕，而是想从"为什么要用这俩方法"、"电力系统模型怎么搭"、"Matlab代码到底怎么写"以及"实际调参会遇到哪些坑"这几个角度，把这条技术路线完完整整地拆开讲一遍。我会把我在实际项目里跑通的代码逻辑、参数设置和踩过的坑都放出来，适合正在做毕业设计、或者刚接手动态状态估计相关课题的同学作为参考。

1. 为什么"动态"状态估计绕不开EKF和UKF

1.1 静态估计的短板：只能看照片，不能看视频

传统电力系统状态估计，绝大多数是用加权最小二乘法（WLS）去解一个静态非线性优化问题。输入是SCADA系统采集的遥测数据——有功、无功、电压幅值这些，输出是系统当前断面的状态量（母线电压幅值和相角）。这个方法非常成熟，但在原理上有一个硬伤：它默认系统是稳态的，估计结果只是某个时刻的"快照"。

如果电网负荷平稳、拓扑不变，WLS的表现是够用的。但问题是，现在的电网真不是这个样子。风电和光伏出力随机波动，电动汽车充电负荷忽高忽低，还有各种电力电子设备导致的快速暂态过程。这时候系统状态每时每刻都在变化，你拿一个断面去当"真相"，下一秒钟这个值就已经过时了。就像拍照只能得到静态照片，而你要的是看清一段连续运动的视频。

动态状态估计正是为了补上这个缺口。它的本质是递推贝叶斯估计，利用系统的动态模型（比如发电机转子运动方程）做一步预测，再用新的量测值去修正预测。这个过程跟目标跟踪里的卡尔曼滤波思想完全同源，所以也叫"基于模型的实时状态追踪"。

1.2 标准卡尔曼滤波在电力系统里为什么不好使

如果电力系统动态模型是线性的，那直接用标准卡尔曼滤波（KF）就行了，公式简单、计算量小、还有最优性保证。但电力系统的动态模型，说句实话，几乎没有线性的。

举个最常见的例子，发电机采用二阶经典模型时，转子运动方程里就有功角和电磁功率的正弦函数关系：

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

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

其中电磁功率(P_e)通常写作：

[ P_e = \frac{E'V}{X}\sin\delta ]

好嘛，这方程里带了个(\sin\delta)，整个系统就是个强非线性模型。量测方程同样不省心——如果量测里有功率量测，那就涉及电压幅值与相角的乘积加三角运算，非线性程度只高不低。

标准的卡尔曼滤波要求状态方程和量测方程都是线性的，这一条电力系统天生满足不了。所以实际工程里，大家就从两个方向扩展卡尔曼滤波：一个用泰勒展开做局部线性化，就是EKF；另一个用采样点去数值逼近非线性传播，就是UKF。这两个方法也是目前电力系统动态状态估计里落地最广的算法。

1.3 EKF和UKF的核心思路差异

EKF的套路可以概括为"线性化后套KF"。在每一步滤波时，把非线性函数在当前估计值附近做一阶泰勒展开，忽略高阶项，得到近似的线性模型，然后完全套用标准卡尔曼滤波的五条公式。这样做的好处是思路直白、实现简单，坏处也显而易见：如果系统非线性很强，这个截断误差会被卡尔曼增益放大，导致估计精度下降甚至发散。

UKF则换了个思路：与其把非线性函数线性化，不如对状态变量的分布做文章。它用无迹变换，在原状态分布中选取一组Sigma点（一般取(2n+1)个），将这组点分别通过非线性函数传播，再从变换后的点中计算均值和协方差。这个过程不需要计算雅可比矩阵，而且理论上能精确到二阶，对强非线性系统的适应能力明显强于EKF。

我个人的选型建议是：模型非线性程度不高、算力敏感、希望代码短平快，那就EKF；系统强非线性、量测精度较高、对估计精度要求严格，那就UKF。两条路都不难，后面我会把代码实现的关键差异也讲清楚。

2. 电力系统动态模型搭建：状态方程与量测方程怎么选

2.1 状态变量的选取原则

做动态状态估计，第一步是确定"你要估计什么"。电力系统里，最常用的动态状态变量来自同步发电机的转子运动方程和电磁暂态方程。我在这篇文章里采用最常见的三阶实用模型，状态变量取：

• 功角(\delta)：发电机转子相对同步旋转坐标系的角位移
• 角速度偏差(\Delta\omega)：转子电角速度与同步速度的偏差
• q轴暂态电动势(E'_q)：反映励磁绕组磁链变化

这里不取更复杂的四阶、五阶模型，一方面是为了把EKF和UKF的算法逻辑讲清楚，避免被模型细节淹没；另一方面，三阶模型在暂态稳定分析里已经能反映功角摇摆和励磁动态的基本特征，用于验证动态状态估计算法完全够用。

状态变量列成向量就是：

x = [δ_1, Δω_1, E'q_1, δ_2, Δω_2, E'q_2, ...]

每个发电机对应三个状态量，如果是对IEEE 14节点系统里那几台发电机做估计，状态向量的维数需要根据发电机数量来定。

2.2 连续状态方程的离散化处理

三阶模型中，连续时间的状态方程大致是：

[ \frac{d\delta}{dt} = \omega_0 \Delta\omega ]

[ \frac{d\Delta\omega}{dt} = \frac{1}{M} \left( P_m - \frac{E'_q V}{X'_d}\sin\delta - D\Delta\omega \right) ]

[ \frac{dE'_q}{dt} = \frac{1}{T'd}\left( E{fd} - E'_q - (X_d - X'_d) I_d \right) ]

写代码时不能直接处理连续方程，必须离散化。工程上做EKF/UKF，常用最简单的一阶欧拉离散，或者四阶龙格库塔（RK4）。采样时间取0.01秒到0.02秒就可以，跟PMU的典型采样率（50帧/秒）也对应得上。如果采样时间太大，离散化误差会明显上升；太小则计算负担变大，实时性受影响。

我在项目里实际用的是RK4离散化EKF，虽然单步计算量比欧拉大一点，但在(T_s=0.01s)时数值稳定性明显更好，尤其适合后续要扩展做非线性程度更高的模型。

2.3 量测方程：用PMU还是SCADA

现在做动态状态估计，量测基本都用同步相量测量单元（PMU）的数据。PMU可以高频率（几十到上百Hz）提供带时标的电压相量、电流相量，直接给出相角信息——这对动态估计来说是极其宝贵的。

量测向量我一般取：

• 发电机机端电压幅值(V_t)
• 发电机机端电压相角(\theta_t)
• 有功功率(P_e)
• 无功功率(Q_e)

这样量测方程(z = h(x))写出来，就是根据状态变量反推这几个量。以有功为例，量测方程里同样是带(\sin\delta)的非线性表达式。这一步是EKF雅可比矩阵的核心来源，也是容易写错的地方。

2.4 过程噪声和量测噪声怎么定

卡尔曼类滤波器的性能，很大程度上取决于噪声矩阵Q和R的设置。

• (Q)是过程噪声协方差矩阵，它反映的是模型本身的不确定度。模型越不准，Q的元素就应该越大，让滤波器更信任量测。
• (R)是量测噪声协方差矩阵，它反映的是表计/PMU的测量误差。R越大，表示量测越不可信，滤波器会更依赖模型预测。

这两个矩阵如果设置得太离谱，EKF和UKF都会很快发散。后面我会专门讲我调Q、R的经验，这里先记住一个原则：宁可让Q稍微偏大，也不能让R偏大。因为模型误差是客观存在的，被模型带着走通常比被量测噪声带着走更危险。

3. Matlab代码实现：EKF的雅可比矩阵与UKF的Sigma点

3.1 EKF的五个核心公式和代码骨架

EKF的算法流程，一句话概括就是"预测—修正"。每个采样周期执行两步：

预测步：

• 用当前状态估计值代入状态方程，得到一步预测(\hat{x}_{k|k-1})
• 用雅可比矩阵(F_k)更新状态协方差：(P_{k|k-1} = F_k P_{k-1} F_k^T + Q)

滤波步：

• 计算量测预测(\hat{z}k = h(\hat{x}{k|k-1}))
• 求量测雅可比矩阵(H_k)
• 算卡尔曼增益：(K_k = P_{k|k-1} H_k^T (H_k P_{k|k-1} H_k^T + R)^{-1})
• 更新状态：(\hat{x}k = \hat{x}{k|k-1} + K_k (z_k - \hat{z}_k))
• 更新协方差：(P_k = (I - K_k H_k) P_{k|k-1})

Matlab代码骨架大致是这样：

% 状态转移函数和量测函数的句柄 f_func = @(x) system_dynamics(x, u, Ts); h_func = @(x) measurement_function(x); % EKF主循环 for k = 1:N % 预测步 x_pred = f_func(x_est(:, k-1)); F = compute_jacobian(f_func, x_est(:, k-1)); % 数值雅可比 P_pred = F * P_est * F' + Q; % 滤波步 z_pred = h_func(x_pred); H = compute_jacobian(h_func, x_pred); K = P_pred * H' / (H * P_pred * H' + R); x_est(:, k) = x_pred + K * (z_meas(:, k) - z_pred); P_est = (eye(n) - K * H) * P_pred; end

3.2 雅可比矩阵的数值求法：省心又不容易出错

很多教材里把雅可比矩阵写成解析表达式，看起来很漂亮，实际写代码时非常容易出岔子——因为电力系统量测方程和状态方程的表达式很长，一个符号写错，整个滤波器就发散，排查起来非常痛苦。

我在工程里一直用数值差分求雅可比。Matlab里可以用forward difference或者central difference，代码非常短：

function J = compute_jacobian(f, x) n = length(x); J = zeros(length(f(x)), n); eps_val = 1e-6; for i = 1:n x_plus = x; x_minus = x; x_plus(i) = x(i) + eps_val; x_minus(i) = x(i) - eps_val; J(:, i) = (f(x_plus) - f(x_minus)) / (2 * eps_val); end end

注意这里用的是中心差分，精度比向前差分高一个数量级，实测下来对滤波稳定性很有帮助。同时要注意，(eps_val)的取值不能太大也不能太小，(1e-6)是个经验上比较合适的值，太小会导致数值舍入误差变大。

3.3 UKF的Sigma点生成与权重计算

UKF避免求雅可比矩阵，代价是要生成Sigma点并逐个通过非线性函数传播。生成规则如下（无迹变换）：

假设状态维度为(n)，在(k-1)时刻的状态均值(\hat{x}{k-1})和协方差(P{k-1})已知，生成(2n+1)个Sigma点：

χ^(0) = x̄ χ^(i) = x̄ + sqrt((n + λ) * P) 的第i列, i = 1, ..., n χ^(n+i) = x̄ - sqrt((n + λ) * P) 的第i列, i = 1, ..., n

这里(\lambda = \alpha^2(n + \kappa) - n)，(\alpha)控制Sigma点的分布范围，一般取(1e-3)到(1)之间；(\kappa)是次级缩放参数，通常取(0)或(3-n)。权重则分为均值权重和协方差权重两组，具体公式在标准参考资料里都有。

Matlab里生成Sigma点的代码大概是：

function sigma_points = generate_sigma_points(x, P, lambda) n = length(x); num_sigma = 2 * n + 1; sigma_points = zeros(n, num_sigma); sqrt_matrix = sqrtm((n + lambda) * P); sigma_points(:, 1) = x; for i = 1:n sigma_points(:, i+1) = x + sqrt_matrix(:, i); sigma_points(:, i+n+1) = x - sqrt_matrix(:, i); end end

注意这里用sqrtm而不是逐元素开的sqrt，因为协方差矩阵不是对角阵，必须用矩阵平方根。这是新手最容易踩的坑——用sqrt(P)去生成Sigma点，结果全错。

3.4 UKF的预测与更新流程

生成Sigma点后，后续步骤就是把每个点分别通过状态方程和量测方程传播：

1. 状态传播：(\chi_{k|k-1}^{(i)} = f(\chi_{k-1}^{(i)}))
2. 计算状态预测均值：(\hat{x}{k|k-1} = \sum W_m^{(i)} \chi{k|k-1}^{(i)})
3. 计算预测协方差：(P_{k|k-1} = \sum W_c^{(i)} (\chi_{k|k-1}^{(i)} - \hat{x}_{k|k-1})(...)^T + Q)
4. 对每个状态预测点计算量测预测：(\zeta^{(i)} = h(\chi_{k|k-1}^{(i)}))
5. 量测均值：(\hat{z}_k = \sum W_m^{(i)} \zeta^{(i)})
6. 量测协方差：(P_{zz} = \sum W_c^{(i)} (\zeta^{(i)} - \hat{z}_k)(...)^T + R)
7. 状态与量测互协方差：(P_{xz} = \sum W_c^{(i)} (\chi_{k|k-1}^{(i)} - \hat{x}_{k|k-1})(\zeta^{(i)} - \hat{z}_k)^T)
8. 卡尔曼增益：(K_k = P_{xz} P_{zz}^{-1})
9. 状态更新：(\hat{x}k = \hat{x}{k|k-1} + K_k (z_k - \hat{z}_k))
10. 协方差更新：(P_k = P_{k|k-1} - K_k P_{zz} K_k^T)

UKF的代码量比EKF多不少，但好处是涉及的非线性函数都是直接传入，不用肉眼去算一阶导数。对电力系统这种需要反复修改模型的情况，UKF的可维护性反而更好——改模型的时候，EKF要重新推导雅可比矩阵的解析式，UKF只需要改一下状态方程和量测方程的实现函数即可。

4. 两种算法在同一算例下的实测对比

4.1 仿真算例设置

我用自己的代码在Matlab里做了一个仿真验证，算例采用单机无穷大系统，发电机用经典三阶模型，配一组PMU量测。之所以先用单机无穷大，是因为它的真值方便获得，模型简单、逻辑清晰，适合作为算法验证的起步场景。

仿真时长设为5秒，采样周期(T_s = 0.01s)，一共500个采样点。初始状态加了一定的偏差，模拟滤波器从非精确初始值启动的过程。过程噪声和量测噪声设置如下：

4.2 收敛速度对比

从仿真的第一秒来看，EKF和UKF都能在约20到30个采样点（即0.2到0.3秒）内把初始偏差修正到真实值附近。但两者的收敛轨迹明显不同：UKF的误差曲线下降更平滑，几乎没有超调；EKF在初始阶段会有明显的一两个振荡峰。这个现象的根本原因还是EKF的线性化误差——初始偏差大的时候，工作点离真实值远，泰勒展开的截断误差也随之变大，滤波器做的事就不是"最优修正"而是"边错边纠"。

如果初始偏差更大（比如+50%真实值），EKF会有一定概率直接发散，而UKF依然能稳定收敛。在工程上，这意味着UKF对状态初值不敏感的鲁棒性，确实更适应电网实时估计这种启动条件不可控的场景。

4.3 稳态精度对比

进入稳态后，两种算法的估计结果都能围绕真值波动，但波动幅度有差异。我统计了最后2秒的均方根误差（RMSE）以及最大绝对偏差，结果如下：

从数据上看，UKF的RMSE大约比EKF低了30%到40%。原因其实也很清楚：功角方程和量测方程里的正弦项在正常工作点附近仍然有不可忽略的高阶项，EKF把这些高阶信息全部截断了；UKF通过Sigma点传播，把非线性函数的真实分布特征保留到了二阶矩，精度自然更高。

4.4 计算耗时对比

精度打不过，EKF在计算速度上扳回一城。在相同仿真条件下，我用tic/toc记录了单步滤波的平均耗时：

• EKF单步耗时：约0.8毫秒
• UKF单步耗时：约2.5毫秒

UKF单步计算量大约是EKF的三倍，主要开销在Sigma点逐个通过非线性函数的传播上。不过这里要强调一下，这个差距是在单机无穷大系统（3维状态、4维量测）这种小规模算例下测得的；如果系统规模增大到几十台发电机，状态维度上百，UKF的Sigma点数量会线性增加，计算开销会进一步放大。对于在线应用，需要结合具体的实时性要求来做权衡。

从实用角度我的建议是：如果是做科研对比、离线分析或者状态维度不高的场景，直接上UKF，精度和鲁棒性都更好；如果是做实时在线估计，并且系统规模很大，EKF仍然有不可替代的计算优势。

5. 调参三个月总结出的避坑经验

5.1 初始协方差矩阵：宁大勿小

滤波器的初始协方差矩阵(P_0)直接决定了它对初始状态误差的"信任程度"。我曾经在测试时把(P_0)设得特别小，结果滤波器认为初始状态非常准确，对量测的修正几乎不响应，导致整个估计曲线拖了很久才收敛。

经验法则：(P_0)的对角元素可以设置成初始状态不确定度的平方，但宁可偏大不要偏小。如果完全不知道初值精度，取真实值相关数量级的平方就差不多。比如功角的初始不确定度设为0.1 rad，那(P_0)对应的对角元素就取(0.1^2 = 0.01)。

5.2 Q矩阵和R矩阵的比值更重要

新手最容易犯的错是把Q和R分开一个个调，调了半天毫无头绪。实际上，卡尔曼滤波的稳态表现主要取决于两者的比值，不是绝对值。你把Q和R同时放大10倍，稳态增益基本不变；但如果你只调其中一个，滤波器行为就会产生剧烈变化。

我的调参流程是：

1. 先根据表计精度确定R矩阵，这部分是物理量，PMU的说明书中通常有精度等级
2. 再反复调整Q矩阵的数值，观察估计曲线的平滑度和响应速度
3. Q太小，估计曲线会追着量测噪声跑，毛刺多；Q太大，估计曲线过于平滑，真实动态会被抹掉

有一种"量化检查"的方法：在稳态段，如果量测残差（innovation）的均值明显不为零，说明Q可能偏大或者模型有偏；如果量测残差的标准差显著小于R矩阵对应的量测噪声标准差，说明Q偏小，滤波器在过度信任模型。

5.3 协方差矩阵失去正定性：发散的前兆

EKF和UKF在数值实现时，一个非常隐蔽的问题是在连续迭代中协方差矩阵逐渐失去对称性、甚至失去正定性，导致后续步骤中出现负方差、滤波发散。这类问题在长时间仿真时特别容易出现，因为浮点舍入误差会不断累积。

我养成了一个习惯：每步滤波结束之后，对称化处理协方差矩阵：

P_est = (P_est + P_est') / 2;

如果发现某一步(P)出现了负特征值，还可以加一个对角修饰（比如用Matlab里的nearestSPD函数找到最近的正定矩阵），保证滤波器鲁棒性。这个处理看似朴素，但它救了我好几次长时间仿真的结果。

5.4 状态方程和量测方程的一致性检查

EKF最怕的是状态方程里面写错了物理量纲，导致量测方程里反推的状态值和真值对不上——这种Bug表面上看是"滤波器不收敛"，实际上是模型本身自洽性出了问题。

我做代码验证时的必做动作是：先用没有任何噪声和初始误差的"完美数据"去跑一遍滤波器。如果在这种情况下估计误差还很大，那一定是状态方程或量测方程写错了——因为理论上来讲，在零噪声、零初值误差的前提下，卡尔曼滤波器应该给出零误差的估计结果。通过这个测试，能把模型层面的错误和算法层面的问题分离开来。

5.5 采样时间的敏感性

采样时间(T_s)的选取也会影响滤波器性能。(T_s)太大，离散化误差会直接变成"模型误差"，相当于变相增大了过程噪声Q；(T_s)太小，计算量显著上升，而且量测之间的时间窗口太短，可能让量测对状态修正的贡献变小。

我建议做动态估计时先把采样时间固定为与PMU数据帧率一致的0.01到0.02秒，后续用同一套噪声参数在不同(T_s)下做对比敏感性测试，这样能明确看出离散化误差的比重。我在仿真中发现(T_s)从0.01秒放大到0.05秒时，EKF的RMSE会增大约3倍，而UKF只增大约1.5倍——UKF的Sigma点传播方式对离散化误差的敏感度天然更低。

5.6 从单机到多机的扩展思路

最后说一句扩展经验。单机无穷大系统验证通过之后，往多机系统扩展时，最需要注意的就是状态方程里发电机的互联项——每台发电机的电磁功率表达式里会包含相邻发电机的功角差信息，状态方程比单机模型复杂很多。

从代码层面，推荐把每台发电机的状态方程封装成独立的子函数，再用循环统一调用，这样既方便跟单机模型对照验证，也为后续引入励磁系统、调速器等动态模型留好接口。EKF在多机下的雅可比矩阵如果用手推解析式，工作量会暴涨；数值雅可比配合封装好的状态方程，代码几乎不需要改动。这也是我在EKF和UKF之间，后期越来越倾向UKF的原因——对模型迭代更友好。

我在实际项目中，最开始就是老老实实用EKF跑通了单机模型，后来扩展到一个5机测试系统后，EKF的调参和维护开始变得吃力。换成UKF之后，同样规模和场景下开发效率明显上去了一截。倒不是说EKF不行，而是"不用推导雅可比矩阵"这个优势，在你反复改模型、试参数的时候，真的能省下太多时间。

最后再分享一个小技巧：别把滤波器的流程代码跟模型代码混在一起写。哪怕最初只是做一个简单的单机估计，也把状态方程函数、量测方程函数、EKF核心循环、UKF核心循环、参数配置、数据绘图分开成独立的文件。一旦后续模型复杂度上来了，你会发现这些早期花在代码结构上的时间，会在调试和扩展阶段成倍地还回来。我现在所有动态估计相关项目都还沿用着那套最初的函数拆分习惯，回头来看，这个决定比选EKF还是UKF本身更值钱。