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

电力系统动态状态估计里，EKF和UKF属于那种"名字人人都听过，但真正能把Matlab代码跑通并解释清楚每一步在干什么"的人其实不多的工作。我刚把整套仿真完整实现了一遍，从发电机动态模型、量测方程到两种滤波器的递推代码、再到故障扰动场景下的对比实验，中间踩了不少坑。这篇文章就把这条线从头到尾串起来，给同样在搞这个方向的同学一个可以直接参考的完整版本。

这篇内容适合两类人：一类是正在做课程设计或毕业设计，需要从零搭建基于EKF/UKF的电力系统动态状态估计仿真框架；另一类是刚接触调度自动化，想搞明白动态状态估计和传统静态估计到底差在哪里的工程技术人员。我会把数学模型、核心代码、参数设置思路和实际操作中容易翻车的细节都写出来，代码尽量只用Matlab基础函数，R2020之后的版本基本直接能跑。

1. 为什么动态状态估计比潮流计算的"升级版"更值得关注

1.1 静态估计的局限：一个时间断面解决不了全过程

很多刚接触这个领域的人会有一个疑问：电网调度中心不是早就有状态估计了吗？为什么还要研究动态状态估计？确实，传统EMS里的状态估计已经运行了很多年，但它的核心是基于加权最小二乘（WLS）在一个时间断面上求解。简单说，就是拿某一时刻所有量测数据，去拟合出一个满足潮流方程的系统状态解。

这种做法在稳态运行下没有问题，但它的短板也很明显：每个断面独立计算，完全不利用时间维度的信息。一旦系统进入动态过程，比如负荷突增、发电机跳闸、新能源出力快速波动，WLS只能给你一个又一个孤立的"快照"，而且每次都需要重新迭代求解。更麻烦的是，当量测存在不良数据或者某几条线路量测缺失时，静态估计的结果会明显畸变，又没有历史轨迹可以参照纠正。

1.2 动态估计的定位：在时间轴上"预报"加"校正"

动态状态估计解决的是另一个层面的问题：它不再把每个时刻当成独立断面，而是把系统状态看作一条随时间连续演化的轨迹，用上一拍的状态加上系统动态模型来预测下一拍的状态，再用当前量测对预测值进行修正。这就是经典的"预测-校正"闭环结构。

这种思路带来两个直接好处：第一，量测短暂丢失时，估计器还能依靠模型预报顶着走一段，不会马上失效；第二，量测噪声被模型平滑掉很大一部分，估计结果比单断面拟合更干净。而真正让动态状态估计在工程上变得可行，靠的是PMU同步相量量测的普及，高帧率同步数据让时间维度的建模有了用武之地。

1.3 为什么选EKF和UKF而不直接上粒子滤波

电力系统模型本质上是非线性的，所以线性卡尔曼滤波没法直接用。非线性系统滤波最经典的三个路线是EKF、UKF和粒子滤波（PF）。粒子滤波精度高但计算量巨大，在PMU百赫兹量测速率下做实时估计基本不现实。EKF和UKF在中小规模系统上可以满足计算实时性要求，Matlab实现也相对成熟，所以科研和工程验证里最常见的就是它们两个。

拿我自己实现的场景来说，单机无穷大系统状态只有三维，UKF的计算量比EKF大三四倍，但依然在毫秒级，完全不影响实时性。这也是为什么我最终选定EKF和UKF作为对比方案——两者在原理上刚好是两个方向：一个靠解析求导做线性化，一个靠采样点传播非线性，放在一起对比，恰好能说明很多滤波器的本质问题。

2. 发电机动态模型与量测方程：先把被估计对象搞干净

2.1 三阶实用模型：从转子运动方程到暂态电势

动态状态估计的第一步不是写滤波器，而是把被估计对象的数学模型写清楚。我选的是最常用的发电机三阶实用模型，状态向量取三个物理量：功角δ、转速偏差Δω、q轴暂态电势Eq'。写成连续时间状态方程就是：

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

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

