企业官网建设流程全解析

电力系统动态状态估计，说具体点就是把电网里的关键状态量——比如发电机功角、转速、暂态电动势——在一连串带噪声的量测背景下实时“猜”出来。EKF（扩展卡尔曼滤波）和UKF（无迹卡尔曼滤波）是这个领域最常用的两种非线性滤波算法，也是很多人做Matlab仿真时的标配开局。最近我用Matlab把这两套滤波器的对比仿真完整跑了一遍，从状态方程搭建到滤波循环、再到调参踩坑，一路折腾下来收获不小，这里把思路、代码细节和踩坑经验都摊开讲讲。

这篇内容适合正在做电力系统动态状态估计课题的同学，以及想上手卡尔曼滤波、需要一个具体工程模型练手的读者。你只需要有一点基础：知道状态空间模型，能看懂矩阵运算，Matlab能跑通简单脚本就够了。我会用一套单机无穷大系统的算例贯穿全文，把EKF和UKF的原理、代码实现、参数整定、常见问题一次说清。

1. 电力系统动态状态估计：问题背景与核心思路

1.1 为什么电网需要“动态”状态估计

传统电力系统状态估计的主力是静态状态估计，基于加权最小二乘（WLS）做断面估计，解决的是“某一时刻的稳态潮流是多少”这类问题。但在故障后的暂态过程、新能源出力快速波动、负荷突变等时间尺度以毫秒到秒计的场景里，系统状态一直在变化，静态估计的刷新速度跟不上，而且它没有包含发电机转子运动等动态约束，无法外推未来轨迹。

动态状态估计的价值在于：用状态方程描述系统的动态演化规律，再用量测滤波不断修正状态，相当于“模型驱动加量测驱动”双轮定位。特别是在广域测量系统（WAMS）中，可以直接把PMU的高频量测喂给滤波器，实时输出功角和各动态变量的轨迹估计，为暂态稳定在线评估、自适应保护、紧急控制提供状态依据。这也是最近几年EKF、UKF以及更强的粒子滤波在电力系统方向频繁出现的原因。

1.2 动态状态估计的数学模型

动态状态估计用的框架和通用非线性滤波一致，状态方程与量测方程写成：

x(k+1) = f(x(k), u(k)) + w(k)，w ～ N(0, Q)
z(k) = h(x(k)) + v(k)，v ～ N(0, R)

在发电机动态估计里，状态量最标准的三要素是功角δ、角速度ω、暂态电动势E'q。u一般是机械功率Pm或励磁电压Efd，量测z来自PMU，通常是发电机机端电压相角、电压幅值、有功功率和无功功率。

f的具体形式来自发电机的转子运动方程和励磁绕组动态，本质上是连续时间微分方程，定好采样步长后用一阶欧拉离散就能转成上式里的f。h则把状态量映射到量测向量，电磁功率Pe = E'q V / X'd sinδ就是h里的一个典型分量。后面我会给出一个完整可跑的简化版模型，手把手带着写出来。

2. EKF与UKF原理拆解：两种截然不同的思路

2.1 EKF：用一阶泰勒展开对付非线性

EKF的基本假设是“非线性不太强时，一阶近似够用”。它对状态方程和量测方程分别在当前估计附近做泰勒展开，只保留一阶项，得到雅可比矩阵F和H，然后完全套用标准卡尔曼滤波的预测-更新流程。

预测步：

x_pred = f(x_hat(k-1))
P_pred = F P(k-1) F^T + Q

更新步：

K = P_pred H^T (H P_pred H^T + R)^(-1)
x_hat(k) = x_pred + K (z(k) - h(x_pred))
P(k) = (I - K H) P_pred

EKF实际用起来，最大的麻烦在于必须算雅可比矩阵。对电力系统这样耦合项较多的模型，手推雅可比容易漏项，尤其方程里有Pe = E'q V / X'd sinδ这种项，需要分别对三个状态量求导，还要仔细核对量纲。

EKF的优势也很明显：程序结构简单，每步只有一个矩阵求逆，实时性在同级别的非线性滤波里算优秀。弱点在于：系统进入强非线性区（比如故障瞬间功角快速摆开）时，一阶近似会带来明显截断误差，滤波轨迹可能滞后真实状态。

2.2 UKF：不线性化，直接算统计量

UKF换了一条路：不试图求导数，而是假设当前状态x服从高斯分布，按确定性的方式取2n+1个Sigma点（n为状态维数），让这些点逐个通过真实的非线性函数f和h传播，再对传播后的点做加权平均，算出下一步的均值和协方差。

