1. 动态状态估计为什么传统的静态估计在电网调度中开始“跟不上”了1.1 从断面快照到动态递推状态估计的两次思路跃迁做电力系统状态估计的同行都知道传统调度自动化系统里跑的基本都是静态状态估计最常见的就是加权最小二乘WLS。它的思路是拿当前一个时间断面的所有量测构建一组方程迭代求解出系统各节点的电压幅值和相角。好处是理论成熟、程序健壮实际运行了这么多年稳定可靠。但它的本质是一个“静态拍照”过程——拿到一张断面图算出此刻系统最可能的状态然后这张图就作废了下一个时刻又重新拍一张。这种思路在稳态工况下没太大问题但一旦系统进入动态过程问题就来了。比如新能源大规模接入之后出力波动频繁母线电压和线路功率的短时变化比过去剧烈得多又比如发生扰动或故障后的暂态过程中功角、角速度这些量以秒甚至毫秒级的尺度在演化这时候静态估计的“断面快照”模式就明显不够用。它没有利用系统自身的演化规律每一拍都从零开始算既浪费了历史信息也对量测误差的抑制不够充分。动态状态估计Dynamic State Estimation, DSE则是换了一个思路把电力系统当成一个动态系统来对待用状态方程描述它的演化规律再用量测方程把状态和观测量联系起来。每一时刻的估计结果一方面来自上一时刻状态的预测递推另一方面来自当前时刻的量测修正两者按各自的置信度加权融合。这就有点像开车时导航的做法车本身的惯性模型告诉你“我大概开到哪了”GPS量测告诉你“我实际观测到在哪”两者结合出来的位置远比单独依赖某一个要准。正是这个“预测修正”的递推结构让动态状态估计具备了静态估计没有的两个能力一是能输出状态量的预测值这在预警和稳定控制中非常有用二是能在量测质量差、甚至短暂缺失时依靠动态模型维持一段时间有效估计。EKF和UKF就是在这条思路上最常用的两种非线性滤波实现手段。1.2 动态状态估计的数学框架状态方程、量测方程与噪声建模要把动态状态估计讲清楚首先得把它的数学表达式立起来。电力系统中我们关心的是发电机的功角δ、角速度ω这类动态状态变量。以最经典、也是最常用来做滤波算法验证的发电机二阶模型摇摆方程为例状态方程δ(k1) δ(k) Δω(k)·Δt w₁(k)Δω(k1) Δω(k) (Δt/M)·[P_m - P_e(δ(k)) - D·Δω(k)] w₂(k)其中M是发电机惯性时间常数D是阻尼系数P_m是机械功率P_e是电磁功率。电磁功率和功角的关系是非线性的P_e(δ) (E′_q·V_s / X′_d)·sin(δ)这里的E′_q是发电机暂态电动势V_s是无穷大母线电压X′_d是暂态电抗。是不是已经闻到非线性的味道了对功角和功率之间差着一个正弦函数这个非线性是躲不掉的。量测方程则根据实际装设的PMU量测来定常见的有z₁(k) δ(k) v₁(k)z₂(k) Δω(k) v₂(k)z₃(k) P_e(δ(k)) v₃(k)w和v分别是过程噪声和量测噪声工程上通常假设为相互独立的高斯白噪声协方差矩阵分别为Q和R。动态状态估计的任务就是基于这些带噪声的量测序列实时估计出真正的状态量δ和ω。刚才说的是单机无穷大OMIB系统的模型。扩展到多机系统时状态向量会变成所有发电机的功角和角速度拼在一起的大向量量测向量也相应变成多个PMU断面的拼接但整体的“状态方程量测方程”框架完全一致。这也是为什么读懂单机系统的滤波流程之后往多机系统扩展是水到渠成的事。1.3 非线性是动态估计绕不开的坎讲到这里有一个问题自然而然浮出来既然状态方程中已经有δ和ω的线性递推关系为什么不能直接套用线性卡尔曼滤波答案就出在量测方程上。P_e(δ)和δ之间的正弦关系、以及多机系统中更复杂的功率方程都是不可忽略的非线性函数。如果把非线性直接忽略用线性卡尔曼滤波硬算状态估计结果大概率要偏离严重时滤波器直接发散。所以问题的核心就变成在非线性条件下如何完成卡尔曼滤波框架中的预测步和更新步EKF和UKF给出的答案代表了两种典型的技术路线EKF的思路是“把非线性问题线性化”——对状态方程和量测方程在当前估计点附近做一阶泰勒展开然后套用标准卡尔曼滤波公式UKF的思路是“对概率分布做采样逼近”——用一组精心挑选的sigma点去逼近状态分布经过非线性变换后的真实均值和协方差完全绕开了雅可比矩阵的计算。这两种思路各自的适用场景、实现难度、数值稳定性都有差异下面两章我展开讲。2. EKF与UKF处理非线性滤波问题的两条路线2.1 先从标准卡尔曼滤波说起为什么它不能直接用在电力系统里要说清EKF和UKF的差异必须先回到标准卡尔曼滤波本身。标准KF针对的是线性系统x(k1) F·x(k) w(k)z(k) H·x(k) v(k)其中F和H都是常矩阵。它的核心思想特别朴素用状态方程做一步预测得到状态预测量x_pred和预测协方差P_pred然后计算卡尔曼增益K用量测残差对预测值做修正。整个过程只需要矩阵乘法和求逆运算计算量小而且在线性高斯条件下是理论最优估计器。问题是电力系统的量测方程根本不是线性的。H不是一个常数矩阵而是状态量的函数H(x)。这时候如果依旧套用线性KF的公式预测协方差和量测协方差的传递就会出现系统性偏差。事实上我见过不少初学者把P_e(δ)和δ之间的关系做线性近似直接用常数增益算结果滤波轨迹和真值之间总是隔着一道“鸿沟”怎么调Q和R都补不回来。本质上就是因为丢弃了非线性项带来的信息。正因为如此工程中必须引入专门处理非线性滤波的算法。EKF和UKF的出发点相同都是把标准KF的思想推广到非线性场景但推广的手段南辕北辙。2.2 EKF一阶泰勒展开思路简单但代价是精度和稳定性的妥协扩展卡尔曼滤波的做法是在每个时刻把非线性函数在当前状态估计点附近做一阶泰勒展开。具体来说对于状态方程x(k1) f(x(k)) w(k)会计算状态转移雅可比矩阵F(k) ∂f/∂x | x x̂(k)对于量测方程z(k) h(x(k)) v(k)计算量测雅可比矩阵H(k) ∂h/∂x | x x_pred(k)有了这两个雅可比矩阵就可以套用和标准KF几乎一样的形式预测步x_pred(k1) f(x̂(k))P_pred(k1) F(k)·P̂(k)·F(k)ᵀ Q更新步K(k1) P_pred(k1)·H(k1)ᵀ·[H(k1)·P_pred(k1)·H(k1)ᵀ R]⁻¹x̂(k1) x_pred(k1) K(k1)·[z(k1) - h(x_pred(k1))]P̂(k1) [I - K(k1)·H(k1)]·P_pred(k1)单看公式EKF确实很好理解代码也不难写。但它的短板同样明显第一一阶泰勒展开本身就是一种近似只保留了函数的线性主部丢掉了二阶及以上的信息。当系统非线性程度较强时比如功角在暂态过程中大幅摆动、P_e(δ)的导数变化剧烈这种一阶近似带来的误差会累积严重时导致滤波发散。第二雅可比矩阵的计算是个麻烦事。对单机二阶模型手动推导F和H还算轻松一旦扩展到多机系统、状态量几十维、量测方程包含各种功率表达式时手动推导雅可比矩阵就是一场灾难。很多人会改用数值差分来替代解析求导但数值差分本身又引入步长选择的困扰步长太大误差大步长太小数值不稳定。第三EKF对初值和噪声统计的敏感性比较高。初值偏差太大时线性化点离真实点太远一阶近似的有效性大打折扣滤波器容易在启动阶段就“带偏”后面再怎么调都很难拉回来。那是不是说EKF就没用了不是。在弱非线性、且对计算速度要求高的场景下EKF依然是性价比很高的选择。毕竟它每次递推只需要做两次矩阵乘法、一次矩阵求逆计算开销远小于后面要讲的UKF。如果系统模型本身特性比较“温和”状态演化始终在平衡点附近小范围波动那么EKF的一阶近似误差完全在可接受范围内。2.3 UKF不碰雅可比用sigma点直接逼近概率分布无迹卡尔曼滤波换了一个完全不同的角度切入。它不去线性化非线性函数而是去逼近状态分布经过非线性变换后的统计特性。这个思路的实现载体就是无迹变换Unscented Transform, UT。UT的核心思想是对于n维状态向量x取其均值点x̄再加上2n个sigma点这些点按照一定规则分布在均值周围并且被赋予不同的权重。然后把每个sigma点分别通过非线性函数传播从传播后的sigma点集合中加权计算输出状态的均值和协方差。这样做的好处是不需要计算任何雅可比矩阵只需要把非线性函数当“黑箱”调用若干次就能以较高的精度逼近真实分布。sigma点选取的经典公式如下χ₀ x̄χᵢ x̄ (√((nλ)·P))ᵢi 1, ..., nχᵢ₊ₙ x̄ - (√((nλ)·P))ᵢi 1, ..., n对应的权重为W₀ᵐ λ/(nλ)W₀ᶜ λ/(nλ) (1 - α² β)Wᵢᵐ Wᵢᶜ 1/[2(nλ)]其中λ α²(nκ) - nα是一个控制sigma点散布范围的参数通常取一个较小的正数比如1e-3κ是次级缩放参数通常取0β用来融入分布先验信息高斯分布时取2。UT完成之后后续的预测步和更新步就是在sigma点层面上展开。预测步中把sigma点通过状态方程传播加权算出状态预测均值和预测协方差更新步中把预测的sigma点通过量测方程传播算出量测预测均值、量测协方差以及状态和量测的互协方差然后计算卡尔曼增益并更新状态。整个过程全部避开了雅可比矩阵算法对非线性的适应能力大幅度提升。理论上UT变换可以精确捕捉到非线性函数传播后的二阶矩信息在精度上相当于二阶泰勒展开而EKF只保留了一阶项这正是UKF在强非线性场景下误差更小的根本原因。不过UKF也不是没有代价。它的计算量明显大于EKF——每步递推需要2n1次非线性函数传播状态维数越高sigma点越多计算开销越大。另外sigma点选取涉及多个超参数虽然经验默认值能覆盖大多数情况但真到具体问题时也需要根据状态尺度做适当调整。2.4 EKF vs UKF一张表看透两者的实质差异我把两种算法放在一起做对比。这张表是我在实际项目里反复对照后梳理出来的比单纯看公式直观得多对比维度EKFUKF核心思想将非线性函数在估计点做一阶泰勒展开用sigma点逼近非线性变换后的概率分布是否需要雅可比矩阵需要手动推导或数值差分不需要对非线性的适应能力弱只保留一阶项强相当于二阶精度计算复杂度低适合高维大系统在线计算偏高需2n1次非线性传播对初值偏差的敏感性敏感初值差容易发散更稳健sigma点覆盖范围大实现难度公式直观但雅可比推导繁琐无需求导但参数和采样规则需要理解适用场景弱非线性、计算资源紧张、在线实时强非线性、对精度要求高、离线或在线均可用注意这张表不是绝对的。实际项目中我见过用EKF跑得稳如泰山的多机系统也见过UKF因为超参数没调好而对量测噪声特别敏感的反例。核心结论是算法的选择必须结合具体系统的非线性程度和算力约束不能简单说谁一定更好。3. Matlab仿真框架搭建与核心代码拆解3.1 整体架构设计模块化思路让滤波对比更清晰在用Matlab实现EKF和UKF之前我建议先想清楚代码的整体组织方式。我看到不少初学者把滤波算法、系统建模、结果绘图全塞在一个脚本里几百行代码铺下来出问题的时候根本不知道是模型写错了还是滤波逻辑写错了。正确的做法是模块化拆分。我给一个建议的工程目录结构dse_project/ ├── main_RunDSE.m # 主入口设置参数、调用滤波算法、绘制结果 ├── system_model.m # 状态方程与量测方程定义 ├── generate_true_states.m # 生成真值轨迹模拟系统的真实动态 ├── generate_measurements.m # 基于真值叠加噪声生成量测序列 ├── ekf_estimation.m # EKF滤波主函数 ├── ukf_estimation.m # UKF滤波主函数 ├── compute_rmse.m # 计算均方根误差 └── plot_results.m # 绘图脚本main脚本负责把其他函数串起来这样换系统参数、换滤波算法、换初始条件都很方便。比如想对比EKF和UKF只需要在main里分别调用两个函数把结果存下来即可。想换系统模型只改system_model.m一个文件其他函数无需动。这种结构的另一个好处是方便查错。滤波代码里的一个符号错误可能只在特定参数组合下才暴露如果全堆在一个脚本里排查的难度会几何级上升。3.2 系统建模单机无穷大系统的状态方程与量测方程实现为了把代码讲透我用单机无穷大OMIB系统作为演示对象。它虽然简化但包含了动态状态估计的全部核心要素而且状态量只有二维结果可视化非常清楚。先把模型参数列出来% 系统参数 M 10; % 惯性时间常数秒 D 1.5; % 阻尼系数 Xd_prime 0.3; % 暂态电抗标幺值 Eq 1.2; % 发电机暂态电动势标幺值 Vs 1.0; % 无穷大母线电压标幺值 Pm 1.0; % 机械功率标幺值 dt 0.01; % 采样间隔秒 T 20; % 仿真时长秒状态向量x [δ, Δω]ᵀ其中δ是功角弧度Δω是角速度偏差标幺值ω - ω_s。状态方程和量测方程在Matlab里定义如下function [x_next, y] system_model(x, dt, params) % 状态方程 f(x) delta x(1); domega x(2); Pe (params.Eq * params.Vs / params.Xd_prime) * sin(delta); x_next zeros(2,1); x_next(1) delta domega * dt; x_next(2) domega (dt / params.M) * (params.Pm - Pe - params.D * domega); % 量测方程 h(x)这里选了3个量测功角、角速度和电磁功率 y zeros(3,1); y(1) delta; y(2) domega; y(3) Pe; end这里有个细节值得注意量测方程中同时包含直接测量量功角、角速度和间接测量量电磁功率。加入功率量测的目的是让滤波器有更多信息来校正状态。PMU可以直接测功角但由于相量测量单元本身的算法和通信环节实际测到的功角依然有噪声角速度和功率量测提供了不同的误差特性多量测融合能显著提升估计精度。3.3 真值数据生成与量测模拟仿真研究的必备一步在仿真验证滤波器性能时我们不能直接拿真实电网数据开跑而是先构造一个“已知答案”的场景生成一条确定性的状态演化轨迹作为真值再叠加噪声模拟量测。这样后续才能定量地比较滤波结果和真值的偏差。真值轨迹的生成代码很简单N round(T / dt); x_true zeros(2, N); y_true zeros(3, N); x_true(:,1) [0.3; 0.0]; % 初始功角0.3rad角速度偏差0 for k 1:N-1 [x_next, y] system_model(x_true(:,k), dt, params); x_true(:,k1) x_next; y_true(:,k) y; end y_true(:,N) system_model(x_true(:,N), dt, params);这里设定了一个初值功角0.3rad角速度偏差0。系统会从这组初值开始演化因为机械功率和电磁功率不平衡功角会出现一段动态过渡过程。这段动态过程恰恰是检验滤波器跟踪能力的“试金石”。量测数据则在真值上加高斯白噪声R diag([1e-4, 1e-4, 1e-3]); % 量测噪声协方差 z_meas y_true sqrt(R) * randn(3, N);注意这里R的含义功角量测噪声的标准差为0.01弧度约0.57度角速度为0.01标幺值功率量测噪声标准差约为0.032标幺值。这个量级比较贴近PMU的实际测量误差范围既不会太理想化也不会让滤波变得过于困难。3.4 EKF核心实现预测-更新循环与雅可比矩阵的处理方式EKF的实现在原理上很直接但实际写代码时有两个容易出错的点一是雅可比矩阵的计算二是预测步和更新步中雅可比矩阵所使用的状态点不同预测步用的是上一时刻的估计值更新步用的是当前时刻的预测值。我用数值差分的方式计算雅可比矩阵这样耦合到多机系统时不需要手动推公式function F compute_F(x, dt, params) % 状态方程对x的雅可比矩阵数值差分法 n length(x); eps 1e-8; F zeros(n, n); [x_base, ~] system_model(x, dt, params); for j 1:n x_pert x; x_pert(j) x_pert(j) eps; [x_pert_out, ~] system_model(x_pert, dt, params); F(:, j) (x_pert_out - x_base) / eps; end end function H compute_H(x, dt, params) % 量测方程对x的雅可比矩阵数值差分法 m 3; n length(x); eps 1e-8; H zeros(m, n); [~, y_base] system_model(x, dt, params); for j 1:n x_pert x; x_pert(j) x_pert(j) eps; [~, y_pert] system_model(x_pert, dt, params); H(:, j) (y_pert - y_base) / eps; end endEKF主循环function [x_est, P_est] ekf_estimation(z_meas, x0, P0, Q, R, dt, params) N size(z_meas, 2); n length(x0); x_est zeros(n, N); P_est zeros(n, n, N); x x0; P P0; for k 1:N % 预测步 F compute_F(x, dt, params); [x_pred, ~] system_model(x, dt, params); P_pred F * P * F Q; % 更新步 H compute_H(x_pred, dt, params); [~, y_pred] system_model(x_pred, dt, params); S H * P_pred * H R; K P_pred * H / S; x x_pred K * (z_meas(:,k) - y_pred); P (eye(n) - K * H) * P_pred; x_est(:,k) x; P_est(:,:,k) P; end end这段代码需要注意的细节是计算H时用的是x_pred预测值而不是上一时刻的x。我刚学EKF的时候在这里栽过跟头把H和F用同一个点的雅可比矩阵算结果滤波结果始终差一口气。原理上预测步的线性化点应该是上一时刻的后验估计而更新步的线性化点应该是当前时刻的先验预测这个顺序不能反。3.5 UKF核心实现sigma点生成、权重计算与状态更新UKF的代码稍微复杂一些但只要理解了sigma点和权重的概念写起来并不困难。核心函数如下function [x_est, P_est] ukf_estimation(z_meas, x0, P0, Q, R, dt, params) N size(z_meas, 2); n length(x0); alpha 1e-3; kappa 0; beta 2; lambda alpha^2 * (n kappa) - n; % 权重 Wm zeros(2*n1, 1); Wc zeros(2*n1, 1); Wm(1) lambda / (n lambda); Wc(1) lambda / (n lambda) (1 - alpha^2 beta); for i 1:n Wm(i1) 1 / (2*(n lambda)); Wm(in1) 1 / (2*(n lambda)); Wc(i1) 1 / (2*(n lambda)); Wc(in1) 1 / (2*(n lambda)); end x_est zeros(n, N); P_est zeros(n, n, N); x x0; P P0; for k 1:N % 生成sigma点 sqrt_P sqrtm((n lambda) * P); chi zeros(n, 2*n1); chi(:,1) x; for i 1:n chi(:,i1) x sqrt_P(:,i); chi(:,in1) x - sqrt_P(:,i); end % 状态传播 chi_pred zeros(n, 2*n1); gamma_pred zeros(3, 2*n1); for i 1:2*n1 [chi_pred(:,i), gamma_pred(:,i)] system_model(chi(:,i), dt, params); end % 计算预测均值和协方差 x_pred zeros(n, 1); for i 1:2*n1 x_pred x_pred Wm(i) * chi_pred(:,i); end P_pred Q; for i 1:2*n1 diff chi_pred(:,i) - x_pred; P_pred P_pred Wc(i) * (diff * diff); end % 量测传播 y_pred zeros(3, 1); for i 1:2*n1 y_pred y_pred Wm(i) * gamma_pred(:,i); end P_yy R; P_xy zeros(n, 3); for i 1:2*n1 dy gamma_pred(:,i) - y_pred; dx chi_pred(:,i) - x_pred; P_yy P_yy Wc(i) * (dy * dy); P_xy P_xy Wc(i) * (dx * dy); end % 更新 K P_xy / P_yy; x x_pred K * (z_meas(:,k) - y_pred); P P_pred - K * P_yy * K; x_est(:,k) x; P_est(:,:,k) P; end end这段代码里我选择用sqrtm函数来求解矩阵平方根。对于协方差矩阵P理论上它是正定对称的sqrtm的结果也是对称的。不过在数值计算中如果P出现非正定这在EKF里比较容易见到sqrtm会报错或给出复数结果。实际使用时可以先用chol函数做Cholesky分解速度更快但对矩阵正定性要求更高。我更推荐在每一步更新后对P做一次对称化处理即P (P P) / 2避免累积的舍入误差破坏对称性。UKF相比EKF还有一个细节优势这里的P_pred和P_yy计算都直接用了Wc权重对差分外积做加权求和整个过程不需要任何求导。当系统模型改动时比如从二阶模型换成三阶模型只需要修改system_model函数UKF部分一个字符都不用动这对做研究和对比实验来说特别友好。3.6 评估指标RMSE才是衡量滤波效果的唯一标尺滤波器跑完了如何评价好坏总不能拿“看起来差不多”这种主观判断来下结论。业内最常用的指标就是均方根误差RMSERMSE sqrt((1/N) · Σ(k1 to N) (x̂(k) - x_true(k))²)对于多状态量可以分别计算每个状态的RMSE也可以做一个综合的加权平均。我的习惯是分开算因为功角和角速度的量纲不同数值尺度差异很大合在一起算会掩盖单个状态的性能差异。function rmse compute_rmse(x_est, x_true) diff x_est - x_true; rmse sqrt(mean(diff.^2, 2)); end这个函数简单得不能再简单但它是横向对比EKF和UKF性能的核心依据。在第四章的仿真结果中我会重点关注两个问题一是RMSE的数值二是滤波轨迹是否在暂态过程中出现明显偏差。4. 仿真结果对比EKF与UKF在同一测试场景下的实战表现4.1 测试场景设计初值偏差和暂态过程是检验滤波器最好的“试金石”为了让两种算法的差异真正显现出来我特意把测试场景设计得“刁钻”一些。首先滤波器的初始状态不设成真值而是设一个明显偏差的初值x0 [0.5; 0.05]; % 滤波器初值功角0.5rad角速度偏差0.05而真值是从功角0.3rad、角速度偏差0开始演化的。也就是说滤波器在启动瞬间就带着一个不小的偏差无论是EKF还是UKF都必须依靠量测信息一点点把状态“拉”回正确轨道。这个过程正好可以观察两种算法的收敛速度。过程噪声和量测噪声的设置也很有讲究Q diag([1e-6, 1e-4]); % 过程噪声协方差 R diag([1e-4, 1e-4, 1e-3]); % 量测噪声协方差Q的量级必须和模型的不确定性匹配。这里把功角的Q设得比角速度小一个量级是因为功角方程本质上是对角速度的积分积分本身会平滑高频噪声而角速度方程直接受到功率不平衡的影响不确定性更大。4.2 状态估计轨迹从图形上看两种算法的跟踪能力差异我跑完仿真之后的典型输出是三条曲线真值轨迹、EKF估计轨迹、UKF估计轨迹。在功角图中三条曲线在前0.5秒内会有明显分离因为滤波器正在从初始偏差中恢复。但之后的走势会有分化EKF的功角估计曲线整体能跟上真值的趋势但在暂态过程的快速阶段功角变化速率最大的区间会出现明显的滞后。这个滞后不是滤波器参数调得不好造成的而是EKF一阶线性化和真实非线性函数之间的偏差在起作用。具体来说当功角变化剧烈时P_e(δ)的导数和二阶导数都很大一阶泰勒展开丢掉的二阶项不可忽略导致预测步的误差偏大。UKF的曲线则更贴近真值尤其在暂态中期和后期两者的偏差肉眼几乎看不出来。原因不难理解UKF通过sigma点传播不需要在某个点做线性化近似它对强非线性函数的均值和协方差传播更准确。特别是在系统从暂态向稳态过渡的阶段UKF表现出更平稳的收敛过程没有EKF那种“過冲”现象。角速度的估计差异同样明显。角速度本身数值较小绝对值在零点零几的范围内波动EKF在快速变化区间的估计抖动明显更大UKF则平稳得多。这种差异在RMSE指标上会体现得非常直接。4.3 量化对比RMSE数值揭示的真相我这里给出一个具有代表性的RMSE对比结果具体数值会因噪声随机种子和参数设置而略有浮动但相对趋势是稳定的状态量EKF RMSEUKF RMSEUKF相对误差改善功角δ弧度0.01180.0062约47%角速度Δω标幺值0.00410.0023约44%从这个结果能看出来UKF在这个测试场景下的RMSE大致是EKF的一半左右。这个差距的绝对值很小但放到电力系统动态状态估计的语境里意义完全不同。功角RMSE从0.012弧度降到0.006弧度意味着对发电机状态不确定性的掌握精度提升了近一倍在功角稳定裕度判断中这会直接影响预警阈值设定。不过我要强调一个容易被忽略的点这个RMSE对比是在“量测信息不足”的时候最悬殊。如果我把量测噪声R调得更小意味着量测更可靠两种算法的差距会缩小反之如果R调大量测更不可靠模型预测的作用增强UKF在模型传播上的精度优势会更加放大。理解这个规律比单纯记住“UKF比EKF好”重要得多。4.4 发散风险与数值稳定性那些“看起来对”的代码是怎么翻车的在仿真过程中我遇到过几个典型问题写出来给各位避坑。第一个是EKF的协方差矩阵P出现非正定。根因是浮点舍入误差在多次递推中累积破坏了P的对称正定性。一旦P非正定后续计算里会蹦出负方差、甚至NaN滤波器直接报废。解决方法是每步更新后强制对称化并且可以用特征值分解把P的很小负特征值“截断”到零。UKF因为用sigma点传播协方差对P的正定性要求同样严格但我的实测经验是UKF出现非正定的概率比EKF低一些因为它没有显式的雅可比矩阵乘法。第二个是数值差分计算雅可比矩阵时的步长选择。我一开始用eps 1e-6结果在状态量数值本身就很小的情况下比如Δω只有0.01量级微扰步长和状态值处于同一量级差分的截断误差和舍入误差都快到极限了算出来的雅可比矩阵噪声很大。后来把步长调到1e-8效果明显改善。如果状态量尺度差异很大还可以逐状态设置差分步长。不过最稳妥的做法还是解析推导雅可比矩阵只是多机系统下推导工作量确实不小。第三个是UKF的参数α。我最早直接用默认的1e-3在单机系统上没问题但后来把算法搬到另一个状态量尺度差异更大的系统时发现sigma点生成的矩阵平方根在数值上出现退化。适当增大α到1e-2问题就消失了。这里的经验是α控制sigma点离均值的距离状态量尺度小的时候α太大会让sigma点覆盖范围过大状态量尺度大的时候α太小又会让sigma点挤在一起需要根据实际情况微调。5. 工程落地的实用建议与调参经验5.1 过程噪声矩阵Q的调参逻辑不是随便给个对角阵就行我在做动态状态估计项目时发现最容易被忽视的参数就是Q矩阵。很多初学者直接把Q设成一个很小的对角阵理由是“我认为模型很准”。但这样做往往会导致滤波器对量测的信任度过低量测信息没有充分发挥作用反过来Q设得太大滤波器又会过度相信量测导致估计结果在噪声中剧烈抖动。一个实用的调参方法是先根据系统模型的不确定性源来估计Q的下界。比如状态方程里忽略了多少高阶动态、模型参数M、D的辨识误差有多大、离散化带来的截断误差有多少把这些不确定性的方差加起来作为Q的初始估计。然后在这个基础上用一段真实或仿真数据反复试验观察滤波残差innovation的自相关特性如果残差存在明显的时间相关性说明Q设小了如果残差的方差远大于R的预测值说明Q设大了。对于单机系统我用的Q diag([1e-6, 1e-4])算是比较可靠的经验值。多机系统下各台发电机之间可能存在动态耦合更严谨的做法是建立完整的关联Q矩阵不过实际项目中为了简化通常只保留对角项。这种简化在大多数情况下是合理的只要量测冗余度足够滤波结果对Q的非对角项的敏感性不高。5.2 量测异常与时延动态估计在真实系统中面临的两大“暗礁”仿真环境里量测是理想化的——每个时刻都有数据、没有丢包、没有时延。但真实的PMU数据流远没有这么干净。我遇到过两种情况让滤波器几乎崩溃。一是量测野值。电网运行中PMU通道偶尔会跳出明显异常的数据点比如幅值突变到其他量级。如果滤波器直接把这个野值当作有效量测更新状态估计值会被瞬间“拽”到错误位置之后需要好几个周期才能恢复。工程上常用的处理办法是增加一个“创新限幅”逻辑计算量测残差zk - h(x_pred)如果残差超过某个阈值比如3倍量测噪声标准差就认为该量测是野值在本步更新中将其丢弃或降权。这个逻辑实现起来不难但对滤波器的稳定性帮助极大。二是量测时延。PMU从采样到数据上送中间要经过合并单元、相量数据集中器等多级环节时延在几十到几百毫秒不等。而动态状态估计的递推周期通常只有10到20毫秒如果不对时延做处理滤波算法拿到的“当前量测”其实是几十个周期前的数据状态已经往前跑了一截了。这种情况下一个常见的妥协做法是采用“延迟量测处理”在k时刻收到的是k-d时刻的量测时先用状态方程把滤波器的状态递推到当前时刻再用滞后量测做更新最后再递推回当前时刻。另一种更工程化的思路是在滤波器设计时预留一个量测缓冲并让滤波周期和PMU的数据周期保持一致从根本上规避时延问题。5.3 从单机到多机代码扩展路径和需要小心的坑单机无穷大系统的滤波跑通之后向多机系统扩展几乎是必然的路径。这个扩展在算法层面并不需要大改EKF和UKF的框架完全不变需要改的是系统模型状态向量从[δ₁, Δω₁]变成[δ₁, Δω₁, δ₂, Δω₂, ..., δₙ, Δωₙ]量测向量相应扩充到所有PMU的断面量测。量测方程会变成更复杂的网络功率方程这时EKF的雅可比矩阵就必须依赖数值差分或者专门的解析推导工具而UKF的优势更加凸显——只要系统模型函数能正确算出量测输出算法就能工作。一个经常被忽视的问题是发电机模型的阶数。前面用的是二阶经典模型只包含功角和角速度。如果要做更精确的动态估计发电机的励磁系统动态不能忽略这时需要引入三阶模型状态方程中再加入暂态电动势E′q的动态量测方程也要对应增加励磁电压相关的量测。状态维数从2维变成3维甚至更高UKF的sigma点数量随之增加计算量会明显上升。但如果你要估计的是功角稳定性相关的状态二阶模型往往已经够用需要分析励磁系统行为时再升级到三阶不必一开始就追求高保真模型。5.4 一个被低估的细节PMU量测的参考相角问题最后说一个我踩过坑、但很多文献里都不提的细节PMU测得的功角是相对时间同步信号的绝对相角而动态状态估计中的功角通常是相对系统参考机或者相对自身初始值的相对量。这两者之间如果不做换算滤波结果会产生一个常数偏置。更麻烦的是参考相角本身会随着系统频率的变化而漂移尤其是系统偏离额定频率时这个漂移会污染整个估计。实际处理办法是先对PMU功角量测做基准变换把绝对相角转换为相对系统参考轴的相对值。通常可以用系统惯量中心COI作为参考或者用一台参考机组的功角作为基准。这个预处理步骤看似不起眼但没有它滤波器的估计误差里会始终带着一个缓慢变化的系统误差RMSE指标再漂亮也掩盖不了这个隐患。5.5 我的总体体会没有“最好”的滤波器只有“适合”的滤波器经过这几次仿真对比我对EKF和UKF的选择有了一个更务实的判断。如果你要做的是数十台以上发电机的在线动态状态估计计算资源有限系统运行点相对平稳弱非线性假设大体成立EKF在今天依然是值得优先考虑的方案。它的公式结构清晰计算量小工程实现成熟出现问题时排查也容易。如果你更看重估计精度且系统可能经历较大的动态过程——比如暂态稳定分析、扰动后的轨迹估计——UKF的多花的那点计算成本完全值得。一次暂态过程的仿真时长也就十几秒哪怕UKF单步耗时是EKF的三倍总时长依旧在毫秒级这个代价换来的精度提升远远抵得过那点计算开销。在Matlab里把这两个算法放在同一个框架下对比测试是我觉得效率最高的做法。先跑通一个简单模型把滤波框架和评估流程固定下来再去扩展模型复杂度、替换数据源、调试噪声参数整个流程的迭代速度会快很多。这也是我写这篇内容的原因——希望更多人能快速迈过“滤波代码跑通”这道门槛把精力放在真正有意思的系统级问题上。 动态状态估计中的EKF与UKF对比：从原理到Matlab仿真
1. 动态状态估计为什么传统的静态估计在电网调度中开始“跟不上”了1.1 从断面快照到动态递推状态估计的两次思路跃迁做电力系统状态估计的同行都知道传统调度自动化系统里跑的基本都是静态状态估计最常见的就是加权最小二乘WLS。它的思路是拿当前一个时间断面的所有量测构建一组方程迭代求解出系统各节点的电压幅值和相角。好处是理论成熟、程序健壮实际运行了这么多年稳定可靠。但它的本质是一个“静态拍照”过程——拿到一张断面图算出此刻系统最可能的状态然后这张图就作废了下一个时刻又重新拍一张。这种思路在稳态工况下没太大问题但一旦系统进入动态过程问题就来了。比如新能源大规模接入之后出力波动频繁母线电压和线路功率的短时变化比过去剧烈得多又比如发生扰动或故障后的暂态过程中功角、角速度这些量以秒甚至毫秒级的尺度在演化这时候静态估计的“断面快照”模式就明显不够用。它没有利用系统自身的演化规律每一拍都从零开始算既浪费了历史信息也对量测误差的抑制不够充分。动态状态估计Dynamic State Estimation, DSE则是换了一个思路把电力系统当成一个动态系统来对待用状态方程描述它的演化规律再用量测方程把状态和观测量联系起来。每一时刻的估计结果一方面来自上一时刻状态的预测递推另一方面来自当前时刻的量测修正两者按各自的置信度加权融合。这就有点像开车时导航的做法车本身的惯性模型告诉你“我大概开到哪了”GPS量测告诉你“我实际观测到在哪”两者结合出来的位置远比单独依赖某一个要准。正是这个“预测修正”的递推结构让动态状态估计具备了静态估计没有的两个能力一是能输出状态量的预测值这在预警和稳定控制中非常有用二是能在量测质量差、甚至短暂缺失时依靠动态模型维持一段时间有效估计。EKF和UKF就是在这条思路上最常用的两种非线性滤波实现手段。1.2 动态状态估计的数学框架状态方程、量测方程与噪声建模要把动态状态估计讲清楚首先得把它的数学表达式立起来。电力系统中我们关心的是发电机的功角δ、角速度ω这类动态状态变量。以最经典、也是最常用来做滤波算法验证的发电机二阶模型摇摆方程为例状态方程δ(k1) δ(k) Δω(k)·Δt w₁(k)Δω(k1) Δω(k) (Δt/M)·[P_m - P_e(δ(k)) - D·Δω(k)] w₂(k)其中M是发电机惯性时间常数D是阻尼系数P_m是机械功率P_e是电磁功率。电磁功率和功角的关系是非线性的P_e(δ) (E′_q·V_s / X′_d)·sin(δ)这里的E′_q是发电机暂态电动势V_s是无穷大母线电压X′_d是暂态电抗。是不是已经闻到非线性的味道了对功角和功率之间差着一个正弦函数这个非线性是躲不掉的。量测方程则根据实际装设的PMU量测来定常见的有z₁(k) δ(k) v₁(k)z₂(k) Δω(k) v₂(k)z₃(k) P_e(δ(k)) v₃(k)w和v分别是过程噪声和量测噪声工程上通常假设为相互独立的高斯白噪声协方差矩阵分别为Q和R。动态状态估计的任务就是基于这些带噪声的量测序列实时估计出真正的状态量δ和ω。刚才说的是单机无穷大OMIB系统的模型。扩展到多机系统时状态向量会变成所有发电机的功角和角速度拼在一起的大向量量测向量也相应变成多个PMU断面的拼接但整体的“状态方程量测方程”框架完全一致。这也是为什么读懂单机系统的滤波流程之后往多机系统扩展是水到渠成的事。1.3 非线性是动态估计绕不开的坎讲到这里有一个问题自然而然浮出来既然状态方程中已经有δ和ω的线性递推关系为什么不能直接套用线性卡尔曼滤波答案就出在量测方程上。P_e(δ)和δ之间的正弦关系、以及多机系统中更复杂的功率方程都是不可忽略的非线性函数。如果把非线性直接忽略用线性卡尔曼滤波硬算状态估计结果大概率要偏离严重时滤波器直接发散。所以问题的核心就变成在非线性条件下如何完成卡尔曼滤波框架中的预测步和更新步EKF和UKF给出的答案代表了两种典型的技术路线EKF的思路是“把非线性问题线性化”——对状态方程和量测方程在当前估计点附近做一阶泰勒展开然后套用标准卡尔曼滤波公式UKF的思路是“对概率分布做采样逼近”——用一组精心挑选的sigma点去逼近状态分布经过非线性变换后的真实均值和协方差完全绕开了雅可比矩阵的计算。这两种思路各自的适用场景、实现难度、数值稳定性都有差异下面两章我展开讲。2. EKF与UKF处理非线性滤波问题的两条路线2.1 先从标准卡尔曼滤波说起为什么它不能直接用在电力系统里要说清EKF和UKF的差异必须先回到标准卡尔曼滤波本身。标准KF针对的是线性系统x(k1) F·x(k) w(k)z(k) H·x(k) v(k)其中F和H都是常矩阵。它的核心思想特别朴素用状态方程做一步预测得到状态预测量x_pred和预测协方差P_pred然后计算卡尔曼增益K用量测残差对预测值做修正。整个过程只需要矩阵乘法和求逆运算计算量小而且在线性高斯条件下是理论最优估计器。问题是电力系统的量测方程根本不是线性的。H不是一个常数矩阵而是状态量的函数H(x)。这时候如果依旧套用线性KF的公式预测协方差和量测协方差的传递就会出现系统性偏差。事实上我见过不少初学者把P_e(δ)和δ之间的关系做线性近似直接用常数增益算结果滤波轨迹和真值之间总是隔着一道“鸿沟”怎么调Q和R都补不回来。本质上就是因为丢弃了非线性项带来的信息。正因为如此工程中必须引入专门处理非线性滤波的算法。EKF和UKF的出发点相同都是把标准KF的思想推广到非线性场景但推广的手段南辕北辙。2.2 EKF一阶泰勒展开思路简单但代价是精度和稳定性的妥协扩展卡尔曼滤波的做法是在每个时刻把非线性函数在当前状态估计点附近做一阶泰勒展开。具体来说对于状态方程x(k1) f(x(k)) w(k)会计算状态转移雅可比矩阵F(k) ∂f/∂x | x x̂(k)对于量测方程z(k) h(x(k)) v(k)计算量测雅可比矩阵H(k) ∂h/∂x | x x_pred(k)有了这两个雅可比矩阵就可以套用和标准KF几乎一样的形式预测步x_pred(k1) f(x̂(k))P_pred(k1) F(k)·P̂(k)·F(k)ᵀ Q更新步K(k1) P_pred(k1)·H(k1)ᵀ·[H(k1)·P_pred(k1)·H(k1)ᵀ R]⁻¹x̂(k1) x_pred(k1) K(k1)·[z(k1) - h(x_pred(k1))]P̂(k1) [I - K(k1)·H(k1)]·P_pred(k1)单看公式EKF确实很好理解代码也不难写。但它的短板同样明显第一一阶泰勒展开本身就是一种近似只保留了函数的线性主部丢掉了二阶及以上的信息。当系统非线性程度较强时比如功角在暂态过程中大幅摆动、P_e(δ)的导数变化剧烈这种一阶近似带来的误差会累积严重时导致滤波发散。第二雅可比矩阵的计算是个麻烦事。对单机二阶模型手动推导F和H还算轻松一旦扩展到多机系统、状态量几十维、量测方程包含各种功率表达式时手动推导雅可比矩阵就是一场灾难。很多人会改用数值差分来替代解析求导但数值差分本身又引入步长选择的困扰步长太大误差大步长太小数值不稳定。第三EKF对初值和噪声统计的敏感性比较高。初值偏差太大时线性化点离真实点太远一阶近似的有效性大打折扣滤波器容易在启动阶段就“带偏”后面再怎么调都很难拉回来。那是不是说EKF就没用了不是。在弱非线性、且对计算速度要求高的场景下EKF依然是性价比很高的选择。毕竟它每次递推只需要做两次矩阵乘法、一次矩阵求逆计算开销远小于后面要讲的UKF。如果系统模型本身特性比较“温和”状态演化始终在平衡点附近小范围波动那么EKF的一阶近似误差完全在可接受范围内。2.3 UKF不碰雅可比用sigma点直接逼近概率分布无迹卡尔曼滤波换了一个完全不同的角度切入。它不去线性化非线性函数而是去逼近状态分布经过非线性变换后的统计特性。这个思路的实现载体就是无迹变换Unscented Transform, UT。UT的核心思想是对于n维状态向量x取其均值点x̄再加上2n个sigma点这些点按照一定规则分布在均值周围并且被赋予不同的权重。然后把每个sigma点分别通过非线性函数传播从传播后的sigma点集合中加权计算输出状态的均值和协方差。这样做的好处是不需要计算任何雅可比矩阵只需要把非线性函数当“黑箱”调用若干次就能以较高的精度逼近真实分布。sigma点选取的经典公式如下χ₀ x̄χᵢ x̄ (√((nλ)·P))ᵢi 1, ..., nχᵢ₊ₙ x̄ - (√((nλ)·P))ᵢi 1, ..., n对应的权重为W₀ᵐ λ/(nλ)W₀ᶜ λ/(nλ) (1 - α² β)Wᵢᵐ Wᵢᶜ 1/[2(nλ)]其中λ α²(nκ) - nα是一个控制sigma点散布范围的参数通常取一个较小的正数比如1e-3κ是次级缩放参数通常取0β用来融入分布先验信息高斯分布时取2。UT完成之后后续的预测步和更新步就是在sigma点层面上展开。预测步中把sigma点通过状态方程传播加权算出状态预测均值和预测协方差更新步中把预测的sigma点通过量测方程传播算出量测预测均值、量测协方差以及状态和量测的互协方差然后计算卡尔曼增益并更新状态。整个过程全部避开了雅可比矩阵算法对非线性的适应能力大幅度提升。理论上UT变换可以精确捕捉到非线性函数传播后的二阶矩信息在精度上相当于二阶泰勒展开而EKF只保留了一阶项这正是UKF在强非线性场景下误差更小的根本原因。不过UKF也不是没有代价。它的计算量明显大于EKF——每步递推需要2n1次非线性函数传播状态维数越高sigma点越多计算开销越大。另外sigma点选取涉及多个超参数虽然经验默认值能覆盖大多数情况但真到具体问题时也需要根据状态尺度做适当调整。2.4 EKF vs UKF一张表看透两者的实质差异我把两种算法放在一起做对比。这张表是我在实际项目里反复对照后梳理出来的比单纯看公式直观得多对比维度EKFUKF核心思想将非线性函数在估计点做一阶泰勒展开用sigma点逼近非线性变换后的概率分布是否需要雅可比矩阵需要手动推导或数值差分不需要对非线性的适应能力弱只保留一阶项强相当于二阶精度计算复杂度低适合高维大系统在线计算偏高需2n1次非线性传播对初值偏差的敏感性敏感初值差容易发散更稳健sigma点覆盖范围大实现难度公式直观但雅可比推导繁琐无需求导但参数和采样规则需要理解适用场景弱非线性、计算资源紧张、在线实时强非线性、对精度要求高、离线或在线均可用注意这张表不是绝对的。实际项目中我见过用EKF跑得稳如泰山的多机系统也见过UKF因为超参数没调好而对量测噪声特别敏感的反例。核心结论是算法的选择必须结合具体系统的非线性程度和算力约束不能简单说谁一定更好。3. Matlab仿真框架搭建与核心代码拆解3.1 整体架构设计模块化思路让滤波对比更清晰在用Matlab实现EKF和UKF之前我建议先想清楚代码的整体组织方式。我看到不少初学者把滤波算法、系统建模、结果绘图全塞在一个脚本里几百行代码铺下来出问题的时候根本不知道是模型写错了还是滤波逻辑写错了。正确的做法是模块化拆分。我给一个建议的工程目录结构dse_project/ ├── main_RunDSE.m # 主入口设置参数、调用滤波算法、绘制结果 ├── system_model.m # 状态方程与量测方程定义 ├── generate_true_states.m # 生成真值轨迹模拟系统的真实动态 ├── generate_measurements.m # 基于真值叠加噪声生成量测序列 ├── ekf_estimation.m # EKF滤波主函数 ├── ukf_estimation.m # UKF滤波主函数 ├── compute_rmse.m # 计算均方根误差 └── plot_results.m # 绘图脚本main脚本负责把其他函数串起来这样换系统参数、换滤波算法、换初始条件都很方便。比如想对比EKF和UKF只需要在main里分别调用两个函数把结果存下来即可。想换系统模型只改system_model.m一个文件其他函数无需动。这种结构的另一个好处是方便查错。滤波代码里的一个符号错误可能只在特定参数组合下才暴露如果全堆在一个脚本里排查的难度会几何级上升。3.2 系统建模单机无穷大系统的状态方程与量测方程实现为了把代码讲透我用单机无穷大OMIB系统作为演示对象。它虽然简化但包含了动态状态估计的全部核心要素而且状态量只有二维结果可视化非常清楚。先把模型参数列出来% 系统参数 M 10; % 惯性时间常数秒 D 1.5; % 阻尼系数 Xd_prime 0.3; % 暂态电抗标幺值 Eq 1.2; % 发电机暂态电动势标幺值 Vs 1.0; % 无穷大母线电压标幺值 Pm 1.0; % 机械功率标幺值 dt 0.01; % 采样间隔秒 T 20; % 仿真时长秒状态向量x [δ, Δω]ᵀ其中δ是功角弧度Δω是角速度偏差标幺值ω - ω_s。状态方程和量测方程在Matlab里定义如下function [x_next, y] system_model(x, dt, params) % 状态方程 f(x) delta x(1); domega x(2); Pe (params.Eq * params.Vs / params.Xd_prime) * sin(delta); x_next zeros(2,1); x_next(1) delta domega * dt; x_next(2) domega (dt / params.M) * (params.Pm - Pe - params.D * domega); % 量测方程 h(x)这里选了3个量测功角、角速度和电磁功率 y zeros(3,1); y(1) delta; y(2) domega; y(3) Pe; end这里有个细节值得注意量测方程中同时包含直接测量量功角、角速度和间接测量量电磁功率。加入功率量测的目的是让滤波器有更多信息来校正状态。PMU可以直接测功角但由于相量测量单元本身的算法和通信环节实际测到的功角依然有噪声角速度和功率量测提供了不同的误差特性多量测融合能显著提升估计精度。3.3 真值数据生成与量测模拟仿真研究的必备一步在仿真验证滤波器性能时我们不能直接拿真实电网数据开跑而是先构造一个“已知答案”的场景生成一条确定性的状态演化轨迹作为真值再叠加噪声模拟量测。这样后续才能定量地比较滤波结果和真值的偏差。真值轨迹的生成代码很简单N round(T / dt); x_true zeros(2, N); y_true zeros(3, N); x_true(:,1) [0.3; 0.0]; % 初始功角0.3rad角速度偏差0 for k 1:N-1 [x_next, y] system_model(x_true(:,k), dt, params); x_true(:,k1) x_next; y_true(:,k) y; end y_true(:,N) system_model(x_true(:,N), dt, params);这里设定了一个初值功角0.3rad角速度偏差0。系统会从这组初值开始演化因为机械功率和电磁功率不平衡功角会出现一段动态过渡过程。这段动态过程恰恰是检验滤波器跟踪能力的“试金石”。量测数据则在真值上加高斯白噪声R diag([1e-4, 1e-4, 1e-3]); % 量测噪声协方差 z_meas y_true sqrt(R) * randn(3, N);注意这里R的含义功角量测噪声的标准差为0.01弧度约0.57度角速度为0.01标幺值功率量测噪声标准差约为0.032标幺值。这个量级比较贴近PMU的实际测量误差范围既不会太理想化也不会让滤波变得过于困难。3.4 EKF核心实现预测-更新循环与雅可比矩阵的处理方式EKF的实现在原理上很直接但实际写代码时有两个容易出错的点一是雅可比矩阵的计算二是预测步和更新步中雅可比矩阵所使用的状态点不同预测步用的是上一时刻的估计值更新步用的是当前时刻的预测值。我用数值差分的方式计算雅可比矩阵这样耦合到多机系统时不需要手动推公式function F compute_F(x, dt, params) % 状态方程对x的雅可比矩阵数值差分法 n length(x); eps 1e-8; F zeros(n, n); [x_base, ~] system_model(x, dt, params); for j 1:n x_pert x; x_pert(j) x_pert(j) eps; [x_pert_out, ~] system_model(x_pert, dt, params); F(:, j) (x_pert_out - x_base) / eps; end end function H compute_H(x, dt, params) % 量测方程对x的雅可比矩阵数值差分法 m 3; n length(x); eps 1e-8; H zeros(m, n); [~, y_base] system_model(x, dt, params); for j 1:n x_pert x; x_pert(j) x_pert(j) eps; [~, y_pert] system_model(x_pert, dt, params); H(:, j) (y_pert - y_base) / eps; end endEKF主循环function [x_est, P_est] ekf_estimation(z_meas, x0, P0, Q, R, dt, params) N size(z_meas, 2); n length(x0); x_est zeros(n, N); P_est zeros(n, n, N); x x0; P P0; for k 1:N % 预测步 F compute_F(x, dt, params); [x_pred, ~] system_model(x, dt, params); P_pred F * P * F Q; % 更新步 H compute_H(x_pred, dt, params); [~, y_pred] system_model(x_pred, dt, params); S H * P_pred * H R; K P_pred * H / S; x x_pred K * (z_meas(:,k) - y_pred); P (eye(n) - K * H) * P_pred; x_est(:,k) x; P_est(:,:,k) P; end end这段代码需要注意的细节是计算H时用的是x_pred预测值而不是上一时刻的x。我刚学EKF的时候在这里栽过跟头把H和F用同一个点的雅可比矩阵算结果滤波结果始终差一口气。原理上预测步的线性化点应该是上一时刻的后验估计而更新步的线性化点应该是当前时刻的先验预测这个顺序不能反。3.5 UKF核心实现sigma点生成、权重计算与状态更新UKF的代码稍微复杂一些但只要理解了sigma点和权重的概念写起来并不困难。核心函数如下function [x_est, P_est] ukf_estimation(z_meas, x0, P0, Q, R, dt, params) N size(z_meas, 2); n length(x0); alpha 1e-3; kappa 0; beta 2; lambda alpha^2 * (n kappa) - n; % 权重 Wm zeros(2*n1, 1); Wc zeros(2*n1, 1); Wm(1) lambda / (n lambda); Wc(1) lambda / (n lambda) (1 - alpha^2 beta); for i 1:n Wm(i1) 1 / (2*(n lambda)); Wm(in1) 1 / (2*(n lambda)); Wc(i1) 1 / (2*(n lambda)); Wc(in1) 1 / (2*(n lambda)); end x_est zeros(n, N); P_est zeros(n, n, N); x x0; P P0; for k 1:N % 生成sigma点 sqrt_P sqrtm((n lambda) * P); chi zeros(n, 2*n1); chi(:,1) x; for i 1:n chi(:,i1) x sqrt_P(:,i); chi(:,in1) x - sqrt_P(:,i); end % 状态传播 chi_pred zeros(n, 2*n1); gamma_pred zeros(3, 2*n1); for i 1:2*n1 [chi_pred(:,i), gamma_pred(:,i)] system_model(chi(:,i), dt, params); end % 计算预测均值和协方差 x_pred zeros(n, 1); for i 1:2*n1 x_pred x_pred Wm(i) * chi_pred(:,i); end P_pred Q; for i 1:2*n1 diff chi_pred(:,i) - x_pred; P_pred P_pred Wc(i) * (diff * diff); end % 量测传播 y_pred zeros(3, 1); for i 1:2*n1 y_pred y_pred Wm(i) * gamma_pred(:,i); end P_yy R; P_xy zeros(n, 3); for i 1:2*n1 dy gamma_pred(:,i) - y_pred; dx chi_pred(:,i) - x_pred; P_yy P_yy Wc(i) * (dy * dy); P_xy P_xy Wc(i) * (dx * dy); end % 更新 K P_xy / P_yy; x x_pred K * (z_meas(:,k) - y_pred); P P_pred - K * P_yy * K; x_est(:,k) x; P_est(:,:,k) P; end end这段代码里我选择用sqrtm函数来求解矩阵平方根。对于协方差矩阵P理论上它是正定对称的sqrtm的结果也是对称的。不过在数值计算中如果P出现非正定这在EKF里比较容易见到sqrtm会报错或给出复数结果。实际使用时可以先用chol函数做Cholesky分解速度更快但对矩阵正定性要求更高。我更推荐在每一步更新后对P做一次对称化处理即P (P P) / 2避免累积的舍入误差破坏对称性。UKF相比EKF还有一个细节优势这里的P_pred和P_yy计算都直接用了Wc权重对差分外积做加权求和整个过程不需要任何求导。当系统模型改动时比如从二阶模型换成三阶模型只需要修改system_model函数UKF部分一个字符都不用动这对做研究和对比实验来说特别友好。3.6 评估指标RMSE才是衡量滤波效果的唯一标尺滤波器跑完了如何评价好坏总不能拿“看起来差不多”这种主观判断来下结论。业内最常用的指标就是均方根误差RMSERMSE sqrt((1/N) · Σ(k1 to N) (x̂(k) - x_true(k))²)对于多状态量可以分别计算每个状态的RMSE也可以做一个综合的加权平均。我的习惯是分开算因为功角和角速度的量纲不同数值尺度差异很大合在一起算会掩盖单个状态的性能差异。function rmse compute_rmse(x_est, x_true) diff x_est - x_true; rmse sqrt(mean(diff.^2, 2)); end这个函数简单得不能再简单但它是横向对比EKF和UKF性能的核心依据。在第四章的仿真结果中我会重点关注两个问题一是RMSE的数值二是滤波轨迹是否在暂态过程中出现明显偏差。4. 仿真结果对比EKF与UKF在同一测试场景下的实战表现4.1 测试场景设计初值偏差和暂态过程是检验滤波器最好的“试金石”为了让两种算法的差异真正显现出来我特意把测试场景设计得“刁钻”一些。首先滤波器的初始状态不设成真值而是设一个明显偏差的初值x0 [0.5; 0.05]; % 滤波器初值功角0.5rad角速度偏差0.05而真值是从功角0.3rad、角速度偏差0开始演化的。也就是说滤波器在启动瞬间就带着一个不小的偏差无论是EKF还是UKF都必须依靠量测信息一点点把状态“拉”回正确轨道。这个过程正好可以观察两种算法的收敛速度。过程噪声和量测噪声的设置也很有讲究Q diag([1e-6, 1e-4]); % 过程噪声协方差 R diag([1e-4, 1e-4, 1e-3]); % 量测噪声协方差Q的量级必须和模型的不确定性匹配。这里把功角的Q设得比角速度小一个量级是因为功角方程本质上是对角速度的积分积分本身会平滑高频噪声而角速度方程直接受到功率不平衡的影响不确定性更大。4.2 状态估计轨迹从图形上看两种算法的跟踪能力差异我跑完仿真之后的典型输出是三条曲线真值轨迹、EKF估计轨迹、UKF估计轨迹。在功角图中三条曲线在前0.5秒内会有明显分离因为滤波器正在从初始偏差中恢复。但之后的走势会有分化EKF的功角估计曲线整体能跟上真值的趋势但在暂态过程的快速阶段功角变化速率最大的区间会出现明显的滞后。这个滞后不是滤波器参数调得不好造成的而是EKF一阶线性化和真实非线性函数之间的偏差在起作用。具体来说当功角变化剧烈时P_e(δ)的导数和二阶导数都很大一阶泰勒展开丢掉的二阶项不可忽略导致预测步的误差偏大。UKF的曲线则更贴近真值尤其在暂态中期和后期两者的偏差肉眼几乎看不出来。原因不难理解UKF通过sigma点传播不需要在某个点做线性化近似它对强非线性函数的均值和协方差传播更准确。特别是在系统从暂态向稳态过渡的阶段UKF表现出更平稳的收敛过程没有EKF那种“過冲”现象。角速度的估计差异同样明显。角速度本身数值较小绝对值在零点零几的范围内波动EKF在快速变化区间的估计抖动明显更大UKF则平稳得多。这种差异在RMSE指标上会体现得非常直接。4.3 量化对比RMSE数值揭示的真相我这里给出一个具有代表性的RMSE对比结果具体数值会因噪声随机种子和参数设置而略有浮动但相对趋势是稳定的状态量EKF RMSEUKF RMSEUKF相对误差改善功角δ弧度0.01180.0062约47%角速度Δω标幺值0.00410.0023约44%从这个结果能看出来UKF在这个测试场景下的RMSE大致是EKF的一半左右。这个差距的绝对值很小但放到电力系统动态状态估计的语境里意义完全不同。功角RMSE从0.012弧度降到0.006弧度意味着对发电机状态不确定性的掌握精度提升了近一倍在功角稳定裕度判断中这会直接影响预警阈值设定。不过我要强调一个容易被忽略的点这个RMSE对比是在“量测信息不足”的时候最悬殊。如果我把量测噪声R调得更小意味着量测更可靠两种算法的差距会缩小反之如果R调大量测更不可靠模型预测的作用增强UKF在模型传播上的精度优势会更加放大。理解这个规律比单纯记住“UKF比EKF好”重要得多。4.4 发散风险与数值稳定性那些“看起来对”的代码是怎么翻车的在仿真过程中我遇到过几个典型问题写出来给各位避坑。第一个是EKF的协方差矩阵P出现非正定。根因是浮点舍入误差在多次递推中累积破坏了P的对称正定性。一旦P非正定后续计算里会蹦出负方差、甚至NaN滤波器直接报废。解决方法是每步更新后强制对称化并且可以用特征值分解把P的很小负特征值“截断”到零。UKF因为用sigma点传播协方差对P的正定性要求同样严格但我的实测经验是UKF出现非正定的概率比EKF低一些因为它没有显式的雅可比矩阵乘法。第二个是数值差分计算雅可比矩阵时的步长选择。我一开始用eps 1e-6结果在状态量数值本身就很小的情况下比如Δω只有0.01量级微扰步长和状态值处于同一量级差分的截断误差和舍入误差都快到极限了算出来的雅可比矩阵噪声很大。后来把步长调到1e-8效果明显改善。如果状态量尺度差异很大还可以逐状态设置差分步长。不过最稳妥的做法还是解析推导雅可比矩阵只是多机系统下推导工作量确实不小。第三个是UKF的参数α。我最早直接用默认的1e-3在单机系统上没问题但后来把算法搬到另一个状态量尺度差异更大的系统时发现sigma点生成的矩阵平方根在数值上出现退化。适当增大α到1e-2问题就消失了。这里的经验是α控制sigma点离均值的距离状态量尺度小的时候α太大会让sigma点覆盖范围过大状态量尺度大的时候α太小又会让sigma点挤在一起需要根据实际情况微调。5. 工程落地的实用建议与调参经验5.1 过程噪声矩阵Q的调参逻辑不是随便给个对角阵就行我在做动态状态估计项目时发现最容易被忽视的参数就是Q矩阵。很多初学者直接把Q设成一个很小的对角阵理由是“我认为模型很准”。但这样做往往会导致滤波器对量测的信任度过低量测信息没有充分发挥作用反过来Q设得太大滤波器又会过度相信量测导致估计结果在噪声中剧烈抖动。一个实用的调参方法是先根据系统模型的不确定性源来估计Q的下界。比如状态方程里忽略了多少高阶动态、模型参数M、D的辨识误差有多大、离散化带来的截断误差有多少把这些不确定性的方差加起来作为Q的初始估计。然后在这个基础上用一段真实或仿真数据反复试验观察滤波残差innovation的自相关特性如果残差存在明显的时间相关性说明Q设小了如果残差的方差远大于R的预测值说明Q设大了。对于单机系统我用的Q diag([1e-6, 1e-4])算是比较可靠的经验值。多机系统下各台发电机之间可能存在动态耦合更严谨的做法是建立完整的关联Q矩阵不过实际项目中为了简化通常只保留对角项。这种简化在大多数情况下是合理的只要量测冗余度足够滤波结果对Q的非对角项的敏感性不高。5.2 量测异常与时延动态估计在真实系统中面临的两大“暗礁”仿真环境里量测是理想化的——每个时刻都有数据、没有丢包、没有时延。但真实的PMU数据流远没有这么干净。我遇到过两种情况让滤波器几乎崩溃。一是量测野值。电网运行中PMU通道偶尔会跳出明显异常的数据点比如幅值突变到其他量级。如果滤波器直接把这个野值当作有效量测更新状态估计值会被瞬间“拽”到错误位置之后需要好几个周期才能恢复。工程上常用的处理办法是增加一个“创新限幅”逻辑计算量测残差zk - h(x_pred)如果残差超过某个阈值比如3倍量测噪声标准差就认为该量测是野值在本步更新中将其丢弃或降权。这个逻辑实现起来不难但对滤波器的稳定性帮助极大。二是量测时延。PMU从采样到数据上送中间要经过合并单元、相量数据集中器等多级环节时延在几十到几百毫秒不等。而动态状态估计的递推周期通常只有10到20毫秒如果不对时延做处理滤波算法拿到的“当前量测”其实是几十个周期前的数据状态已经往前跑了一截了。这种情况下一个常见的妥协做法是采用“延迟量测处理”在k时刻收到的是k-d时刻的量测时先用状态方程把滤波器的状态递推到当前时刻再用滞后量测做更新最后再递推回当前时刻。另一种更工程化的思路是在滤波器设计时预留一个量测缓冲并让滤波周期和PMU的数据周期保持一致从根本上规避时延问题。5.3 从单机到多机代码扩展路径和需要小心的坑单机无穷大系统的滤波跑通之后向多机系统扩展几乎是必然的路径。这个扩展在算法层面并不需要大改EKF和UKF的框架完全不变需要改的是系统模型状态向量从[δ₁, Δω₁]变成[δ₁, Δω₁, δ₂, Δω₂, ..., δₙ, Δωₙ]量测向量相应扩充到所有PMU的断面量测。量测方程会变成更复杂的网络功率方程这时EKF的雅可比矩阵就必须依赖数值差分或者专门的解析推导工具而UKF的优势更加凸显——只要系统模型函数能正确算出量测输出算法就能工作。一个经常被忽视的问题是发电机模型的阶数。前面用的是二阶经典模型只包含功角和角速度。如果要做更精确的动态估计发电机的励磁系统动态不能忽略这时需要引入三阶模型状态方程中再加入暂态电动势E′q的动态量测方程也要对应增加励磁电压相关的量测。状态维数从2维变成3维甚至更高UKF的sigma点数量随之增加计算量会明显上升。但如果你要估计的是功角稳定性相关的状态二阶模型往往已经够用需要分析励磁系统行为时再升级到三阶不必一开始就追求高保真模型。5.4 一个被低估的细节PMU量测的参考相角问题最后说一个我踩过坑、但很多文献里都不提的细节PMU测得的功角是相对时间同步信号的绝对相角而动态状态估计中的功角通常是相对系统参考机或者相对自身初始值的相对量。这两者之间如果不做换算滤波结果会产生一个常数偏置。更麻烦的是参考相角本身会随着系统频率的变化而漂移尤其是系统偏离额定频率时这个漂移会污染整个估计。实际处理办法是先对PMU功角量测做基准变换把绝对相角转换为相对系统参考轴的相对值。通常可以用系统惯量中心COI作为参考或者用一台参考机组的功角作为基准。这个预处理步骤看似不起眼但没有它滤波器的估计误差里会始终带着一个缓慢变化的系统误差RMSE指标再漂亮也掩盖不了这个隐患。5.5 我的总体体会没有“最好”的滤波器只有“适合”的滤波器经过这几次仿真对比我对EKF和UKF的选择有了一个更务实的判断。如果你要做的是数十台以上发电机的在线动态状态估计计算资源有限系统运行点相对平稳弱非线性假设大体成立EKF在今天依然是值得优先考虑的方案。它的公式结构清晰计算量小工程实现成熟出现问题时排查也容易。如果你更看重估计精度且系统可能经历较大的动态过程——比如暂态稳定分析、扰动后的轨迹估计——UKF的多花的那点计算成本完全值得。一次暂态过程的仿真时长也就十几秒哪怕UKF单步耗时是EKF的三倍总时长依旧在毫秒级这个代价换来的精度提升远远抵得过那点计算开销。在Matlab里把这两个算法放在同一个框架下对比测试是我觉得效率最高的做法。先跑通一个简单模型把滤波框架和评估流程固定下来再去扩展模型复杂度、替换数据源、调试噪声参数整个流程的迭代速度会快很多。这也是我写这篇内容的原因——希望更多人能快速迈过“滤波代码跑通”这道门槛把精力放在真正有意思的系统级问题上。