[ \frac{dE_q'}{dt} = \frac{1}{T_{d0}'}(E_{fd} - E_q' - (x_d - x_d')I_d) ]

这里面几个参数的物理含义要清楚：M是发电机惯性时间常数折算后的标幺惯性，D是阻尼系数，Td0'是励磁绕组时间常数，x_d和x_d'分别是同步电抗和暂态电抗，E_fd是励磁电压。用生活化类比来说，第一行方程是说功角的变化率就是转速偏差；第二行方程就是牛顿第二定律在转子上的体现——机械功率和电磁功率不平衡就产生加速度；第三行描述的是励磁绕组这个"大电感"里的磁链不能突变，所以暂态电势是逐渐变化的。

这套模型抓取了转子运动动态和励磁动态两个关键时间尺度，在动态状态估计的入门和中等精度场景下都非常适用。

2.2 单机无穷大系统下的量测方程

状态方程决定了状态怎么演化，量测方程决定了我们从外部能观测到什么。在我的仿真中，量测量取PMU最容易提供的两个电气量：发电机输出的有功功率P_e和无功功率Q_e。在单机无穷大系统的简化模型下，它们和状态变量之间有明确的解析关系：

[ P_e = \frac{E_q' V_b}{X_{\Sigma}} \sin\delta ]

[ Q_e = \frac{E_q'^2 - E_q' V_b \cos\delta}{X_{\Sigma}} ]

其中V_b是无穷大母线电压，X_Σ是暂态电抗加上变压器和线路的总电抗。这两个式子是从发电机注入无穷大母线的复功率推导出来的，推导过程不复杂：先写出电流相量，再做共轭相乘取实部虚部就能得到。关键是它告诉了我们一件事——量测方程是强非线性的，δ通过正弦函数进入P_e，这让滤波器的设计无法回避非线性处理。

我在实现时把这两个量测都加了独立高斯噪声，标准差取0.005标幺值。PMU的典型幅值测量精度一般在0.1%到0.5%之间，这个噪声水平大致对应中等偏严苛的工况，不会让滤波变成纯粹的"拟合噪声"，也不会让量测完美得看不出滤波器的差别。

2.3 连续方程的离散化：离散步长决定滤波器性能上限

扩展卡尔曼滤波和无迹卡尔曼滤波处理的对象是离散时间系统，所以连续状态方程必须离散化。最常见的做法是一阶欧拉近似：

[ x_{k+1} \approx x_k + T_s f(x_k) ]

还有一种更精细的做法是二阶Runge-Kutta（也叫Heun方法）：

[ k_1 = f(x_k), \quad k_2 = f(x_k + T_s k_1) ]

[ x_{k+1} = x_k + \frac{T_s}{2}(k_1 + k_2) ]

采样步长Ts的选择直接影响滤波器性能上限。我用的是Ts=10ms，对应PMU常见的100Hz量测帧率。在这个步长下，一阶欧拉也能跑出不错的效果，但我在滤波器内部统一用了二阶Runge-Kutta来传播状态，这样可以把数值积分误差从二阶降到三阶，避免离散化误差被误认为是滤波器性能差。

2.4 噪声矩阵Q、R和初始协方差P0的物理标定

卡尔曼滤波器的调参核心就是三个矩阵：过程噪声协方差Q、量测噪声协方差R、初始状态协方差P0。很多新手在这里随便给一组数，结果滤波效果一塌糊涂还不知道问题出在哪。我的经验是从物理意义出发给初始值。

量测噪声R最容易定，它直接对应PMU的测量精度，我取P_e和Q_e的标准差都是0.005标幺，所以R=diag([2.5e-5, 2.5e-5])。过程噪声Q反映的是模型一步预测的不确定度：比如δ在两个采样点之间的变化量大致是Δω×Ts，如果Δω正常波动范围在0.05 rad/s附近，那么δ的预测不确定度约0.0005 rad，方差量级就是2.5e-7；但实际工程中还有负荷随机波动等未建模动态，所以Q不能取得太小，我给的是Q=diag([2.5e-5, 1e-4, 2.5e-5])，比纯理论值放宽一到两个数量级。P0则根据初始状态误差的估计来定，如果功角初始偏差可能达到0.2 rad，那么P0(1,1)取0.04量级。

3. EKF实现细节：每一步都在算雅可比矩阵

3.1 EKF的标准递推流程

扩展卡尔曼滤波的本质就是把非线性模型在当前状态附近做一阶泰勒展开，然后用线性卡尔曼滤波的框架做预测和更新。它的标准流程是五步：状态预测、协方差预测、量测预测、滤波增益计算、状态与协方差更新。用公式写出来就是：

[ x_{k|k-1} = f(x_{k-1}) ]

[ P_{k|k-1} = F_k P_{k-1} F_k^T + Q ]

[ K_k = P_{k|k-1} H_k^T (H_k P_{k|k-1} H_k^T + R)^{-1} ]

[ x_k = x_{k|k-1} + K_k (z_k - h(x_{k|k-1})) ]

[ P_k = (I - K_k H_k) P_{k|k-1} ]

这里F_k和H_k分别是状态方程和量测方程在最近估计点处的雅可比矩阵。注意一个顺序问题：F_k在上一时刻的滤波值x_{k-1}处计算，而H_k应该在预测值x_{k|k-1}处计算。这个顺序我一开始搞反过，导致初始阶段误差明显偏大——线性化点选的不一样，结果会差很多。

3.2 雅可比矩阵的解析推导：F和H到底长什么样

EKF实现最繁琐的部分就是求雅可比矩阵。以我用的三阶发电机模型为例，连续状态方程的雅可比矩阵A是3×3矩阵，在点x=[δ, Δω, Eq']处求值为：

[ A = \begin{bmatrix} 0 & 1 & 0 \ -\frac{E_q' V_b \cos\delta}{M X_\Sigma} & -\frac{D}{M} & -\frac{V_b \sin\delta}{M X_\Sigma} \ -\frac{(x_d - x_d') V_b \sin\delta}{T_{d0}' X_\Sigma} & 0 & -\frac{1}{T_{d0}'} - \frac{x_d - x_d'}{T_{d0}' X_\Sigma} \end{bmatrix} ]

离散化后的状态转移矩阵F用一阶近似就是F≈I+Ts×A。量测雅可比H是2×3矩阵，由P_e和Q_e对状态求偏导得到：

[ H = \begin{bmatrix} \frac{E_q' V_b \cos\delta}{X_\Sigma} & 0 & \frac{V_b \sin\delta}{X_\Sigma} \ \frac{E_q' V_b \sin\delta}{X_\Sigma} & 0 & \frac{2E_q' - V_b \cos\delta}{X_\Sigma} \end{bmatrix} ]

这几个偏导数推导不难但容易出错，尤其是Q_e对Eq'的偏导，很容易丢掉2Eq'项。建议推完先用数值差分验证一遍，再放进滤波器。

3.3 EKF的Matlab核心代码骨架

下面是我用的EKF主循环核心代码，注释写在关键位置：

nx = 3; % 状态维数 Ts = 0.01; % 采样周期 10ms N = length(t); % 总仿真步数 % 连续状态方程 f f = @(x) [ x(2); (Pm - x(3)*Vb*sin(x(1))/Xsum - D*x(2)) / M; (Efd - x(3) - (xd - xdp)*(x(3) - Vb*cos(x(1)))/Xsum) / Tdo ]; % 量测方程 h h = @(x) [ x(3)*Vb*sin(x(1))/Xsum; (x(3)^2 - x(3)*Vb*cos(x(1)))/Xsum ]; % 状态转移雅可比 A A = @(x) [ 0, 1, 0; -x(3)*Vb*cos(x(1))/(M*Xsum), -D/M, -Vb*sin(x(1))/(M*Xsum); -(xd-xdp)*Vb*sin(x(1))/(Tdo*Xsum), 0, -1/Tdo-(xd-xdp)/(Tdo*Xsum) ]; % 量测雅可比 H Hfun = @(x) [ x(3)*Vb*cos(x(1))/Xsum, 0, Vb*sin(x(1))/Xsum; x(3)*Vb*sin(x(1))/Xsum, 0, (2*x(3)-Vb*cos(x(1)))/Xsum ]; x_est = x0_est; P = P0; for k = 1:N-1 % 当前时刻量测（仿真中由真值加噪声生成） z = h(x_true(:,k+1)) + sqrt(R)*randn(2,1); % 预测 x_pred = x_est + Ts*f(x_est); % 欧拉预测 F = eye(nx) + Ts*A(x_est); % 离散化状态转移矩阵 P_pred = F*P*F' + Q; % 更新 H = Hfun(x_pred); % 在线性化点 x_pred 处求 H S = H*P_pred*H' + R; K = P_pred*H'/S; x_est = x_pred + K*(z - h(x_pred)); % 状态更新 P = (eye(nx) - K*H)*P_pred; % 协方差更新 end

注意上面代码在生成量测时用了sqrt(R)*randn(2,1)，这样就不需要额外的统计工具箱函数。另外，我在代码里故意用欧拉法做状态预测，这是为了展示最基本结构；实际在项目里建议把x_pred那行换成二阶Runge-Kutta传播函数。

3.4 EKF的三个经典翻车点

第一是线性化误差。EKF只保留一阶项，在功角摆动幅度大的时候，正弦函数在非线性最强处（δ在π/2附近）被线性近似后误差会被放大，如果此时过程噪声Q又给得小，滤波器会过于相信模型预测，导致估计偏差迟迟消不掉。

第二是协方差矩阵破坏对称性。标准更新公式P=(I-KH)P_pred在数值上每次运算都会引入微小不对称，长期运行后P可能变成非正定矩阵，Cholesky分解直接报错。解决方法是改用Joseph形式：

[ P_k = (I - K_k H_k)P_{k|k-1}(I - K_k H_k)^T + K_k R K_k^T ]

这个形式在数值上更稳定，代价只是多几次矩阵乘加。

第三是初值给得太离谱时EKF很容易发散。如果功角初始偏差超过30度，EKF的迭代初始段会出现明显震荡甚至直接发散，这件事在后面做对比实验时体现得特别明显。

4. UKF实现细节：用sigma点把非线性"原样搬过去"

4.1 无迹变换的核心思想

UKF和EKF的根本差异在于如何处理非线性。EKF的做法是"我先把非线性函数线性化，再用线性高斯框架算"；UKF的做法是"我不去近似非线性函数，而是去近似状态的分布"。具体手段就是无迹变换——选取一组精心设计的采样点（sigma点），让这些点的均值和协方差严格等于当前状态的均值和协方差，然后把这组点原封不动地通过非线性函数传播，再从传播后的点中重新统计出新的均值和协方差。

这就像评估一支队伍翻过一座山后的集合位置：EKF是拿队伍中心点的高度和坡度估算整支队伍翻山后的位置，UKF则是让几个有代表性的队员实际翻一遍山，再用他们的分布推断整支队伍的情况。显然，后者对强非线性更稳健。

4.2 sigma点生成与权重参数选择

对于n维状态，UKF需要2n+1个sigma点。生成方式如下：先对协方差矩阵做Cholesky分解，得到下三角矩阵S，令：

[ \chi_0 = x, \quad \chi_i = x + \sqrt{n+\lambda}, S_i, \quad \chi_{i+n} = x - \sqrt{n+\lambda}, S_i ]

这里的下标i表示矩阵的第i列。缩放参数λ由三个参数共同决定：α、β、κ。我用的经典组合是α=1e-3，β=2，κ=0，于是λ=α²(n+κ)−n≈−3。α控制sigma点离中心点的距离，取小值让采样点贴近中心，适合大多数平滑非线性系统；β是状态分布的先验信息项，对高斯分布取2是经验最优；κ一般取0或3−n，保证协方差半正定。

权重按照以下公式计算：

[ W_0^m = \frac{\lambda}{n+\lambda}, \quad 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 ]

4.3 UKF的Matlab核心代码骨架

UKF的代码结构比EKF直观不少，因为不需要手动求雅可比。核心预测和更新代码如下：

n = 3; m = 2; alpha = 1e-3; beta = 2; kappa = 0; lambda = alpha^2*(n+kappa) - n; % 权重 Wm = zeros(1,2*n+1); Wc = zeros(1,2*n+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 x_est = x0_est; P = P0; for k = 1:N-1 z = h(x_true(:,k+1)) + sqrt(R)*randn(2,1); % 生成sigma点 S = chol((n+lambda)*P, 'lower'); chi = [x_est, x_est + S, x_est - S]; % 通过状态方程传播sigma点 Xpred = zeros(n, 2*n+1); for i = 1:2*n+1 Xpred(:,i) = rk2_step(f, chi(:,i), Ts); % 二阶Runge-Kutta end % 加权统计预测均值和协方差 x_pred = Xpred * Wm'; P_pred = zeros(n,n); for i = 1:2*n+1 dx = Xpred(:,i) - x_pred; P_pred = P_pred + Wc(i)*(dx*dx'); end P_pred = P_pred + Q; % 通过量测方程传播sigma点 Zpred = zeros(m, 2*n+1); for i = 1:2*n+1 Zpred(:,i) = h(Xpred(:,i)); end z_pred = Zpred * Wm'; % 计算互协方差和滤波增益 Pzz = zeros(m,m); Pxz = zeros(n,m); for i = 1:2*n+1 dz = Zpred(:,i) - z_pred; dx = Xpred(:,i) - x_pred; Pzz = Pzz + Wc(i)*(dz*dz'); Pxz = Pxz + Wc(i)*(dx*dz'); end Pzz = Pzz + R; K = Pxz / Pzz; % 更新状态与协方差 x_est = x_pred + K*(z - z_pred); P = P_pred - K*Pzz*K'; end

上面代码里rk2_step是实现二阶Runge-Kutta传播的子函数，本质上就是把前面提到的k1、k2算一遍。量测更新部分尤其要注意Pzz必须叠加R，否则滤波增益会偏大，导致估计结果过度依赖量测，噪声滤不干净。

4.4 UKF和EKF的工程取舍对比

两种滤波器放在一起对比，没有绝对的优劣，只有适不适合。下面这个表是我在单机无穷大系统上实际对比后的直观感受：

实际项目里，如果模型比较温和、初始状态有把握、实时性要求极高，EKF完全够用。但如果系统会经历大扰动、模型非线性强、或者你不想每次改模型都痛苦地推雅可比矩阵，UKF是更稳的选择。

5. 仿真实验设计：单机无穷大系统上的正面交锋

5.1 仿真场景与扰动设置

我用单机无穷大系统做对比实验，发电机参数按经典值设置：H=5.0s，即M=2H/ωs≈0.0318（50Hz系统，ωs=314.16），x_d=1.8，x_d'=0.3，Td0'=6.0s，X_Σ=0.6，V_b=1.0，D=2。系统初始稳态在P_m=0.8运行，0.05s处P_m从0.8阶跃到1.2，模拟负荷突增约50%。这个扰动会让功角从初始约0.3 rad摆开到1 rad以上，发电机经历一个明显的动态振荡过程。

仿真时长取4s，采样步长Ts=10ms，总步数400步。真值轨迹用ode45生成，但在离散时刻采样；滤波器只拿到带噪声的量测，不直接看到真值。这里有个细节：真值用高精度求解器生成，而滤波器内部用二阶Runge-Kutta传播，两者之间存在离散化失配，这其实是刻意为之——回顾实际工程，模型总是不完美的，与其让滤波器"完美匹配"数据生成模型，不如让它在轻微失配下接受检验。

滤波器初始值故意给偏：真实初始状态是[0.3, 0, 1.0]，滤波器初始估计设为[0.45, 0.02, 1.15]，P0取diag([0.04, 0.0025, 0.04])，相当于一开始对功角的把握误差在0.15 rad左右。两组滤波器使用完全相同的随机量测序列，保证对比公平。

5.2 正常扰动下的估计结果对比

在P_m阶跃场景下，EKF和UKF都能跟上真实轨迹，但细节上差异明显。初始收敛阶段，EKF的功角估计误差峰值大约0.12 rad，UKF大约0.06 rad，UKF收敛到稳态误差带的时间比EKF早了将近0.1s。在功角摆到最大位置附近（δ超过1.0 rad），EKF出现了明显的滞后跟踪——这是线性化误差在正弦非线性最强处的直接体现；而UKF在这个区间的跟踪误差只有EKF的一半左右。

我统计了两种滤波器在完整仿真时长内的RMSE（均方根误差），趋势是：功角RMSE方面，EKF约为0.045 rad，UKF约为0.021 rad；转速偏差RMSE方面，EKF约为0.018 rad/s，UKF约为0.009 rad/s；暂态电势RMSE方面，EKF约为0.035 p.u.，UKF约为0.014 p.u.。UKF在各个状态量上的估计误差总体小一半左右。

5.3 把初始误差拉大：EKF和UKF的鲁棒性差距

为了测鲁棒性，我把滤波器初始功角估计改为0.6 rad，偏差达到0.3 rad，同时把P0相应放大。这个条件下EKF的前50步出现了明显的波动，功角估计一度被拉向错误方向，然后才慢慢拉回来；UKF则从第一步起就稳定朝真值收敛，没有出现大幅摆动。

我的实际体会是：在模型失配或量测短暂异常的情况下，UKF的sigma点传播机制天然具备更强的"纠错"能力。因为EKF把分布压缩成了一条线，在线性化点上做局部近似，一旦线性化点本身偏离真值较远，后续估计很容易朝错误方向惯性滑行；而UKF用一组分布点去捕捉非线性映射，即使中心点偏了，点的整体分布也能提供更充分的梯度信息。

5.4 滤波器调参的实战心得

调参是这个项目里最耗时间的部分。一个最典型的坑是：Q给得越小，稳态轨迹越平滑，但一旦发生真实扰动，滤波器反应越迟钝，甚至出现发散。我试过Q=diag([1e-6, 1e-6, 1e-6])的极端情况，P_m阶跃后估计值需要将近0.5s才跟上真值，明显滞后；反过来把Q放大到Q=diag([1e-1, 1e-1, 1e-1])，跟踪是快了，但滤波后的轨迹满是量测噪声的毛刺，状态曲线完全不平滑。

调参建议是从物理量级出发，先给合理的Q和R，然后固定R只调Q，观察均方根误差和轨迹平滑度的平衡。如果跟踪慢、误差先大后小，说明Q偏小；如果曲线毛糙、稳态误差波动大，说明Q偏大。另外要始终保证R的量级和量测噪声标准差的平方匹配，R给大了就是"不信量测"，给小了就是"过度相信量测"，两种极端在仿真结果里都一眼能看出来。

6. 实战中容易踩的坑与扩展方向

6.1 协方差矩阵数值恶化：最容易被忽视的隐性故障

卡尔曼滤波在长期运行时最容易出现的隐性故障就是协方差矩阵失去正定性。症状是程序跑着跑着突然报错，说Cholesky分解失败，但之前一切正常。EKF里更常见，因为标准更新公式P=(I−KH)P_pred本身不能保证对称，浮点运算的微小不对称在几百个步长后会被放大。

我的处理习惯是两件事同时做：第一，每次更新后强制对称化，P=(P+P')/2；第二，把协方差更新改成Joseph形式。虽然增加了一点计算量，但长期运行稳定很多。还有一个更简单的兜底办法是给P_pred加一个极小的对角增量εI，数值上托住正定性，ε取1e-12量级就够。

6.2 滤波发散的检测与应对：别等到误差爆炸才发现

滤波发散的典型表现是：估计误差不收敛，持续增大，而协方差P却不断缩小，滤波器越来越"自信地错误"。这是因为量测更新长期不被采纳，或者模型预测和真实动态之间出现结构性偏差。最简单的检测方法是监控新息序列（残差）的统计特性。理论上归一化新息应该服从标准正态分布，如果残差均值明显偏离0、方差持续超过预期，基本可以判定发散或即将发散。

应对策略有三种：一是适当增大Q，承认模型不确定；二是引入自适应Q调度，比如检测到残差增大时按比例放大Q，这也是自适应卡尔曼滤波的常见做法；三是检查模型和量测里有没有未建模的突变，比如我仿真里P_m阶跃导致的模型失配，靠滤波本身是消化不了的，必须从源头修正。

6.3 从单机到多机系统的扩展思路

单机无穷大系统的代码跑通后，扩展到多机系统是自然的方向。多机系统的状态向量是所有发电机的功角、转速、暂态电势拼接起来的高维向量，维度从3变成3N_g。这里有两个现实问题：第一，EKF的状态转移雅可比矩阵从3×3变成3N_g×3N_g，推导量和计算量都显著增加；第二，UKF的sigma点数量是2n+1，状态维数一旦超过20，每步要传播的sigma点数量就超过41个，实时性开始吃紧。

工程上比较常用的折中是降维处理：把互联网络等效处理，忽略部分厂站间的动态耦合，或者采用分区分布式估计，每个区域只估计本地几台机，区域间交换边界信息。这种做法在调度自动化中有明显的工程价值，但实现复杂度比单机Demo高一个量级，建议先把单机的原理吃透再动手。

6.4 进阶方向：自适应、鲁棒化和多算法融合

如果只停留在EKF/UKF本身，理解深度会受限。我建议下一步往三个方向走：第一是自适应卡尔曼滤波，通过协方差匹配在线调整Q和R，解决固定噪声参数无法适应工况变化的问题；第二是鲁棒滤波，比如H∞滤波、基于M估计的鲁棒卡尔曼滤波，专门对付量测中的非高斯异常值；第三是粒子滤波与UKF的混合结构，用粒子滤波处理重尾噪声场景下的状态突变，用UKF加速局部采样。

这三个方向的Matlab代码都可以在当前框架上逐步叠加，不会把前面的工作推翻重来。我的经验是，每扩展一个方向，都要先在新息序列的可视化上下功夫——把残差、P对角元素、估计误差画在同一张图里，你能看到很多理论推导里看不到的动态行为，这也是调整滤波器参数最直接的依据。

最后聊一点实际操作层面的体会：这套EKF/UKF仿真跑通并不难，难点在于理解每一个矩阵更新背后的物理含义。EKF的雅可比矩阵里每一项都有对应的电气量变化率，UKF的sigma点里每一个都是发电机的一个可能的运行状态。把这些对应关系理清楚之后，你会发现不管是改成多机系统，还是换一种滤波器，代码结构和调试思路基本都是相通的。如果读完这篇文章你想自己动手，建议先用代码把单机无穷大系统的完整轨迹画出来，再慢慢往里面加扰动、加对比指标、加自适应模块，这条路走通之后，动态状态估计里最核心的那些直觉你就都有了。