这个过程本质上叫无迹变换（Unscented Transform），核心是用一组采样点去近似高斯随机变量经过非线性变换后的分布。通过调节α、β、κ三个参数，可以让均值精确到三阶、协方差精确到二阶（对高斯分布而言），精度上高于EKF的一阶展开。

Sigma点选完后，后面的结构和卡尔曼滤波几乎一样：预测均值、预测协方差、量测均值、量测协方差与互协方差，然后求增益矩阵并更新。可以理解成“用一群点替你算积分”的卡尔曼滤波器，这是它被称为“无迹”的原因——整个流程里根本不需要雅可比矩阵。

2.3 选型逻辑：具体场景怎么选

从我的实际对比来看，EKF和UKF在电力系统动态状态估计里没有谁绝对碾压谁。如果量测质量高、系统模型相对简单、实时性要求强，EKF是最稳的开局选择；如果系统进入强非线性段，或者状态维数达到五六维以上，UKF通常能把误差压得更低，代价是每个滤波周期要额外算2n+1次函数传播，计算量是EKF的2到3倍。

刚开始做这个方向的话，我建议两条路都写一遍。先用EKF把流程跑通，把模型和雅可比调对，再切到UKF做对比。这样既理解了两类算法的差异，也不会一上来就被Sigma点和权重参数搞晕。两种算法的主循环框架高度一致，切换成本并不高。

3. Matlab实现：模型搭建与滤波器代码

3.1 单机无穷大模型的状态方程与量测函数

我用的算例是一台发电机经线路接到无穷大母线，忽略励磁调节器动态，状态量取x = [δ; ω; E'q]，用经典三阶模型：

dδ/dt = ω - ωs
dω/dt = (Pm - Pe - D(ω - ωs)) / (2H)
dE'q/dt = (Efd - E'q - (Xd - X'd) V cosδ / X'd) / T'do

其中Pe = E'q V sinδ / X'd。

为了后面代码复用，我把参数放在结构体param里，这样函数引用方便，不容易串数据。连续动态函数可以写成：

function dx = gen_dyn(x, u, param) % x = [delta; omega; Eq] delta = x(1); omega = x(2); Eq = x(3); Pm = u(1); ws = 2*pi*50; Pe = Eq * param.V / param.Xdp * sin(delta); dx = zeros(3,1); dx(1) = omega - ws; dx(2) = (Pm - Pe - param.D * (omega - ws)) / (2 * param.H); dx(3) = (param.Efd - Eq - (param.Xd - param.Xdp) * param.V * cos(delta) / param.Xdp) ... / param.Tdo; end

离散化直接采用一阶欧拉：

function x_next = f_disc(x, u, param, dt) x_next = x + dt * gen_dyn(x, u, param); end

量测取两个量：发电机机端电压相角（以无穷大母线为参考，等于δ）和电磁功率Pe：

function z = h_func(x, param) delta = x(1); Eq = x(3); Pe = Eq * param.V / param.Xdp * sin(delta); z = [delta; Pe]; end

这里有个基础但非常关键的细节：角度单位全程用弧度。不要在代码里混用度数和弧度，否则状态方程里的sin相位和量测数值会错位，滤波基本必发散。另外D和H的单位也要统一，工程上常用标幺制，我这里直接用标幺值体系，转速偏差也用标幺量表示。

3.2 EKF的雅可比矩阵到底怎么求

EKF需要状态转移雅可比F和量测雅可比H。实际写代码时，我见过手推、符号求导、数值差分三种做法。我的建议是三种都过一遍，最终以数值差分为验证标准。

方法一，手推。对上面这个三阶模型，F的解析形式并不难写，但状态数一旦超过五六个，手推几乎必然出错。

方法二，用Matlab符号工具箱。把符号表达式定义好，用jacobian函数求导，再用matlabFunction转成数值函数句柄。这个方法适合一次性验证，缺点是符号表达式在复杂系统里很长，计算偏慢。

方法三，用中心差分数值法。这是我推荐写入工程代码里的方案：

function F = num_jac_f(x, u, param, dt) n = length(x); h = 1e-6; F = zeros(n, n); for i = 1:n xp = x; xm = x; xp(i) = x(i) + h; xm(i) = x(i) - h; F(:, i) = (f_disc(xp, u, param, dt) - f_disc(xm, u, param, dt)) / (2*h); end end

步长h从1e-6到1e-8都可以，太小会引入浮点误差，太大会让近似失真。我用1e-6对照符号求导的结果，误差通常在1e-10量级。量测雅可比H也用同样的差分方式求。写完之后，把解析结果和数值结果打印出来比一次，能对上，EKF大概率就没大问题了。

3.3 UKF的Sigma点与权重代码

UKF第一步是生成Sigma点。设当前状态估计为x_pred、协方差P_pred、状态维数n，以及缩放参数lambda（lambda = alpha^2 * (n + kappa) - n）：

P_root = chol((n + lambda) * P_pred, 'lower'); X_sig = zeros(n, 2*n+1); X_sig(:, 1) = x_pred; for i = 1:n X_sig(:, i+1) = x_pred + P_root(:, i); X_sig(:, n+i+1) = x_pred - P_root(:, i); end

权重按如下公式计算：

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) = 1 / (2*(n + lambda)); end

