做电力系统动态状态估计这个方向我前前后后折腾了一个多月从一开始对着 EKF 公式发懵到后来把 UKF 也跑通再到在 Matlab 里一套流程完整复现中间踩了不少坑也积累了不少能直接用的经验。这篇博文就围绕“基于扩展卡尔曼滤波(EKF)和无迹卡尔曼滤波(UKF) 的电力系统动态状态估计Matlab 代码实现”这个项目展开把从模型搭建、算法推导、代码实现到参数调试的完整过程都讲清楚。如果你正准备做电力系统动态状态估计或者刚接触卡尔曼滤波、想在 Matlab 里跑通一个可复现的实验这篇内容应该能帮你省下不少走弯路的时间。1. 先把问题讲明白电力系统动态状态估计到底在估什么1.1 为什么静态状态估计不够用传统电力系统状态估计绝大多数是基于 SCADA 量测的加权最小二乘WLS静态状态估计。它的思路是拿到当前断面的遥测数据然后解一个优化问题得到系统各节点电压幅值和相角的“最优估计值”。这个方法在能量管理系统中运行了非常多年但它本质上只能给出一个静态断面的状态没办法反映系统的动态演化过程。随着广域测量系统WAMS和同步相量测量装置PMU的普及量测数据的刷新率从秒级提高到了毫秒级。我们手里有了高速率、带时标的同步相量数据如果还只做静态估计就浪费了这些数据的时间维信息。更重要的是在暂态过程、低频振荡、电压失稳等动态场景下调度和控制系统需要的是发电机功角、转速、暂态电动势等状态量的实时轨迹而不仅仅是稳态断面。这就催生了动态状态估计Dynamic State Estimation, DSE的需求。动态状态估计的核心是把电力系统动态方程通常是发电机转子运动方程和电磁暂态方程作为状态转移模型把 PMU/SCADA 量测作为观测模型然后用非线性滤波算法从带噪声的量测中递推估计出系统状态随时间的变化。和静态估计相比DSE 相当于在时间轴上做了一次“跟踪”能提前发现状态的趋势性变化这也是它在暂态稳定预测和在线安全评估中备受关注的原因。1.2 动态状态估计的数学模型怎么搭要做动态状态估计第一步是把电力系统的动态过程写成状态空间形式。标准形式如下[ \begin{aligned} x_{k1} f(x_k, u_k) w_k \ z_k h(x_k) v_k \end{aligned} ]其中(x_k) 是第 (k) 个采样时刻的状态向量(u_k) 是输入向量(z_k) 是量测向量(w_k) 和 (v_k) 分别是过程噪声和量测噪声通常假设它们服从零均值高斯分布协方差矩阵分别为 (Q) 和 (R)。具体到电力系统状态量通常选取发电机动态状态比如功角 (\delta)、转速偏差 (\omega)、暂态电动势 (E_q)、(E_d) 等。以最常用的三阶单机模型为例连续时间动态方程可以写成[ \begin{aligned} \frac{d\delta}{dt} \omega - \omega_s \ \frac{d\omega}{dt} \frac{\omega_s}{2H}(P_m - P_e - K_D(\omega - \omega_s)) \ \frac{dEq}{dt} \frac{1}{T{d0}}(-Eq E{fd} - (x_d - x_d)I_d) \end{aligned} ]这里 (\delta) 是发电机功角(\omega) 是电角速度(H) 是惯性时间常数(P_m)、(P_e) 分别是机械功率和电磁功率(T_{d0}) 是励磁绕组时间常数(x_d)、(x_d) 是同步电抗和暂态电抗。如果要做更精细的建模还可以扩展到四阶、五阶甚至六阶模型阶数越高状态量越多滤波器的计算负担也越大。在实际 Matlab 仿真中我不会直接用连续时间模型而是先对它做离散化。最粗暴的方式是欧拉法[ x_{k1} x_k T_s \cdot f(x_k, u_k) ]其中 (T_s) 是采样周期。如果模型非线性程度比较强我会用四阶龙格库塔RK4做离散化精度更高代码也写不了几行。后面第 3 章会给出具体代码实现。1.3 EKF 和 UKF 怎么选扩展卡尔曼滤波EKF和无迹卡尔曼滤波UKF是两种最常见的非线性滤波算法。EKF 的思路很简单在每一个时刻把非线性函数 (f) 和 (h) 在当前状态附近做一阶泰勒展开用雅可比矩阵代替线性系统中的状态转移矩阵和观测矩阵然后套用标准卡尔曼滤波的预测—更新框架。EKF 的好处是实现简单、计算量小但缺点是线性化误差在强非线性场景下会被放大而且雅可比矩阵的推导容易出错。UKF 则走了一条完全不同的路子它不对方程做线性化而是通过无迹变换Unscented Transform, UT选取一组称为 sigma 点的采样点让这些点经过非线性函数传播再根据传播后的点集计算均值和协方差。因为 sigma 点能捕捉到非线性函数的高阶统计信息UKF 在强非线性下的估计精度通常更高而且不需要求雅可比矩阵。从工程实现角度两者最主要的差别可以归纳为下面这个表对比项EKFUKF核心思想一阶泰勒展开线性化sigma 点无迹变换是否需要雅可比矩阵需要推导烦琐不需要计算量较小略大但可接受强非线性场景精度一般线性化误差明显较高能捕捉二阶以上信息实现难度中低公式简洁中需要处理权重和采样逻辑参数敏感性对 Q、R 敏感对 Q、R 和 sigma 参数敏感我个人的习惯是如果是刚入门、先想跑通一个简单模型用 EKF 打底把流程搞清楚等模型非线性程度变高或者精度不够时再切换到 UKF。两种滤波器在 Matlab 里都可以用矩阵运算批量实现代码量差异不大后面我会把两者的实现骨架都贴出来。2. EKF 的推导与实现雅可比矩阵才是真正的坎2.1 EKF 的预测—更新流程EKF 的基本流程和线性卡尔曼滤波很相似只是把状态转移矩阵 (F) 和观测矩阵 (H) 换成了雅可比矩阵。预测步[ \begin{aligned} \hat{x}{k1|k} f(\hat{x}{k|k}) \ P_{k1|k} F_k P_{k|k} F_k^T Q_k \end{aligned} ]更新步[ \begin{aligned} K_{k1} P_{k1|k} H_{k1}^T (H_{k1} P_{k1|k} H_{k1}^T R_{k1})^{-1} \ \hat{x}{k1|k1} \hat{x}{k1|k} K_{k1}(z_{k1} - h(\hat{x}{k1|k})) \ P{k1|k1} (I - K_{k1} H_{k1}) P_{k1|k} \end{aligned} ]其中 (F_k \frac{\partial f}{\partial x}\big|{\hat{x}{k|k}})(H_{k1} \frac{\partial h}{\partial x}\big|{\hat{x}{k1|k}})。用开车导航来类比这个过程特别好理解。预测步就是“你根据上一秒的车速和方向猜测这一秒大概开到了哪里”协方差矩阵代表你对自己猜测有多大的把握——猜测越久、把握越低所以 (P) 会变大。更新步则是“GPS 报了一个位置同时你瞄了一眼路边的参照物综合两者获得一个更精确的位置”。卡尔曼增益 (K) 就是用来权衡“我预测的结果”和“传感器测量结果”各自可信度的一个系数。2.2 雅可比矩阵怎么求才不容易错EKF 最麻烦的一步就是求雅可比矩阵。很多人在推导阶段就卡住了矩阵维度一多就晕。我自己的做法分三种情况第一种情况模型简单我直接手推。比如状态向量是 (\begin{bmatrix}\delta \omega E_q\end{bmatrix}^T)观测方程是 (h(x) \begin{bmatrix}\delta \omega\end{bmatrix}^T)那么观测矩阵 (H) 就是[ H \begin{bmatrix} 1 0 0 \ 0 1 0 \end{bmatrix} ]这个很简单因为观测就是直接测量状态量。第二种情况模型复杂我直接用 Matlab 的 Symbolic Math Toolbox 求偏导然后把符号表达式转换成可调用的函数句柄。比如syms delta omega Eq real syms Pm Pe KD H ws Ts real % 状态方程 f1 omega - ws; f2 ws / (2*H) * (Pm - Pe - KD*(omega - ws)); f3 (-Eq 1) / (Td0); % 占位示例实际按模型写 F_sym jacobian([f1; f2; f3], [delta, omega, Eq]);用jacobian求出矩阵后再通过matlabFunction转成数值函数放进 EKF 循环里反复调用。这样做的好处是公式修改后不用手推但注意符号运算在循环里不能频繁调用否则速度会慢到怀疑人生。正确做法是提前生成好函数再在滤波循环外一次性编译。第三种情况实在不想求导就直接上 UKF绕开雅可比。这也是我在项目后期切换 UKF 的最大动力——模型单改一个参数EKF 的雅可比就得跟着重新推UKF 完全不用管。2.3 离散化与噪声协方差的一个隐蔽细节EKF 实现里有一个特别容易忽略的问题我们从物理模型写出的是连续时间状态方程但数字实现是离散时间的。最简单的做法是直接欧拉离散状态转移函数变为[ f_d(x_k) x_k T_s f(x_k) ]但过程噪声协方差 (Q) 也需要从连续形式转换过来。连续系统的过程噪声协方差是 (Q_c)离散化后的近似关系是 (Q_d \approx Q_c \cdot T_s)。如果直接沿用连续时间下的 (Q)滤波器的增益会偏大或偏小导致估计轨迹出现奇怪的高频抖动。在我做的算例里系统采样周期 (T_s 0.01s)连续过程噪声标准差设成 (0.001)那么离散 (Q_d) 的对角元素我一开始直接用了 (0.001^2)结果滤波曲线抖得厉害。后来把 (Q_d Q_c \cdot T_s) 换成 (10^{-8}) 量级曲线立刻平滑了很多。这种细节不写进教材但实际调试过程中遇到的概率极高。3. UKF 的实现不取导数也能传播协方差3.1 UT 变换到底做了什么UKF 的核心是 UT 变换。它的想法非常直观既然直接求非线性函数的均值和协方差很困难那我就巧妙地选几个点让这些点经过非线性函数传播后能还原出足够精确的统计特性。具体做法是对于 (n) 维状态向量生成 (2n1) 个 sigma 点[ \begin{aligned} \chi_0 \bar{x} \ \chi_i \bar{x} (\sqrt{(n\lambda)P})i, \quad i1,\dots,n \ \chi{in} \bar{x} - (\sqrt{(n\lambda)P})_i, \quad i1,\dots,n \end{aligned} ]其中 ((\sqrt{(n\lambda)P})_i) 表示矩阵 ((n\lambda)P) 的 Cholesky 分解后第 (i) 列。每个 sigma 点对应的权重为[ \begin{aligned} W_0^{(m)} \frac{\lambda}{n\lambda} \ W_0^{(c)} \frac{\lambda}{n\lambda} (1-\alpha^2\beta) \ W_i^{(m)} W_i^{(c)} \frac{1}{2(n\lambda)}, \quad i1,\dots,2n \end{aligned} ]参数 (\lambda \alpha^2(n\kappa) - n)。其中 (\alpha) 决定 sigma 点的散布程度一般取 (1e-3) 到 (1) 之间(\beta) 用于融入先验分布信息高斯分布下取 2 最优(\kappa) 是次级缩放参数通常取 0 或 (3-n)。把这些 sigma 点逐个扔进状态转移函数 (f(\cdot))得到传播后的点集然后加权计算预测均值和协方差。量测更新同理。整个过程完全没有求导运算所以对非线性函数形式没有任何平滑性要求这也是 UKF 对“硬非线性”更友好的原因。3.2 sigma 点参数怎么给在 Matlab 里做 UKF参数取值是决定成败的关键。我试过几组典型参数经验如下(\alpha 1e-3)sigma 点离均值很近适合状态量变化平缓的场景。(\alpha 1)sigma 点散布范围大状态突变时跟踪更快但稳态噪声可能偏大。(\beta 2)高斯分布下最优。(\kappa 0) 或 (3-n)如果 (n) 较大(\kappa) 取负数可能导致协方差矩阵失去半正定性这时我会改用 (\kappa0)。有一个常见坑计算 ((n\lambda)P) 的时候如果 (n\lambda \le 0)Cholesky 分解会直接失败。所以实现里必须先检查 (n\lambda 0)否则要给 (\kappa) 加一个正数偏移。这个报错在 Matlab 里通常会以Matrix must be positive definite的形式出现我一开始遇到还以为是自己协方差更新公式写错了查了半天才发现是参数问题。3.3 UKF 的 Matlab 骨架代码下面这段代码是 UKF 的核心循环骨架适合直接改成自己的模型function [x_est, P_est] ukf_predict_update(x_est, P_est, z, Q, R, f, h, param) n length(x_est); alpha param.alpha; beta param.beta; kappa param.kappa; lambda alpha^2 * (n kappa) - n; % 权重 Wm [lambda/(nlambda), 0.5/(nlambda)*ones(1, 2*n)]; Wc [lambda/(nlambda) (1-alpha^2beta), 0.5/(nlambda)*ones(1, 2*n)]; % 生成 sigma 点 A chol((nlambda) * P_est, lower); chi zeros(n, 2*n1); chi(:,1) x_est; for i 1:n chi(:,i1) x_est A(:,i); chi(:,ni1) x_est - A(:,i); end % 预测步 chi_pred zeros(n, 2*n1); for i 1:2*n1 chi_pred(:,i) f(chi(:,i)); end x_pred chi_pred * Wm; P_pred Q; for i 1:2*n1 d chi_pred(:,i) - x_pred; P_pred P_pred Wc(i) * (d * d); end % 量测更新 Z_pred zeros(size(z,1), 2*n1); for i 1:2*n1 Z_pred(:,i) h(chi_pred(:,i)); end z_pred Z_pred * Wm; Pzz R; Pxz zeros(n, size(z,1)); for i 1:2*n1 dz Z_pred(:,i) - z_pred; dx chi_pred(:,i) - x_pred; Pzz Pzz Wc(i) * (dz * dz); Pxz Pxz Wc(i) * (dx * dz); end K Pxz / Pzz; x_est x_pred K * (z - z_pred); P_est P_pred - K * Pzz * K; end这段代码里有一点需要注意P_est P_pred - K * Pzz * K是 Joseph 形式的等价写法比(I - KH)P_pred数值稳定性更好能尽量避免协方差矩阵负定。实际使用时f和h要定义成函数句柄形如f (x) x Ts * model_dyn(x, u);。如果模型状态维数不高、循环次数不多这种逐点传播的写法够用如果状态维数和采样点都很大可以考虑用矩阵批量传播来加速但代码可读性会下降。4. Matlab 实操从模型搭建到仿真出图的全流程4.1 仿真系统怎么搭我的做法是在 Matlab 里自己构造一个“数字孪生”系统先选定一个发电机模型用已知的状态初值和输入序列通过状态方程生成一条“真实状态轨迹”然后叠加高斯白噪声生成模拟量测。这样做的好处是真实状态已知可以计算估计误差并验证滤波器性能。这里我用一个简化两阶模型做演示状态量取功角 (\delta) 和转速偏差 (\omega)。状态方程为[ \begin{aligned} \delta_{k1} \delta_k T_s(\omega_k - \omega_s) \ \omega_{k1} \omega_k T_s \cdot \frac{\omega_s}{2H}(P_m - P_e(\delta_k) - K_D(\omega_k - \omega_s)) \end{aligned} ]其中电磁功率 (P_e(\delta_k) EV/x \sin(\delta_k))这是一个强非线性环节量测设置为直接测量 (\delta) 和 (\omega)并加入方差 (R) 的高斯噪声。系统参数大体如下参数数值说明(\omega_s)2pi50同步角速度H5.0 s惯性时间常数KD2.0阻尼系数Ts0.01 s采样周期/滤波步长R_diag(0.01)^2量测噪声方差Q_diag(0.001)^2过程噪声方差4.2 主程序流程整个主程序按照“生成真实轨迹 → 生成量测 → 初始化滤波器 → 滤波循环 → 误差分析”的顺序组织。核心循环部分如下% 参数初始化 n 2; % 状态维数 x_true zeros(n, N); x_true(:,1) x_init; z_meas zeros(n, N); z_meas(:,1) h(x_init) sqrt(R)*randn(n,1); % 生成真实轨迹和量测 for k 1:N-1 x_true(:,k1) f_dt(x_true(:,k), u, Ts); z_meas(:,k1) h(x_true(:,k1)) sqrt(R)*randn(n,1); end % 滤波器初始化 x_est zeros(n, N); x_est(:,1) x_init 0.1*randn(n,1); P_est eye(n); % 滤波循环 for k 1:N-1 % 预测 [x_pred, P_pred] ukf_predict(x_est(:,k), P_est, (x) f_dt(x,u,Ts), Q, param); % 更新 [x_est(:,k1), P_est] ukf_update(x_pred, P_pred, z_meas(:,k1), (x) h(x), R, param); end % 误差计算 rmse_delta sqrt(mean((x_est(1,:) - x_true(1,:)).^2));这里f_dt是离散化后的状态转移函数h是观测函数。EKF 版本只需要把ukf_predict和ukf_update换成对应的 EKF 函数即可主程序结构完全一样。4.3 Q、R 怎么调才能出好效果在动态状态估计里Q 和 R 的设置直接决定滤波效果。我的方法分三步第一步R 从量测的物理噪声水平出发。PMU 的幅值测量误差典型范围在 0.1% 到 1% 之间角度误差在 0.01 到 0.1 弧度级。据此设置 R 的对角元素这是有物理依据的不建议为了效果盲目调小。第二步Q 先给一个小值然后根据观测残差调整。如果滤波曲线滞后于真实轨迹说明 Q 太小——模型预测的“信任度”太高滤波器不跟测量走如果滤波曲线抖动剧烈说明 Q 太大——噪声被当成真实变化滤波器被测量噪声牵着鼻子走。第三步用 RMSE 做定量评估配合几个典型工况反复试。我常用的经验公式是先设 (Q q \cdot I)跑一次得到 RMSE然后 (q) 按 10 倍步长上下扫描看 RMSE 的谷值出现在哪个量级。这个方法虽然朴素但比凭感觉调参可靠得多。4.4 结果图怎么画才直观滤波做完后画图是必不可少的。我最常用的画图组合是第一张图把真实状态、带噪声量测、EKF估计值、UKF估计值四条曲线画在一起直观对比跟踪效果第二张图画估计误差并画出正负 3 倍标准差的置信区间第三张图画协方差矩阵对角线元素的变化趋势观察滤波器是否收敛。画图代码大概长这样figure; subplot(2,1,1); plot(t, x_true(1,:), k-, LineWidth, 1.5); hold on; plot(t, z_meas(1,:), b., MarkerSize, 3); plot(t, x_est_ekf(1,:), r--, LineWidth, 1.2); plot(t, x_est_ukf(1,:), g-., LineWidth, 1.2); legend(真实值, 量测值, EKF, UKF); xlabel(时间 (s)); ylabel(功角 (rad)); grid on;为什么要把量测值也画上去因为这样可以直观看出滤波到底“平滑”了多少。如果滤波结果和量测噪声几乎重合那说明滤波增益过大相当于没滤波如果滤波结果和量测偏差太大说明增益过小跟踪跟不上。5. 常见问题与排错实录5.1 滤波器发散怎么排查我遇到过的最典型的发散现象是估计值在第一、第二步就被量测噪声带飞之后再也回不到真实轨迹附近协方差矩阵也快速膨胀甚至变成 NaN。排查顺序我一般按下面这个清单症状可能原因对策估计值完全偏离真实值初值设置离真实值太远用静态估计结果或潮流解初始化前几步就出现 NaN协方差矩阵非正定用 Joseph 形式更新检查 sigma 点参数稳态时误差仍然很大Q 过大或 R 过小减小 Q、增大 R用残差统计调整响应太慢、滞后明显Q 过小或 R 过大增大 Q、减小 R轨迹跟随量测噪声高频抖动Q 过大或 R 过小检查 Q 是否忘了乘采样周期5.2 状态量纲不一致导致的数值问题功角的量纲是弧度转速偏差可能取极标幺值也可能取弧度每秒两者数值量级可能差了 10^3 以上。如果不做归一化协方差矩阵里不同对角元的数量级差太多Cholesky 分解很容易数值不稳定。我后来把所有状态量都做了归一化处理功角用弧度转速偏差用标幺值 ( \Delta\omega / \omega_s )电磁功率用标幺值。全部统一到 0.1 到 10 这个量级范围滤波稳定性和精度都好了很多。这也提醒我卡尔曼滤波虽然对量纲不敏感但数值实现里量级差距过大会引发各种隐蔽问题。5.3 初值到底该怎么给动态状态估计的初值不能随便给。EKF 和 UKF 本质上都是局部递推算法初值离真实状态太远第一轮更新就可能让协方差矩阵失去意义。工程上常用的做法是先用一个常规静态状态估计器得到初始断面或者直接用潮流计算结果初始化。如果实在没有可以先把滤波器放在“稳态工况”下预热跑一段时间用滤波结果反推初值。我在单机模型算例里试过初值偏差在 5% 以内时滤波器都能快速收敛偏差超过 20% 时UKF 还能勉强拉回来EKF 就已经发散得没法看了。5.4 计算效率能不能再快一点Matlab 里做 UKF 容易陷入双重循环导致仿真时长暴增。我的优化策略有两个第一sigma 点传播尽量用矩阵运算代替循环。把 2n1 个 sigma 点拼成一个 (n, 2n1) 的矩阵然后用“一次函数求值”的方式处理整个矩阵但前提是f和h能向量化。第二避免在滤波循环内重复分配大数组。提前把chi,chi_pred,Z_pred等变量初始化好循环内只做赋值。第三个策略是能离线算的都不在线算。比如状态转移函数的某些固定结构、雅可比矩阵的符号表达式都提前算好。之前我犯过在 EKF 循环里每次都调matlabFunction的错那个速度之慢跑一次仿真够我喝三杯咖啡。6. 一次完整的调试经历最后分享一个具体的调试案例。当时我在单机无穷大系统上做 EKF参数全部按论文里的标准值设定结果滤波曲线在第 1 秒就开始剧烈振荡然后直接发散。第一反应是雅可比矩阵求错了。我重新用符号工具求了一次偏导逐项对比没发现错。后来我把过程噪声 (Q) 的对角元素从 (10^{-4}) 一步步往下调调到 (10^{-8}) 之后滤波曲线才稳定下来。原因在于我的状态转移方程里转速偏差的变化率数值本身就很小过程噪声给大了状态预测的协方差膨胀很快导致滤波器过度信任量测而量测噪声又存在于是估计值被量测噪声完全主导。这个教训让我彻底记住了Q 和 R 的绝对数值不重要重要的是它们和状态方程中各项系数的相对关系。这也是为什么很多论文里直接给一个 Q、R 矩阵你照搬过来却跑不出同样效果的原因——你的模型参数、单位、离散化方式都和作者不一样。后来我在同一套模型上把 UKF 也跑通对比 RMSE 发现 UKF 比 EKF 的估计误差平均小 30% 左右尤其在功角突变阶段。虽然 UKF 的计算时间大约是 EKF 的 1.8 倍但在这个算例规模下完全可接受。综合精度和调参成本我现在更推荐在电力系统动态状态估计中优先尝试 UKF。如果后续你有兴趣继续深入可以考虑把量测扩展为“伪量测PMU 混合量测”或者在 EKF/UKF 基础上加入鲁棒性设计应对量测异常和参数不确定性。这些都是很有价值的扩展方向但前提是把基础滤波流程彻底跑通。 基于EKF与UKF的电力系统动态状态估计Matlab实现
做电力系统动态状态估计这个方向我前前后后折腾了一个多月从一开始对着 EKF 公式发懵到后来把 UKF 也跑通再到在 Matlab 里一套流程完整复现中间踩了不少坑也积累了不少能直接用的经验。这篇博文就围绕“基于扩展卡尔曼滤波(EKF)和无迹卡尔曼滤波(UKF) 的电力系统动态状态估计Matlab 代码实现”这个项目展开把从模型搭建、算法推导、代码实现到参数调试的完整过程都讲清楚。如果你正准备做电力系统动态状态估计或者刚接触卡尔曼滤波、想在 Matlab 里跑通一个可复现的实验这篇内容应该能帮你省下不少走弯路的时间。1. 先把问题讲明白电力系统动态状态估计到底在估什么1.1 为什么静态状态估计不够用传统电力系统状态估计绝大多数是基于 SCADA 量测的加权最小二乘WLS静态状态估计。它的思路是拿到当前断面的遥测数据然后解一个优化问题得到系统各节点电压幅值和相角的“最优估计值”。这个方法在能量管理系统中运行了非常多年但它本质上只能给出一个静态断面的状态没办法反映系统的动态演化过程。随着广域测量系统WAMS和同步相量测量装置PMU的普及量测数据的刷新率从秒级提高到了毫秒级。我们手里有了高速率、带时标的同步相量数据如果还只做静态估计就浪费了这些数据的时间维信息。更重要的是在暂态过程、低频振荡、电压失稳等动态场景下调度和控制系统需要的是发电机功角、转速、暂态电动势等状态量的实时轨迹而不仅仅是稳态断面。这就催生了动态状态估计Dynamic State Estimation, DSE的需求。动态状态估计的核心是把电力系统动态方程通常是发电机转子运动方程和电磁暂态方程作为状态转移模型把 PMU/SCADA 量测作为观测模型然后用非线性滤波算法从带噪声的量测中递推估计出系统状态随时间的变化。和静态估计相比DSE 相当于在时间轴上做了一次“跟踪”能提前发现状态的趋势性变化这也是它在暂态稳定预测和在线安全评估中备受关注的原因。1.2 动态状态估计的数学模型怎么搭要做动态状态估计第一步是把电力系统的动态过程写成状态空间形式。标准形式如下[ \begin{aligned} x_{k1} f(x_k, u_k) w_k \ z_k h(x_k) v_k \end{aligned} ]其中(x_k) 是第 (k) 个采样时刻的状态向量(u_k) 是输入向量(z_k) 是量测向量(w_k) 和 (v_k) 分别是过程噪声和量测噪声通常假设它们服从零均值高斯分布协方差矩阵分别为 (Q) 和 (R)。具体到电力系统状态量通常选取发电机动态状态比如功角 (\delta)、转速偏差 (\omega)、暂态电动势 (E_q)、(E_d) 等。以最常用的三阶单机模型为例连续时间动态方程可以写成[ \begin{aligned} \frac{d\delta}{dt} \omega - \omega_s \ \frac{d\omega}{dt} \frac{\omega_s}{2H}(P_m - P_e - K_D(\omega - \omega_s)) \ \frac{dEq}{dt} \frac{1}{T{d0}}(-Eq E{fd} - (x_d - x_d)I_d) \end{aligned} ]这里 (\delta) 是发电机功角(\omega) 是电角速度(H) 是惯性时间常数(P_m)、(P_e) 分别是机械功率和电磁功率(T_{d0}) 是励磁绕组时间常数(x_d)、(x_d) 是同步电抗和暂态电抗。如果要做更精细的建模还可以扩展到四阶、五阶甚至六阶模型阶数越高状态量越多滤波器的计算负担也越大。在实际 Matlab 仿真中我不会直接用连续时间模型而是先对它做离散化。最粗暴的方式是欧拉法[ x_{k1} x_k T_s \cdot f(x_k, u_k) ]其中 (T_s) 是采样周期。如果模型非线性程度比较强我会用四阶龙格库塔RK4做离散化精度更高代码也写不了几行。后面第 3 章会给出具体代码实现。1.3 EKF 和 UKF 怎么选扩展卡尔曼滤波EKF和无迹卡尔曼滤波UKF是两种最常见的非线性滤波算法。EKF 的思路很简单在每一个时刻把非线性函数 (f) 和 (h) 在当前状态附近做一阶泰勒展开用雅可比矩阵代替线性系统中的状态转移矩阵和观测矩阵然后套用标准卡尔曼滤波的预测—更新框架。EKF 的好处是实现简单、计算量小但缺点是线性化误差在强非线性场景下会被放大而且雅可比矩阵的推导容易出错。UKF 则走了一条完全不同的路子它不对方程做线性化而是通过无迹变换Unscented Transform, UT选取一组称为 sigma 点的采样点让这些点经过非线性函数传播再根据传播后的点集计算均值和协方差。因为 sigma 点能捕捉到非线性函数的高阶统计信息UKF 在强非线性下的估计精度通常更高而且不需要求雅可比矩阵。从工程实现角度两者最主要的差别可以归纳为下面这个表对比项EKFUKF核心思想一阶泰勒展开线性化sigma 点无迹变换是否需要雅可比矩阵需要推导烦琐不需要计算量较小略大但可接受强非线性场景精度一般线性化误差明显较高能捕捉二阶以上信息实现难度中低公式简洁中需要处理权重和采样逻辑参数敏感性对 Q、R 敏感对 Q、R 和 sigma 参数敏感我个人的习惯是如果是刚入门、先想跑通一个简单模型用 EKF 打底把流程搞清楚等模型非线性程度变高或者精度不够时再切换到 UKF。两种滤波器在 Matlab 里都可以用矩阵运算批量实现代码量差异不大后面我会把两者的实现骨架都贴出来。2. EKF 的推导与实现雅可比矩阵才是真正的坎2.1 EKF 的预测—更新流程EKF 的基本流程和线性卡尔曼滤波很相似只是把状态转移矩阵 (F) 和观测矩阵 (H) 换成了雅可比矩阵。预测步[ \begin{aligned} \hat{x}{k1|k} f(\hat{x}{k|k}) \ P_{k1|k} F_k P_{k|k} F_k^T Q_k \end{aligned} ]更新步[ \begin{aligned} K_{k1} P_{k1|k} H_{k1}^T (H_{k1} P_{k1|k} H_{k1}^T R_{k1})^{-1} \ \hat{x}{k1|k1} \hat{x}{k1|k} K_{k1}(z_{k1} - h(\hat{x}{k1|k})) \ P{k1|k1} (I - K_{k1} H_{k1}) P_{k1|k} \end{aligned} ]其中 (F_k \frac{\partial f}{\partial x}\big|{\hat{x}{k|k}})(H_{k1} \frac{\partial h}{\partial x}\big|{\hat{x}{k1|k}})。用开车导航来类比这个过程特别好理解。预测步就是“你根据上一秒的车速和方向猜测这一秒大概开到了哪里”协方差矩阵代表你对自己猜测有多大的把握——猜测越久、把握越低所以 (P) 会变大。更新步则是“GPS 报了一个位置同时你瞄了一眼路边的参照物综合两者获得一个更精确的位置”。卡尔曼增益 (K) 就是用来权衡“我预测的结果”和“传感器测量结果”各自可信度的一个系数。2.2 雅可比矩阵怎么求才不容易错EKF 最麻烦的一步就是求雅可比矩阵。很多人在推导阶段就卡住了矩阵维度一多就晕。我自己的做法分三种情况第一种情况模型简单我直接手推。比如状态向量是 (\begin{bmatrix}\delta \omega E_q\end{bmatrix}^T)观测方程是 (h(x) \begin{bmatrix}\delta \omega\end{bmatrix}^T)那么观测矩阵 (H) 就是[ H \begin{bmatrix} 1 0 0 \ 0 1 0 \end{bmatrix} ]这个很简单因为观测就是直接测量状态量。第二种情况模型复杂我直接用 Matlab 的 Symbolic Math Toolbox 求偏导然后把符号表达式转换成可调用的函数句柄。比如syms delta omega Eq real syms Pm Pe KD H ws Ts real % 状态方程 f1 omega - ws; f2 ws / (2*H) * (Pm - Pe - KD*(omega - ws)); f3 (-Eq 1) / (Td0); % 占位示例实际按模型写 F_sym jacobian([f1; f2; f3], [delta, omega, Eq]);用jacobian求出矩阵后再通过matlabFunction转成数值函数放进 EKF 循环里反复调用。这样做的好处是公式修改后不用手推但注意符号运算在循环里不能频繁调用否则速度会慢到怀疑人生。正确做法是提前生成好函数再在滤波循环外一次性编译。第三种情况实在不想求导就直接上 UKF绕开雅可比。这也是我在项目后期切换 UKF 的最大动力——模型单改一个参数EKF 的雅可比就得跟着重新推UKF 完全不用管。2.3 离散化与噪声协方差的一个隐蔽细节EKF 实现里有一个特别容易忽略的问题我们从物理模型写出的是连续时间状态方程但数字实现是离散时间的。最简单的做法是直接欧拉离散状态转移函数变为[ f_d(x_k) x_k T_s f(x_k) ]但过程噪声协方差 (Q) 也需要从连续形式转换过来。连续系统的过程噪声协方差是 (Q_c)离散化后的近似关系是 (Q_d \approx Q_c \cdot T_s)。如果直接沿用连续时间下的 (Q)滤波器的增益会偏大或偏小导致估计轨迹出现奇怪的高频抖动。在我做的算例里系统采样周期 (T_s 0.01s)连续过程噪声标准差设成 (0.001)那么离散 (Q_d) 的对角元素我一开始直接用了 (0.001^2)结果滤波曲线抖得厉害。后来把 (Q_d Q_c \cdot T_s) 换成 (10^{-8}) 量级曲线立刻平滑了很多。这种细节不写进教材但实际调试过程中遇到的概率极高。3. UKF 的实现不取导数也能传播协方差3.1 UT 变换到底做了什么UKF 的核心是 UT 变换。它的想法非常直观既然直接求非线性函数的均值和协方差很困难那我就巧妙地选几个点让这些点经过非线性函数传播后能还原出足够精确的统计特性。具体做法是对于 (n) 维状态向量生成 (2n1) 个 sigma 点[ \begin{aligned} \chi_0 \bar{x} \ \chi_i \bar{x} (\sqrt{(n\lambda)P})i, \quad i1,\dots,n \ \chi{in} \bar{x} - (\sqrt{(n\lambda)P})_i, \quad i1,\dots,n \end{aligned} ]其中 ((\sqrt{(n\lambda)P})_i) 表示矩阵 ((n\lambda)P) 的 Cholesky 分解后第 (i) 列。每个 sigma 点对应的权重为[ \begin{aligned} W_0^{(m)} \frac{\lambda}{n\lambda} \ W_0^{(c)} \frac{\lambda}{n\lambda} (1-\alpha^2\beta) \ W_i^{(m)} W_i^{(c)} \frac{1}{2(n\lambda)}, \quad i1,\dots,2n \end{aligned} ]参数 (\lambda \alpha^2(n\kappa) - n)。其中 (\alpha) 决定 sigma 点的散布程度一般取 (1e-3) 到 (1) 之间(\beta) 用于融入先验分布信息高斯分布下取 2 最优(\kappa) 是次级缩放参数通常取 0 或 (3-n)。把这些 sigma 点逐个扔进状态转移函数 (f(\cdot))得到传播后的点集然后加权计算预测均值和协方差。量测更新同理。整个过程完全没有求导运算所以对非线性函数形式没有任何平滑性要求这也是 UKF 对“硬非线性”更友好的原因。3.2 sigma 点参数怎么给在 Matlab 里做 UKF参数取值是决定成败的关键。我试过几组典型参数经验如下(\alpha 1e-3)sigma 点离均值很近适合状态量变化平缓的场景。(\alpha 1)sigma 点散布范围大状态突变时跟踪更快但稳态噪声可能偏大。(\beta 2)高斯分布下最优。(\kappa 0) 或 (3-n)如果 (n) 较大(\kappa) 取负数可能导致协方差矩阵失去半正定性这时我会改用 (\kappa0)。有一个常见坑计算 ((n\lambda)P) 的时候如果 (n\lambda \le 0)Cholesky 分解会直接失败。所以实现里必须先检查 (n\lambda 0)否则要给 (\kappa) 加一个正数偏移。这个报错在 Matlab 里通常会以Matrix must be positive definite的形式出现我一开始遇到还以为是自己协方差更新公式写错了查了半天才发现是参数问题。3.3 UKF 的 Matlab 骨架代码下面这段代码是 UKF 的核心循环骨架适合直接改成自己的模型function [x_est, P_est] ukf_predict_update(x_est, P_est, z, Q, R, f, h, param) n length(x_est); alpha param.alpha; beta param.beta; kappa param.kappa; lambda alpha^2 * (n kappa) - n; % 权重 Wm [lambda/(nlambda), 0.5/(nlambda)*ones(1, 2*n)]; Wc [lambda/(nlambda) (1-alpha^2beta), 0.5/(nlambda)*ones(1, 2*n)]; % 生成 sigma 点 A chol((nlambda) * P_est, lower); chi zeros(n, 2*n1); chi(:,1) x_est; for i 1:n chi(:,i1) x_est A(:,i); chi(:,ni1) x_est - A(:,i); end % 预测步 chi_pred zeros(n, 2*n1); for i 1:2*n1 chi_pred(:,i) f(chi(:,i)); end x_pred chi_pred * Wm; P_pred Q; for i 1:2*n1 d chi_pred(:,i) - x_pred; P_pred P_pred Wc(i) * (d * d); end % 量测更新 Z_pred zeros(size(z,1), 2*n1); for i 1:2*n1 Z_pred(:,i) h(chi_pred(:,i)); end z_pred Z_pred * Wm; Pzz R; Pxz zeros(n, size(z,1)); for i 1:2*n1 dz Z_pred(:,i) - z_pred; dx chi_pred(:,i) - x_pred; Pzz Pzz Wc(i) * (dz * dz); Pxz Pxz Wc(i) * (dx * dz); end K Pxz / Pzz; x_est x_pred K * (z - z_pred); P_est P_pred - K * Pzz * K; end这段代码里有一点需要注意P_est P_pred - K * Pzz * K是 Joseph 形式的等价写法比(I - KH)P_pred数值稳定性更好能尽量避免协方差矩阵负定。实际使用时f和h要定义成函数句柄形如f (x) x Ts * model_dyn(x, u);。如果模型状态维数不高、循环次数不多这种逐点传播的写法够用如果状态维数和采样点都很大可以考虑用矩阵批量传播来加速但代码可读性会下降。4. Matlab 实操从模型搭建到仿真出图的全流程4.1 仿真系统怎么搭我的做法是在 Matlab 里自己构造一个“数字孪生”系统先选定一个发电机模型用已知的状态初值和输入序列通过状态方程生成一条“真实状态轨迹”然后叠加高斯白噪声生成模拟量测。这样做的好处是真实状态已知可以计算估计误差并验证滤波器性能。这里我用一个简化两阶模型做演示状态量取功角 (\delta) 和转速偏差 (\omega)。状态方程为[ \begin{aligned} \delta_{k1} \delta_k T_s(\omega_k - \omega_s) \ \omega_{k1} \omega_k T_s \cdot \frac{\omega_s}{2H}(P_m - P_e(\delta_k) - K_D(\omega_k - \omega_s)) \end{aligned} ]其中电磁功率 (P_e(\delta_k) EV/x \sin(\delta_k))这是一个强非线性环节量测设置为直接测量 (\delta) 和 (\omega)并加入方差 (R) 的高斯噪声。系统参数大体如下参数数值说明(\omega_s)2pi50同步角速度H5.0 s惯性时间常数KD2.0阻尼系数Ts0.01 s采样周期/滤波步长R_diag(0.01)^2量测噪声方差Q_diag(0.001)^2过程噪声方差4.2 主程序流程整个主程序按照“生成真实轨迹 → 生成量测 → 初始化滤波器 → 滤波循环 → 误差分析”的顺序组织。核心循环部分如下% 参数初始化 n 2; % 状态维数 x_true zeros(n, N); x_true(:,1) x_init; z_meas zeros(n, N); z_meas(:,1) h(x_init) sqrt(R)*randn(n,1); % 生成真实轨迹和量测 for k 1:N-1 x_true(:,k1) f_dt(x_true(:,k), u, Ts); z_meas(:,k1) h(x_true(:,k1)) sqrt(R)*randn(n,1); end % 滤波器初始化 x_est zeros(n, N); x_est(:,1) x_init 0.1*randn(n,1); P_est eye(n); % 滤波循环 for k 1:N-1 % 预测 [x_pred, P_pred] ukf_predict(x_est(:,k), P_est, (x) f_dt(x,u,Ts), Q, param); % 更新 [x_est(:,k1), P_est] ukf_update(x_pred, P_pred, z_meas(:,k1), (x) h(x), R, param); end % 误差计算 rmse_delta sqrt(mean((x_est(1,:) - x_true(1,:)).^2));这里f_dt是离散化后的状态转移函数h是观测函数。EKF 版本只需要把ukf_predict和ukf_update换成对应的 EKF 函数即可主程序结构完全一样。4.3 Q、R 怎么调才能出好效果在动态状态估计里Q 和 R 的设置直接决定滤波效果。我的方法分三步第一步R 从量测的物理噪声水平出发。PMU 的幅值测量误差典型范围在 0.1% 到 1% 之间角度误差在 0.01 到 0.1 弧度级。据此设置 R 的对角元素这是有物理依据的不建议为了效果盲目调小。第二步Q 先给一个小值然后根据观测残差调整。如果滤波曲线滞后于真实轨迹说明 Q 太小——模型预测的“信任度”太高滤波器不跟测量走如果滤波曲线抖动剧烈说明 Q 太大——噪声被当成真实变化滤波器被测量噪声牵着鼻子走。第三步用 RMSE 做定量评估配合几个典型工况反复试。我常用的经验公式是先设 (Q q \cdot I)跑一次得到 RMSE然后 (q) 按 10 倍步长上下扫描看 RMSE 的谷值出现在哪个量级。这个方法虽然朴素但比凭感觉调参可靠得多。4.4 结果图怎么画才直观滤波做完后画图是必不可少的。我最常用的画图组合是第一张图把真实状态、带噪声量测、EKF估计值、UKF估计值四条曲线画在一起直观对比跟踪效果第二张图画估计误差并画出正负 3 倍标准差的置信区间第三张图画协方差矩阵对角线元素的变化趋势观察滤波器是否收敛。画图代码大概长这样figure; subplot(2,1,1); plot(t, x_true(1,:), k-, LineWidth, 1.5); hold on; plot(t, z_meas(1,:), b., MarkerSize, 3); plot(t, x_est_ekf(1,:), r--, LineWidth, 1.2); plot(t, x_est_ukf(1,:), g-., LineWidth, 1.2); legend(真实值, 量测值, EKF, UKF); xlabel(时间 (s)); ylabel(功角 (rad)); grid on;为什么要把量测值也画上去因为这样可以直观看出滤波到底“平滑”了多少。如果滤波结果和量测噪声几乎重合那说明滤波增益过大相当于没滤波如果滤波结果和量测偏差太大说明增益过小跟踪跟不上。5. 常见问题与排错实录5.1 滤波器发散怎么排查我遇到过的最典型的发散现象是估计值在第一、第二步就被量测噪声带飞之后再也回不到真实轨迹附近协方差矩阵也快速膨胀甚至变成 NaN。排查顺序我一般按下面这个清单症状可能原因对策估计值完全偏离真实值初值设置离真实值太远用静态估计结果或潮流解初始化前几步就出现 NaN协方差矩阵非正定用 Joseph 形式更新检查 sigma 点参数稳态时误差仍然很大Q 过大或 R 过小减小 Q、增大 R用残差统计调整响应太慢、滞后明显Q 过小或 R 过大增大 Q、减小 R轨迹跟随量测噪声高频抖动Q 过大或 R 过小检查 Q 是否忘了乘采样周期5.2 状态量纲不一致导致的数值问题功角的量纲是弧度转速偏差可能取极标幺值也可能取弧度每秒两者数值量级可能差了 10^3 以上。如果不做归一化协方差矩阵里不同对角元的数量级差太多Cholesky 分解很容易数值不稳定。我后来把所有状态量都做了归一化处理功角用弧度转速偏差用标幺值 ( \Delta\omega / \omega_s )电磁功率用标幺值。全部统一到 0.1 到 10 这个量级范围滤波稳定性和精度都好了很多。这也提醒我卡尔曼滤波虽然对量纲不敏感但数值实现里量级差距过大会引发各种隐蔽问题。5.3 初值到底该怎么给动态状态估计的初值不能随便给。EKF 和 UKF 本质上都是局部递推算法初值离真实状态太远第一轮更新就可能让协方差矩阵失去意义。工程上常用的做法是先用一个常规静态状态估计器得到初始断面或者直接用潮流计算结果初始化。如果实在没有可以先把滤波器放在“稳态工况”下预热跑一段时间用滤波结果反推初值。我在单机模型算例里试过初值偏差在 5% 以内时滤波器都能快速收敛偏差超过 20% 时UKF 还能勉强拉回来EKF 就已经发散得没法看了。5.4 计算效率能不能再快一点Matlab 里做 UKF 容易陷入双重循环导致仿真时长暴增。我的优化策略有两个第一sigma 点传播尽量用矩阵运算代替循环。把 2n1 个 sigma 点拼成一个 (n, 2n1) 的矩阵然后用“一次函数求值”的方式处理整个矩阵但前提是f和h能向量化。第二避免在滤波循环内重复分配大数组。提前把chi,chi_pred,Z_pred等变量初始化好循环内只做赋值。第三个策略是能离线算的都不在线算。比如状态转移函数的某些固定结构、雅可比矩阵的符号表达式都提前算好。之前我犯过在 EKF 循环里每次都调matlabFunction的错那个速度之慢跑一次仿真够我喝三杯咖啡。6. 一次完整的调试经历最后分享一个具体的调试案例。当时我在单机无穷大系统上做 EKF参数全部按论文里的标准值设定结果滤波曲线在第 1 秒就开始剧烈振荡然后直接发散。第一反应是雅可比矩阵求错了。我重新用符号工具求了一次偏导逐项对比没发现错。后来我把过程噪声 (Q) 的对角元素从 (10^{-4}) 一步步往下调调到 (10^{-8}) 之后滤波曲线才稳定下来。原因在于我的状态转移方程里转速偏差的变化率数值本身就很小过程噪声给大了状态预测的协方差膨胀很快导致滤波器过度信任量测而量测噪声又存在于是估计值被量测噪声完全主导。这个教训让我彻底记住了Q 和 R 的绝对数值不重要重要的是它们和状态方程中各项系数的相对关系。这也是为什么很多论文里直接给一个 Q、R 矩阵你照搬过来却跑不出同样效果的原因——你的模型参数、单位、离散化方式都和作者不一样。后来我在同一套模型上把 UKF 也跑通对比 RMSE 发现 UKF 比 EKF 的估计误差平均小 30% 左右尤其在功角突变阶段。虽然 UKF 的计算时间大约是 EKF 的 1.8 倍但在这个算例规模下完全可接受。综合精度和调参成本我现在更推荐在电力系统动态状态估计中优先尝试 UKF。如果后续你有兴趣继续深入可以考虑把量测扩展为“伪量测PMU 混合量测”或者在 EKF/UKF 基础上加入鲁棒性设计应对量测异常和参数不确定性。这些都是很有价值的扩展方向但前提是把基础滤波流程彻底跑通。