1. 项目到底在解什么题：动态状态估计的工程意义

电力系统的状态估计，说起来是个老话题了。传统的静态状态估计，利用SCADA系统采集的遥测数据，通过加权最小二乘等方法估计系统的节点电压幅值和相角，两三秒甚至更长时间才刷新一次，对于稳态运行分析是够用的。但现在电网的运行环境变化很快，新能源大规模接入、电力电子化设备增多、扰动事件频发，调度人员需要实时看清楚发电机的功角、转速这类动态状态量，才能判断机组是否稳定、是否有失步风险。这时候，静态估计的刷新速度和分析深度就不太够用了，于是动态状态估计（Dynamic State Estimation）被推到了前台——利用高采样率的同步相量测量装置（PMU）数据，在几毫秒到几十毫秒的时间尺度上，持续跟踪发电机等动态元件的状态变化。

动态状态估计解决的核心问题可以概括为一句话：我只有带噪声的电流、电压、功率量测，但我想知道系统内部真正发生了什么——特别是发电机转子这类无法直接测量的内部状态。传统静态估计做不到这一点，而卡尔曼滤波这类递推估计方法天然适合这个任务。扩展卡尔曼滤波（EKF）和无迹卡尔曼滤波（UKF）是其中最经典、最常被拿来做对比的两条路线：EKF靠泰勒展开做局部线性化，UKF靠采样点做统计近似。这个Matlab项目，就是把这两种算法完整实现在一个典型的电力系统动态模型上，跑出估计曲线和误差指标。

这篇博文适合三类人看：一是正在做电力系统动态状态估计课题的研究生，需要一套能跑通、能出图的参考代码；二是学卡尔曼滤波但一直没见过实际工程案例的算法学习者，想知道EKF和UKF到底差在哪；三是做PMU数据应用、宽频测量装置研发的工程师，想评估这两种算法在实际数据上的表现。读完这篇，你应该能搭出一套完整的仿真框架，并且知道每一种参数、每一段代码为什么要这么写。

2. EKF和UKF的原理脉络：两种非线性逼近路线

2.1 EKF的线性化思路与局限

扩展卡尔曼滤波的思路非常直白：既然系统是非线性的，那我就找一个点把它展开成线性的来近似处理。具体做法是对非线性状态转移函数 (f(x)) 和量测函数 (h(x)) 求雅可比矩阵，在每一步递推时用当前状态估计值附近的线性化模型去替代原来的非线性模型，然后套用标准卡尔曼滤波的预测—更新框架。

这个思路在工程上应用极广，IMU姿态解算、GPS/INS组合导航、电池SOC估计，处处都有EKF的身影。但它的短板也很明显：线性化误差的大小完全取决于非线性强度和当前工作点。如果系统在某个环节呈现强非线性——比如电力系统中发电机电磁功率与功角之间的正弦关系，当功角变化幅度较大时，一阶泰勒展开忽略掉的二阶以上项会带来不可忽略的误差。而且雅可比矩阵的计算本身是个麻烦事，推导要花时间，写代码容易出错，某些复杂模型甚至很难解析求导。

2.2 UKF的统计近似思路

无迹卡尔曼滤波换了一条路：不取线性化近似，而是用一定数量的采样点（sigma点）去直接逼近非线性变换后随机变量的统计特性。它背后的理论是无迹变换（Unscented Transform）：假设状态量服从高斯分布，按确定的规则在分布中选取若干个sigma点，把这些点分别通过非线性函数，再对变换后的点加权求均值和协方差。理论上，这种统计近似对非线性函数的逼近精度可以达到三阶泰勒展开的水平，明显优于EKF的一阶近似，而且完全不需要计算雅可比矩阵。

我第一次看UKF的权重公式时也觉得有点绕，但本质上它做的事情很简单：让采样点”跟着分布走“，非线性函数只是把这些点映射到新空间，然后在新空间里重新统计。这就好比你去测一个不规则区域的实际面积，EKF的做法是把区域大致看成一个矩形来估算，UKF则是在区域内选一些有代表性的点，实测坐标后把面积算出来，自然更贴近真实情况。

2.3 EKF vs UKF的特点对比

放在电力系统动态状态估计的背景下对比，两者的差异有几个关键层面：