参数alpha控制Sigma点离均值点的远近，一般取1e-3到1之间；kappa一般取0；beta在高斯分布下取2最优。对于状态维数n=3的系统，如果取alpha=0.01、kappa=0，lambda会接近-2.9997，此时(n+lambda)非常小，chol分解矩阵接近半正定，数值上很容易报错。稳妥做法是不要把alpha取得太小，0.5到1之间更安全。网上很多默认模板不会提示这一点，但实际仿真里这就是chol崩掉的最主要原因。

Sigma点传播时，把X_sig每一列当成一个样本，分别调用f_disc，得到传播后的点集，再加权平均得到预测均值；每个点与预测均值的偏差用来合成预测协方差，最后加上过程噪声Q。量测更新同理：传播后的Sigma点再过h_func，算z_pred、P_zz和P_xz，最终算出增益矩阵K并更新。

3.4 完整滤波循环与两种算法的主流程代码

我把EKF和UKF封装成两个独立函数，返回值都是更新后的状态和协方差。主循环里先按真实动态生成带噪声的量测，再做滤波对比：

for k = 1:N % 真实轨迹与量测 x_true(:, k+1) = f_disc(x_true(:, k), u, param, dt) + sqrt(Q) * randn(n, 1); z(:, k) = h_func(x_true(:, k+1), param) + sqrt(R) * randn(m, 1); % EKF [x_ekf(:, k+1), P_ekf] = ekf_update(x_ekf(:, k), P_ekf, z(:, k), ... param, dt, Q, R); % UKF [x_ukf(:, k+1), P_ukf] = ukf_update(x_ukf(:, k), P_ukf, z(:, k), ... param, dt, Q, R, alpha, beta, kappa); end

EKF里我用数值差分求F和H，UKF里直接调用同一个f_disc和h_func，这样两组结果只反映算法差异，不掺杂模型和量测差异，对比才公平。指标上我习惯统计每个状态量滤波值与真实值的均方根误差（RMSE），这比单看一条曲线更客观。

4. 仿真算例与算法对比分析

4.1 算例参数与数据生成

算例主要参数我列成一张表，方便直接抄走用：

总仿真时长5秒，共250步。状态初值取[0.4; 2pi50; 1.0]，在t=2s时给Pm加上0.2的阶跃扰动，模拟负荷突变场景。噪声矩阵设置为：

Q = diag([1e-6, 1e-4, 1e-6])
R = diag([1e-4, 1e-4])
P0 = diag([1e-3, 1e-2, 1e-3])

真实轨迹生成时，用ODE风格的高精度积分也不过分，但对教学算例，欧拉法加小步长已经足够。若要更严谨，可以用ode45生成参考轨迹，再按0.02s重采样，这一步可以根据个人习惯来。

4.2 跟踪轨迹与误差对比

在上述扰动工况下，EKF和UKF都能收敛并跟踪真实轨迹。但看误差细节，EKF在扰动发生后的头几个周期有明显超调，功角RMSE大概在0.0032 rad附近；UKF的RMSE大约只有EKF的一半，功角误差在0.0017 rad左右。计算耗时上，EKF在我的机器上每步约0.4ms，UKF约1.2ms，两者的实时性都绰绰有余。

这个结果符合理论预期：这个三阶模型维数不高，非线性主要体现在Pe项上，扰动瞬间δ变化快，一阶线性化的截断误差被放大。UKF通过Sigma点传播，相当于“多测几个点再平均”，对非线性段的适应力更强。稳态阶段两者差距不大，说明EKF在缓变区间完全够用。

4.3 扩展到多机系统时的计算量变化

多机系统的状态维数会升到几十甚至上百维。EKF需要构造n×n状态雅可比，UKF每次要传播2n+1个点，两者的开销都会增长，但UKF增幅更明显。实际工程里，可以把全网按“感兴趣区域加外部等值”的思路分解，只对关键发电机节点跑动态估计，或者对弱非线性区域用EKF、强非线性区域用UKF混合处理。我在仿真阶段没有做得这么细，但这是一个值得延伸的方向。

多机系统还有一个需要注意的问题：可观测性。状态维数上去之后，量测数量不够时滤波可能收敛到错误轨迹。建议提前做一个简单的可观测性检查，用线性化模型看看秩是否够，或者加虚拟量测约束估计范围。

5. 常见问题与调试经验实录

5.1 滤波发散：原因定位和应急手段

EKF和UKF在实际仿真里最让人头疼的就是发散。先给一张原因速查表：

定位发散的第一步，是做一个“理想量测”实验：把R设到极小量级，比如1e-10，看滤波值能不能稳定贴住真实轨迹。能贴住，说明模型和算法主体没问题，问题就在噪声参数或量测匹配上；贴不住，优先查f_disc、h_func和雅可比。

应急手段也有两个实用技巧。一是给预测协方差乘一个大于1的衰减因子，比如P_pred = 1.02 * P_pred，相当于给当前量测更高权重，对强非线性暂态段有效。二是引入渐消记忆，把Q按指数形式放大，让滤波器主动遗忘过去。

5.2 Q、R、P0三个矩阵怎么调

调参是我这次仿真里耗时最长的环节。先说R。R通常由量测噪声标准差直接换算。PMU相角测量典型误差在0.01到0.1度量级，换算成弧度后平方，就是R对应元素。先按这个方式定R，不要再额外多调，否则容易过拟合。

P0本质上是“对初始状态的信任程度”。如果初值来自潮流计算或静态状态估计，不会差太远，P0取状态量量级的千分之一到万分之一即可，设太大没有意义，反而会让前几步滤波剧烈摆动。

Q代表模型误差。模型越粗糙、近似越多，Q就应越大；反过来模型足够准，Q可以很小。我一般让Q对角线取P0的1/100到1/10之间起步，再根据轨迹的“粘滞感”微调。滤波值明显滞后真实轨迹，就把Q调大一档；滤波值高频抖动，就把Q调小一档。

调参时建议做一个记录表，把每次的P0、Q、R和对应RMSE记下来。这个操作听起来土，但非常有用。我试过连续调十组，最后发现某几组不同参数组合的RMSE几乎一样，说明在这个算例里滤波对这个参数不敏感，这本身就是很宝贵的判断信息。

5.3 Matlab编码中的四个坑

第一是chol分解失败。协方差矩阵非正定时无法分解。规避方法：打印P的特征值，如果出现负数或非常接近0的数，给P对角线加1e-8级别的正则化项。我代码里专门写了一个辅助函数来检查并修复P的正定性。

第二是角度单位混用。δ在状态方程和量测函数里必须是弧度，但画图时很多人想用度数。我处理的办法是画图时单独转换，计算链路里始终用弧度，避免来回换算导致误差。

第三是randn与sqrt(Q)维度不匹配。状态维数变了，但噪声向量维度还是旧的，这是改模型时最隐蔽的bug。建议在仿真开头用assert强制检查维度：

assert(isequal(size(Q), [n n]), 'Q维度与状态维数不匹配'); assert(isequal(size(R), [m m]), 'R维度与量测维数不匹配');

第四是P对称性丢失。经过多轮更新后，数值误差会让P越来越不对称，最终导致chol崩掉。我在每轮滤波循环末尾都会加一句：

P = 0.5 * (P + P');

这个小操作成本极低，但对数值稳定性帮助非常大。

排查时还有一种高效路径：先故意构造一个全知场景，让真实轨迹从滤波器初值附近出发，且量测噪声设得很小，这种情况下滤波理论上必然收敛。如果连这个场景都不收敛，说明是代码逻辑错误；能收敛，再逐步加大噪声、加扰动，问题来源会被快速隔离。

最后分享一点个人体会。这类仿真里，最耗时间的往往不是滤波算法本身，而是把状态方程、量测方程、雅可比、参数体系调试到四者完全自洽。模型是对的，EKF和UKF都能给出可用结果；模型只要有错，再好的滤波器也是垃圾进垃圾出。所以我强烈建议，拿到一个课题先用EKF把全链路打通，再上UKF做对比，这个顺序会省掉一大半调参的苦功夫。

后续如果想继续做深，我建议往三个方向延伸：一是用平方根UKF提升数值稳定性，彻底绕开chol分解问题；二是在量测里混入少量坏数据，加卡方检验或鲁棒处理，模拟PMU受干扰的真实环境；三是把模型从单机扩展到多机系统，那时候状态维数上来，滤波器的收敛性差异和计算负担差异才会真正体现出来。有条件的同学不妨从今天这个单机版本开始改，上手难度会低很多。