从这张表能看出来，这个题目把两种算法都做出来，本身就是完整的对比实验。项目中不仅会验证“UKF精度更高”的常见结论，还会量化地给出一组RMSE和计算耗时的数据，这是论文和课题报告中非常有说服力的素材。

3. 仿真模型与Matlab代码骨架

3.1 系统模型与基础参数设置

动态状态估计的仿真通常从单机无穷大系统（SMIB）开始，这套模型简单但物理意义清晰——一台发电机经电抗连接到无穷大母线，既能反映机电暂态的核心动态过程，又方便对比算法的性能差异。

发电机的经典二阶模型（摇摆方程）如下，状态变量选取功角 (\delta) 和转子角速度 (\omega)：

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

[ \frac{d\omega}{dt} = \frac{1}{2H}(T_m - \frac{E'V}{X}\sin\delta - D(\omega - \omega_0)) ]

式中的 (E') 是发电机暂态电抗后的电势，(V) 是无穷大母线电压，(X) 是发电机暂态电抗到无穷大母线的总电抗，(H) 是惯性时间常数，(D) 是阻尼系数，(T_m) 是机械功率。量测量选择最容易获得的电磁功率 (P_e = \frac{E'V}{X}\sin\delta)，实际仿真中还可以叠加一个机端电压幅值量测，但功率量测已经足以体现非线性特征。

Matlab代码里把系统参数集中在一个结构体里管理，方便批量调参：

% 系统参数配置 params.omega0 = 2*pi*50; % 同步角速度 rad/s params.H = 6.0; % 惯性时间常数 s params.D = 2.0; % 阻尼系数 p.u. params.Tm = 1.0; % 机械功率 p.u. params.E = 1.2; % 暂态电势 p.u. params.V = 1.0; % 无穷大母线电压 p.u. params.X = 0.5; % 总电抗 p.u. params.dt = 0.01; % 采样周期 s，工频下对应压缩力学意义明确

采样周期的选择这里多说一句。发电机转子运动的时间常数通常在秒级，而动态状态估计要求至少能捕捉到机电振荡（0.1Hz~2Hz）的动态过程，所以采样周期取0.01s（100Hz）是合理的起点。如果你后续要做与故障暂态更紧密的分析，可能要加密到0.001s，但对应的计算负担和量测噪声特性会完全不一样。

3.2 状态方程与量测方程封装

仿真代码通常把连续时间的状态方程封装成一个函数，这样EKF、UKF以及真实轨迹生成可以共用同一套动力学模型：

function dx = generator_dynamics(x, params) % x = [delta; omega] delta = x(1); omega = x(2); Pe = params.E * params.V * sin(delta) / params.X; dx = zeros(2,1); dx(1) = omega - params.omega0; dx(2) = (params.Tm - Pe - params.D*(omega - params.omega0)) / (2*params.H); end

量测方程则单独封装，因为后续EKF和UKF对量测的处理方式有区别：

function z = generator_measurement(x, params) % 量测：电磁功率 Pe，叠加噪声在仿真主循环中处理 delta = x(1); z = params.E * params.V * sin(delta) / params.X; end

真实轨迹的生成用四阶龙格库塔（ode45或手写RK4），不要用简单的欧拉法，尤其是对比EKF和UKF精度的时候，如果真实轨迹本身误差大，算法的对比结果就不干净了。我习惯在仿真循环里手动写RK4步进，因为ode45的变步长机制配合固定采样时刻的量测更新会带来不便。

3.3 EKF核心代码实现

EKF的Matlab实现分预测和更新两步。预测步需要计算状态转移雅可比矩阵 (F)，在每一时刻基于当前状态估计值重新求偏导：

% EKF 预测步 x_pred = x_est + params.dt * generator_dynamics(x_est, params); % 状态转移雅可比矩阵 F % d(delta)/d(delta) = 0, d(delta)/d(omega) = 1 % d(omega)/d(delta) = -E*V*cos(delta)/(2*H*X) % d(omega)/d(omega) = -D/(2*H) F = [0, 1; -params.E*params.V*cos(x_est(1))/(2*params.H*params.X), ... -params.D/(2*params.H)]; P_pred = F * P_est * F' + Q;

更新步的量测雅可比 (H) 对 (P_e = E'V\sin(\delta)/X) 求导，结果很简洁：

% EKF 更新步 H = params.E * params.V * cos(x_pred(1)) / params.X; K = P_pred * H' / (H * P_pred * H' + R); z_pred = generator_measurement(x_pred, params); x_est = x_pred + K * (z_meas - z_pred); P_est = (eye(2) - K * H) * P_pred;

这里有一个容易踩的坑：量测雅可比 (H) 是一行两列的矩阵，但如果你写成行向量或者列向量搞混了，Matlab的维数检查不会报错，但结果会悄悄出错。我在调试时习惯在每一步递推后用size()检查各矩阵维数，并打印前几个时刻的卡尔曼增益，确认它在一段合理区间内变化，而不是忽大忽小甚至为负。

3.4 UKF核心代码实现

UKF的实现比EKF多了一段sigma点生成和权重计算的代码，但好处是一旦写好，换任何非线性模型都能直接复用。

sigma点生成的核心逻辑是先做Cholesky分解得到协方差矩阵的平方根，再沿每个状态维度以正负方向扩展采样点。对二维状态系统，一共生成 (2n+1=5) 个sigma点：

% UKF 参数 alpha = 1e-3; % 决定sigma点散布范围，通常取小值 beta = 2; % 高斯分布的最优取值 kappa = 0; % 次级尺度参数，一般取0或3-n lambda = alpha^2 * (n + kappa) - n; % 生成sigma点：利用Cholesky分解 sqrtP = chol((n + lambda) * P_est, 'lower'); Sigma = zeros(n, 2*n+1); Sigma(:,1) = x_est; for i = 1:n Sigma(:,i+1) = x_est + sqrtP(:,i); Sigma(:,n+i+1) = x_est - sqrtP(:,i); end % 权重计算 Wm = [lambda/(n+lambda), ones(1,2*n)/(2*(n+lambda))]; Wc = [lambda/(n+lambda)+(1-alpha^2+beta), ones(1,2*n)/(2*(n+lambda))];

sigma点构建好后，预测步就是逐个点跑非线性状态方程，再加权汇总：

% 将每个sigma点通过状态方程步进（同样用RK4或欧拉离散） X_pred = zeros(n, 2*n+1); for i = 1:2*n+1 X_pred(:,i) = Sigma(:,i) + params.dt * generator_dynamics(Sigma(:,i), params); end x_pred = Wm * X_pred'; % 加权求状态均值 P_pred = Q; for i = 1:2*n+1 diff = X_pred(:,i) - x_pred; P_pred = P_pred + Wc(i) * (diff * diff'); end

量测更新类似，把预测的sigma点再过一遍量测方程，计算量测均值、量测协方差和互协方差，再求卡尔曼增益。UKF代码常见的坑在chol函数的输入矩阵必须是对称正定的，如果某一步P_est失去正定性，Cholesky分解会直接报错。解决办法有两个方向：数值上给P_est加上一个小的eye(n)*eps，或者在每个递推周期后检查特征值是否都大于0。

3.5 初始化与噪声矩阵设定

滤波器初值和噪声矩阵的设定直接决定收敛速度和稳态精度。我的经验参数如下：

x0 = [0.1; 314.1593]; % 功角初始值0.1 rad，角速度初值接近同步转速 P0 = diag([0.1^2, 1e-3^2]); % 初始协方差：功角不确定度较大 Q = diag([1e-5, 1e-5]); % 过程噪声：模型误差和未知扰动 R = 0.01^2; % 量测噪声：对应PMU功率量测误差1%

对于 (Q) 和 (R) 的比值，需要说明的是它直接决定卡尔曼增益的倾向：(Q) 相对 (R) 越大，滤波器越相信新量测、响应越快，但稳态误差的波动也更大；反之则更平滑但响应变慢。实际仿真中我第一次跑的时候状态估计迟迟不收敛，把 (Q) 调大一个数量级后明显改善。这里提供一个实用口诀：先调 (P_0) 让曲线收得住，再调 (Q/R) 让曲线跟得紧，最后微调 (R) 让估值的噪声幅度合适。

4. 仿真结果解读与调参心得

4.1 典型仿真结果对比

以一典型场景为例：系统稳态运行到0.2s时发生一个持续0.1s的机械功率阶跃扰动（(T_m) 从1.0跳到1.05），随后恢复。这个扰动会让功角产生一个典型的摆动过程，正好用来考验滤波器对动态过程的跟踪能力。

运行两种算法后分别计算功角估计的RMSE，典型结果大致是：EKF的功角RMSE在0.02~0.04 rad量级，UKF可以压到0.005~0.01 rad量级。角速度的估计精度同样UKF占优，约比EKF高一倍以上。这个差距在扰动发生的瞬态窗口内更加明显——EKF在功角变化速率最快的时候会出现明显的跟踪滞后，因为线性化模型在快速变化的工作点附近近似误差被放大，而UKF的统计近似对这个过程适应得更好。当然，不同系统参数下具体数值会有浮动，但UKF精度优于EKF这个定性结论在大多数中等非线性强度场景下都是成立的，这也符合理论预期。

计算量方面，UKF每步递推需要把5个sigma点分别通过状态方程和量测方程，而EKF只需要一个点加上雅可比矩阵运算，所以UKF单步耗时通常是EKF的1.5~2倍。但在2维状态空间下，这个差距完全可以忽略，一步也就零点几毫秒的事。如果扩展到几十维状态的多机系统，这个差距才需要认真对待，届时就可能要考虑降维或者用平方根UKF改善数值稳定性。

4.2 核心调参实战心得

调参是这类仿真项目中最耗时间的环节，我把自己实际踩过的几个规律总结成下面几条：

第一，如果滤波器发散，先怀疑数值稳定性而不是算法本身。我遇到过的EKF发散案例中，一半以上是因为 (P) 矩阵不断迭代后失去对称正定性，导致更新步算出的协方差出现负数对角线元素。解决办法是在每一步递推后强制对称化：P = (P+P')/2，必要时加P = P + 1e-12*eye(n)。这与系统本身无关，但很多人栽在这里很久出不来。

第二，初始协方差 (P_0) 不要取得太大。理论上 (P_0) 反映我们对初始状态的不确定度，取大一点可以加速收敛。但过大的 (P_0) 会让最初的卡尔曼增益异常偏大，导致状态估计值在头几个采样点剧烈跳变，甚至反过来破坏收敛过程。稳妥的做法是先根据工况大致确定功角和角速度的合理范围，取范围上限的1/3平方作为对应的 (P_0) 对角线值。

第三，过程噪声 (Q) 的物理含义要想清楚。很多人把 (Q) 当成一个随意设的调节旋钮，但它实际代表的是模型没有捕捉到的随机扰动功率，比如机械功率波动、电气参数变化。如果你加上系统的量测噪声是白噪声而模型是精确的，(Q) 理论上应该很小，取 (10^{-6}) 左右都不过分。但实际模型和真实系统总存在失配，(Q) 太小时滤波器会过度相信模型预测，导致估计结果跟不上真实的动态变化，表现为持续的固定偏差。遇到这种情况，把 (Q) 逐步调大，直到估计结果对新息有足够响应但噪声又不过分放大为止。

第四，评价指标不要只看RMSE。论文里写RMSE当然没问题，但你自己判断算法好坏时，要看新息序列（量测残差）是否接近白噪声。如果新息序列有明显的周期性或系统性偏置，说明模型动态与实际失配，光调滤波器参数救不回来。我习惯把新息序列画出来，跟量测噪声的标准差对比，如果新息方差远大于量测噪声方差，基本可以判定模型有问题，或者状态维度不够。

5. 常见问题与排错实录

5.1 矩阵维数和Cholesky分解报错

这是UKF新手最容易撞上的问题。Cholesky分解要求输入矩阵必须是实对称正定矩阵，但递推过程中(P)由于数值截断误差可能变成几乎奇异，甚至出现负特征值。典型报错是Matrix must be positive definite。

排查路径我总结过一套固定的顺序：先检查(P)的对称性，norm(P-P')是否在机器精度量级；再检查特征值，eig(P)是否全部大于0；最后看是不是初始化解就有问题，比如把P0设成一个非正定矩阵。大多情况下加一个小的对角矩阵就能绕过，但如果加了1e-12还是不能解决，那大概率是过程噪声(Q)设成了零矩阵，让协方差在预测步的山谷里被压没了。

5.2 估计曲线发散或长期不收敛

EKF和UKF都可能出现发散，但表现不同。EKF发散通常是曲线突然冲向某个极端值，然后再也回不来，原因是线性化误差在某个工作点附近被放大，协方差矩阵快速膨胀。UKF发散相对温和，更常见的是曲线长期振荡却收敛不到真实值附近，尤其是最初的采样点分布不当导致的。

现场排查我推荐一个实用动作：把滤波过程中的新息序列和卡尔曼增益一并画出来。如果卡尔曼增益在应该趋于稳定的时候还在剧烈波动，说明协方差更新可能出了问题；如果新息序列均值明显非零，说明状态偏差没有被修正掉，可能是量测模型或系统模型有系统性错误——比如量测方程里忘了乘以一个比例系数。

5.3 采样时间与量测噪声的不匹配

采样时间如果取得太大会导致离散化误差增大，滤波器估计的曲线会像锯齿一样抖动；取的太小则会导致两个连续量测之间的状态变化幅度小于数值噪声，滤波器变得“懒得更新”，表现是估计曲线过于平滑、对真实变化反应迟钝。量测噪声(R)的设定还要和采样时间匹配：PMU报告频率如果是几十帧每秒，对应的噪声标准差就要按该设备的实际指标来，不能随手填一个0.01。

如果后续准备用真实PMU数据验证代码，还需注意数据帧丢失和时间戳不同步的问题。仿真里我们假设每步量测都准时到达，但实际PMU数据流会有丢帧、乱序甚至同步错误，这时候直接把原始数据灌进滤波器，结果几乎必然会恶化。建议先在仿真环境加入缺失量测机制，测试滤波器在数据不完整时的表现，再过渡到真实数据。

5.4 代码版本兼容性问题

这套代码用到的核心函数（chol、矩阵运算、基本控制流）都是Matlab最基础的语法，从R2016a到R2024b都能正常运行。唯一可能出现兼容性警告的是二维数组索引和某些绘图函数，不涉及任何工具箱依赖。我见过有人引入Control System Toolbox或Simulink来实现，但完全没必要——这个项目用纯基础Matlab脚本就是最干净的做法。如果你用的是新版Matlab，打开脚本时自动转换旧语法后，照常运行即可。

6. 我的几点实操体会

项目做完之后我最大的感受是：EKF和UKF的差距，不是你调几个参数就能抹平的，但也不是所有场景下都值得用UKF。单机无穷大系统里两者的计算量差距可以忽略，但如果你想把这套方法推向实际大电网，几十台发电机、上百个状态变量时UKF的sigma点数量会直接翻倍，计算负担呈线性增长。那时候你可能要考虑基于稀疏矩阵的平方根UKF，或者干脆做混合策略——常规扰动用EKF，强非线性时段才切换UKF。

还有一个在代码之外但同样重要的经验：状态估计的验证不能只看自己仿真里的曲线，要看它对后续控制决策的影响。我后来把估算出的功角和角速度接入一个简单的紧急控制逻辑里，才发现同样的RMSE水平下，UKF的估计值在暂态过程中的“趋势跟随”能力明显更好，功角摆动的拐点抓得更准。这一点在RMSE指标上只有很小的差异，但对暂态稳定控制来说可能是天壤之别。

如果你刚拿到这个项目的代码，建议按这个顺序跑：先运行真实轨迹生成部分，画出功角摆动曲线，确认模型动态合理；再运行EKF，观察估计曲线是否贴住真实轨迹；最后切换UKF，对比误差指标。每一步都确认没毛病了再进行下一步，这样即使出问题也能快速定位。

后续还可以沿三个方向扩展：把发电机模型从二阶提升到四阶或六阶，加入励磁动态和PSS；将单机系统推广到多机系统，检验算法在高维状态下的数值表现；接入真实PMU数据，处理非理想量测条件。每个方向都有现成的公开数据和论文可参考，但起点永远是这一套能跑通、能改、能出图的Matlab框架。