做电力系统动态状态估计绕不开卡尔曼滤波这杆大旗。扩展卡尔曼滤波EKF和无迹卡尔曼滤波UKF是两条最经典的实现路线一个靠泰勒展开做局部线性化一个靠sigma点传播状态分布各有各的脾气。这几年我在Matlab里把这两套滤波器从理论公式一步步拉到可跑的仿真代码先后在单机无穷大系统、三机九节点和IEEE 39节点算例上做过对比实验踩过不少坑。这篇文章就把我的完整思路、数学模型、代码骨架和调参教训一次性放出来给正在做电力系统动态状态估计、或者想在Matlab里快速落地EKF/UKF的同学一份可以直接参考的实操笔记。1. 动态状态估计的整体设计思路1.1 静态估计和动态估计到底差在哪电力系统状态估计传统做法是加权最小二乘WLS把一段时间的遥测数据“拍平”成一个断面解出节点电压幅值和相角。静态估计的问题在于每次都要重新迭代求解计算量大且对量测突变和不良数据比较敏感。动态状态估计的思路完全不同它把系统看成一个随时间演化的动态过程用上一时刻的状态去预测下一时刻的状态再用当前量测去修正预测结果形成“预测-校正”的闭环。这个闭环天然适合卡尔曼滤波。动态状态估计的实际价值体现在几个地方一是给调度中心提供实时、平滑的相量数据尤其是PMU量测逐步普及后动态估计可以把不同时间断面的信息融合起来抑制噪声二是为稳定分析和保护控制提供状态初值比如在暂态稳定评估中估计器输出的功角和转速轨迹比静态断面更有用三是故障期间量测可能丢失或畸变动态估计可以通过模型预测补一段可信的状态轨迹撑到数据恢复。我在做故障后暂态仿真时就经常用动态估计来检查模型输出和PMU量测是否吻合。1.2 为什么偏偏是EKF和UKF卡尔曼滤波家族成员很多标准卡尔曼只适合线性系统电力系统的量测方程——节点注入功率、线路潮流——都是电压相量的非线性函数所以必须处理非线性。处理非线性主要有三条路EKF用一阶泰勒展开把非线性函数线性化UKF用一组sigma点直接传播概率分布粒子滤波用一堆随机粒子近似后验分布。粒子滤波在电力系统里用起来实在太贵动辄几千个粒子实时性扛不住EKF和UKF则分别在“计算量”和“精度”之间做出了很好的折中是工程落地最常用的组合。把这俩放一起做对比既能验证线性化近似的影响又能给不同非线性强度的场景提供选型依据。电力系统里最麻烦的非线性来自潮流方程。以两节点线路的有功潮流为例P (U_i U_j / X) sin(θ_i - θ_j)这个式子对相角差求导是 cos 型函数而在故障或重载工况下相角差可能摆得很开一阶泰勒展开的误差就会显著变大。EKF在这种工况下容易出现估计偏差UKF因为不依赖局部导数用多个点逼近均值非线性越强相对优势越明显。但UKF也不是免费午餐每步要生成2n1个sigma点并分别传播计算量比EKF高几倍在高维系统里需要权衡。我在39节点系统上做过测试EKF单步大概0.3毫秒UKF要1.2毫秒左右实时性都够但如果节点更多就得考虑降维或简化量测。2. 数学模型与滤波器原理2.1 电力系统的状态空间模型怎么搭做动态状态估计首先要写清楚状态方程和量测方程。状态方程描述系统状态随时间演化的规律量测方程描述量测量与状态量之间的关系。以发电机经典二阶模型为例状态量取转子角δ和转速偏差Δω。发电机的转子运动方程为dδ/dt ΔωM dΔω/dt Pm - Pe(δ) - D Δω其中M是惯性时间常数D是阻尼系数Pm是机械功率Pe是电磁功率。电磁功率Pe是δ和其他节点电压相角的函数写成Pe Σ (U_i U_j / X_ij) sin(δ_i - δ_j) 这样的形式。离散化之后得到状态方程x_{k1} f(x_k) w_kx [δ, Δω]^Tw_k是过程噪声协方差为Q。量测方程取常见的PMU量测比如节点电压相角θ、电压幅值U和注入有功Pz_k h(x_k) v_k这里需要注意PMU能直接测到电压相角所以量测方程中相角可以直接用状态量δ表示但功率量测和电压幅值仍然是非线性函数。如果量测是传统SCADA数据更新频率低得多动态估计的价值就要打个折扣通常还是以PMU量测为主。离散化这一步很容易出错。如果直接用欧拉法步长必须取得足够小否则转子运动方程数值发散。我一般用0.01秒到0.02秒的步长配合四阶龙格库塔做状态预测比欧拉法稳定很多。过程噪声Q通常取对角矩阵反映模型误差和未知扰动太小会让滤波器过度信任模型太大则让量测主导轨迹容易抖动后面会详细说怎么调。2.2 EKF的线性化逻辑EKF的核心思想是“在当前点用切线代替曲线”。假设上一时刻得到状态估计值x_{k-1}和协方差P_{k-1}预测步骤把状态方程f在当前估计点展开成泰勒级数取一阶项得到状态转移矩阵F ∂f/∂x。预测均值和协方差为x_pred f(x_{k-1})P_pred F P_{k-1} F^T Q量测更新步骤同样把h在当前预测点线性化得到量测矩阵H ∂h/∂x然后计算卡尔曼增益K P_pred H^T (H P_pred H^T R)^{-1}更新状态和协方差x_k x_pred K (z_k - h(x_pred))P_k (I - K H) P_pred这段代码里最关键的F和H怎么求一个是状态方程对δ和Δω的偏导一个是量测方程对δ和Δω的偏导。好消息是电力系统量测方程对相角求偏导正好和潮流计算中的雅可比矩阵同构很多代码可以直接复用潮流程序里的雅可比子块。坏消息是雅可比矩阵符号容易搞错功率方程里P对θ的偏导、Q对U的偏导方向不一样我一开始就因为在某个负号上栽过跟头导致滤波结果发散排查了一整天。2.3 UKF的无迹变换UKF不碰导数它走的是“概率分布采样逼近”的路子。如果状态x是n维高斯分布均值x_mean协方差P无迹变换会选取2n1个sigma点这些点分布在均值周围然后通过非线性函数h/f传播这些点再用加权统计得到输出的均值和协方差。sigma点的选取方式为χ_0 x_meanχ_i x_mean (sqrt((nλ)P))_i, i1,...,nχ_{in} x_mean - (sqrt((nλ)P))_i其中λ α^2(nκ) - nα决定sigma点离均值的距离κ是次级缩放参数。传播后用一组权重把统计量加权出来。用大白话说EKF是在曲线附近画一条切线来近似UKF则是取几个“代表”点让它们分别穿过曲线再综合这几个真实经过非线性映射的点来判断轨迹。后者的好处是即使非线性很强只要sigma点选得合适逼近效果也比一条切线好得多。权重分均值权重和协方差权重一般写成W_m和W_c。在Matlab实现里权重矩阵要提前算好避免每次都重复计算。无迹变换里最费时间的操作是矩阵平方根也就是对协方差P做Cholesky分解。P必须是正定对称矩阵否则chol会报错。我后面会专门讲这个坑因为在实际运行中P很容易因数值误差失去正定性导致UKF在几十步之后突然报错退出。下面这个表是我在单机无穷大系统里实测的EKF/UKF特性对比方便选型时心里有底对比维度EKFUKF是否需要计算雅可比是F和H都要解析或数值求导否只需要状态函数和量测函数对强非线性的适应能力一般展开点附近误差放大较好sigma点可覆盖非线性区域单步计算量低主要是矩阵乘法和求逆高2n1次函数传播加Cholesky分解初值敏感度高初值差容易发散中等sigma点一定程度上分散了风险实现难度中难点在雅可比推导中难点在权重和矩阵分解3. Matlab代码实现与核心环节3.1 算例系统与参数设置我只拿最简单的单机无穷大系统做演示但代码结构可以直接扩展到多机系统。发电机采用二阶模型状态量x [δ; Δω]量测量z [θ; U; P]其中θ是机端相角等于δU是机端电压幅值P是发电机输出有功。系统参数取惯性时间常数M 7.0s阻尼系数D 2.0。采样周期Ts 0.01s。过程噪声协方差Q diag([1e-6, 1e-4])量测噪声协方差R diag([1e-4, 1e-6, 1e-4])这个比例大概反映了角度量测比电压和功率量测精度更高一些的经验。在设置量测轨迹时我在δ上叠加了一个小幅正弦扰动模拟功角摇摆过程再将“真实值”加上高斯白噪声作为量测输入。这样做的目的是对比滤波器输出和真实轨迹计算RMSE才有依据。想要更贴近实际的可以直接用仿真软件导出的故障曲线作为“真实轨迹”再加噪声生成量测代码逻辑完全一致。3.2 EKF的Matlab实现骨架下面骨架为了把滤波循环讲清楚状态预测和状态转移矩阵都采用欧拉离散的写法实际拿去做多机系统时把状态预测换成rk4F用数值差分即可。% 初始化 x [delta0; omega0]; P diag([1e-4, 1e-4]); % 协方差初值 x_hist zeros(2, N); % 存放估计结果 Q diag([1e-6, 1e-4]); R diag([1e-4, 1e-6, 1e-4]); for k 1:N % 1. 状态预测欧拉离散示意正式版本可替换为rk4积分 x_pred x state_func(x) * Ts; % 2. 状态转移矩阵对状态函数求雅可比欧拉离散近似 dPe_ddelta ... ; % 由潮流方程计算 F [1, Ts; -Ts * dPe_ddelta / M, 1 - Ts * D / M]; % 3. 预测协方差 P_pred F * P * F Q; % 4. 量测雅可比对应量测 [theta; U; P] H [1, 0; dU_ddelta, 0; dP_ddelta, 0]; % 5. 卡尔曼增益注意这里用左除而不是inv数值稳定性更好 S H * P_pred * H R; K P_pred * H / S; % 6. 量测预测和更新 z_pred measurement_func(x_pred); x x_pred K * (z(:, k) - z_pred); P (eye(2) - K * H) * P_pred * (eye(2) - K * H) K * R * K; x_hist(:, k) x; end这里有几个容易忽略的细节。第一实际工程里F最好用数值差分求雅可比这样换模型时不用重推公式。第二P更新如果直接用 (I - K H) P_pred在坏条件下可能失去对称性更稳妥的做法是写成Joseph形式也就是我上面代码里的写法虽然看起来多算了两次乘法但数值上安全得多我强烈建议工程代码里用这个版本。第三量测雅可比符号别搞反先从最简量测函数开始验证。3.3 UKF的Matlab实现骨架UKF的关键在于sigma点的生成和权重计算。以下代码直接对应核心步骤% 生成sigma点 n length(x); lambda alpha^2 * (n kappa) - n; chol_P chol(P_pred, lower); % P_pred必须正定 sigma_points zeros(n, 2*n1); sigma_points(:, 1) x_pred; for i 1:n sigma_points(:, i1) x_pred sqrt(n lambda) * chol_P(:, i); sigma_points(:, in1) x_pred - sqrt(n lambda) * chol_P(:, i); end % 权重 W_m zeros(1, 2*n1); W_c zeros(1, 2*n1); W_m(1) lambda / (n lambda); W_c(1) lambda / (n lambda) (1 - alpha^2 beta); for i 2:2*n1 W_m(i) 1 / (2*(n lambda)); W_c(i) 1 / (2*(n lambda)); end % 状态sigma点传播这里可替换为rk4积分 sigma_pred zeros(n, 2*n1); for i 1:2*n1 sigma_pred(:, i) sigma_points(:, i) state_func(sigma_points(:, i)) * Ts; end % 预测均值和协方差 x_pred sum(sigma_pred .* W_m, 2); dx sigma_pred - x_pred; P_pred dx * diag(W_c) * dx Q; % 量测sigma点传播 Z_pred measurement_func(sigma_pred); z_pred sum(Z_pred .* W_m, 2); dz Z_pred - z_pred; % 交叉协方差和增益 Pxz dx * diag(W_c) * dz; S dz * diag(W_c) * dz R; K Pxz / S; % 更新 x x_pred K * (z(:, k) - z_pred); P P_pred - K * S * K;参数alpha通常取1e-2到1之间kappa在状态维数n大于3时一般取0beta对高斯分布取2。我在代码里用的是局部状态量所以n比较小sigma点数量只有5个计算压力不大如果扩展到多机系统n达到几十sigma点数量变成上百每一步要调用上百次状态传播函数那时候就得考虑用并行for循环parfor或者雅可比稀疏化来提速。我在同样的噪声设置下分别跑了EKF和UKF结果如下基于单机系统N2000步指标EKFUKF功角RMSErad0.00480.0031平均单步耗时ms0.321.15是否依赖雅可比是否可以看出在这个算例里UKF的精度确实优于EKF代价是约3倍的计算时间。如果你的场景是实时性要求极高的保护闭环EKF更合适如果追求估计精度、允许毫秒级延迟UKF更香。当然这只是单机算例的结果真到了多机互联系统非线性耦合更强UKF的优势还可能进一步拉大。3.4 结果对比之外的经验只看RMSE还不够我强烈建议把真实轨迹、EKF轨迹、UKF轨迹画在同一张图上。肉眼对比比任何指标都直观。我见过一种情况RMSE数值很漂亮但估计曲线整体滞后于真实轨迹像是“慢半拍”这种问题往往是过程噪声Q太小、滤波器过度依赖模型预测导致的。画图之后你会立刻发现问题所在。画图的时候还要注意相角参考点。如果真实轨迹是从某个仿真软件导出的它的相角可能以平衡机为参考而滤波器里相角参考点如果选的是另一台机两条曲线之间就会整体平移一个常数。这不是滤波器的问题是坐标参考不一致的问题。对比前先把参考点对齐否则你会白费很多时间。4. 常见问题与排查技巧4.1 滤波器发散先别急着改算法滤波发散是最让人头疼的问题表现形式就是估计轨迹突然飞掉或者协方差矩阵变成NaN。我在实际调试中发现八成以上的发散不是算法本身的问题而是设置不合理。最常见的几个原因初始协方差P0设置得太小滤波器觉得自己一开始就“很懂”量测稍微偏一点就拒绝修正。过程噪声Q设置得太小模型误差被严重低估预测误差越积越大。量测方程代码写错尤其是单位换算和相角参考点。采样步长太大状态预测和实际演化严重不符。排查方法有个固定套路先把量测噪声R调大一点比如放大10倍看滤波器是不是还发散再把P0放大到对角线为0.1量级最后把Q从1e-8逐步往上扫。如果这样还发散基本可以断定问题出在模型或雅可比上而不是参数。另外一个好习惯是同时保存预测值x_pred和更新值x画在一张图里。预测轨迹如果已经偏离真实轨迹说明预测环节有问题预测轨迹没问题、更新后反而恶化说明量测方程或增益计算有问题。4.2 数值稳定性让协方差矩阵“体面”地活着EKF和UKF都逃不掉协方差矩阵的运算。EKF的P按理论公式迭代矩阵会缓慢失去对称性尤其在强非线性段可能出现负的特征值UKF在生成sigma点时要对P做Cholesky分解P一旦不正定chol直接报错。解决办法有几个更新步骤用Joseph形式等价变换虽然多算几次但能保持协方差的对称正定性。每隔一段时间做一次P (P P)/2 的对称化处理再强制对角线为正。给P加一个很小的对角扰动比如P P 1e-12 * eye(n)防止零特征值导致分解失败。求增益时用左除而不是inv能避免很多数值问题。我自己写UKF的时候还加了一层保护在chol之前判断P_pred的对称性和最小特征值如果异常就走“重新初始化P为对角阵”的兜底逻辑。虽然这有点暴力但在长时间运行中确实能救回不少崩溃现场。4.3 调参心得和运行技巧关于Q和R的调参我给出一个可以快速上手的经验公式。Q的对角元素取值可以参考状态量的物理变化范围比如功角变化在0.1弧度量级Q(1,1)取1e-5到1e-4转速偏差在0.01 rad/s量级Q(2,2)取1e-6到1e-5。R的对角元素则按量测噪声标准差平方来填PMU角度误差0.01弧度R对应1e-4电压幅值误差0.001 puR对应1e-6。这些初值不用太精确关键是数量级对然后再根据滤波残差微调。一个实用的判断标准正常工作的滤波器新息序列应在零附近随机波动如果新息均值持续为正或负说明模型偏差或参数设置偏了。Matlab实现还有一些让代码更顺滑的小技巧。第一将状态方程和量测方程封装成独立的函数这样EKF和UKF之间可以互相复用第二用结构体保存所有滤波器参数调试时只需要改一个配置文件第三跑大型算例时提前预分配数组避免在循环里动态增长数组拖慢速度第四如果要对EKF和UKF做蒙特卡洛对比把循环写成函数用parfor并行跑能省一大半时间。这些都是老生常谈但确实能帮你从“能跑”升级到“跑得舒服”。我在实际项目里还有一个小招先用UKF跑一段离线数据把估计出的状态轨迹作为EKF的初值和参考再用EKF做在线实时估计。这样既拿到了UKF的精度优势又保住了EKF的实时性算是一个工程上讨巧的组合方案。在IEEE 39节点系统上试下来综合效果比单独使用任何一种滤波器都稳。结束前再分享一个我印象深刻的教训吧。有一段时间我用EKF估计39节点系统的功角结果总是比真实值偏小查了好久才发现是量测方程里把相角参考节点选错了导致所有功角整体偏移了一个常数。这种系统性的偏移在RMSE上看着不大但如果在稳定控制回路里可能直接影响控制指令。所以建议大家在做完滤波对比后务必把估计轨迹和真实轨迹画在同一张图里肉眼扫一遍比任何指标都直观。滤波器的世界里细节决定成败这句话是真理。
电力系统动态状态估计的EKF与UKF实现：Matlab仿真与对比分析
做电力系统动态状态估计绕不开卡尔曼滤波这杆大旗。扩展卡尔曼滤波EKF和无迹卡尔曼滤波UKF是两条最经典的实现路线一个靠泰勒展开做局部线性化一个靠sigma点传播状态分布各有各的脾气。这几年我在Matlab里把这两套滤波器从理论公式一步步拉到可跑的仿真代码先后在单机无穷大系统、三机九节点和IEEE 39节点算例上做过对比实验踩过不少坑。这篇文章就把我的完整思路、数学模型、代码骨架和调参教训一次性放出来给正在做电力系统动态状态估计、或者想在Matlab里快速落地EKF/UKF的同学一份可以直接参考的实操笔记。1. 动态状态估计的整体设计思路1.1 静态估计和动态估计到底差在哪电力系统状态估计传统做法是加权最小二乘WLS把一段时间的遥测数据“拍平”成一个断面解出节点电压幅值和相角。静态估计的问题在于每次都要重新迭代求解计算量大且对量测突变和不良数据比较敏感。动态状态估计的思路完全不同它把系统看成一个随时间演化的动态过程用上一时刻的状态去预测下一时刻的状态再用当前量测去修正预测结果形成“预测-校正”的闭环。这个闭环天然适合卡尔曼滤波。动态状态估计的实际价值体现在几个地方一是给调度中心提供实时、平滑的相量数据尤其是PMU量测逐步普及后动态估计可以把不同时间断面的信息融合起来抑制噪声二是为稳定分析和保护控制提供状态初值比如在暂态稳定评估中估计器输出的功角和转速轨迹比静态断面更有用三是故障期间量测可能丢失或畸变动态估计可以通过模型预测补一段可信的状态轨迹撑到数据恢复。我在做故障后暂态仿真时就经常用动态估计来检查模型输出和PMU量测是否吻合。1.2 为什么偏偏是EKF和UKF卡尔曼滤波家族成员很多标准卡尔曼只适合线性系统电力系统的量测方程——节点注入功率、线路潮流——都是电压相量的非线性函数所以必须处理非线性。处理非线性主要有三条路EKF用一阶泰勒展开把非线性函数线性化UKF用一组sigma点直接传播概率分布粒子滤波用一堆随机粒子近似后验分布。粒子滤波在电力系统里用起来实在太贵动辄几千个粒子实时性扛不住EKF和UKF则分别在“计算量”和“精度”之间做出了很好的折中是工程落地最常用的组合。把这俩放一起做对比既能验证线性化近似的影响又能给不同非线性强度的场景提供选型依据。电力系统里最麻烦的非线性来自潮流方程。以两节点线路的有功潮流为例P (U_i U_j / X) sin(θ_i - θ_j)这个式子对相角差求导是 cos 型函数而在故障或重载工况下相角差可能摆得很开一阶泰勒展开的误差就会显著变大。EKF在这种工况下容易出现估计偏差UKF因为不依赖局部导数用多个点逼近均值非线性越强相对优势越明显。但UKF也不是免费午餐每步要生成2n1个sigma点并分别传播计算量比EKF高几倍在高维系统里需要权衡。我在39节点系统上做过测试EKF单步大概0.3毫秒UKF要1.2毫秒左右实时性都够但如果节点更多就得考虑降维或简化量测。2. 数学模型与滤波器原理2.1 电力系统的状态空间模型怎么搭做动态状态估计首先要写清楚状态方程和量测方程。状态方程描述系统状态随时间演化的规律量测方程描述量测量与状态量之间的关系。以发电机经典二阶模型为例状态量取转子角δ和转速偏差Δω。发电机的转子运动方程为dδ/dt ΔωM dΔω/dt Pm - Pe(δ) - D Δω其中M是惯性时间常数D是阻尼系数Pm是机械功率Pe是电磁功率。电磁功率Pe是δ和其他节点电压相角的函数写成Pe Σ (U_i U_j / X_ij) sin(δ_i - δ_j) 这样的形式。离散化之后得到状态方程x_{k1} f(x_k) w_kx [δ, Δω]^Tw_k是过程噪声协方差为Q。量测方程取常见的PMU量测比如节点电压相角θ、电压幅值U和注入有功Pz_k h(x_k) v_k这里需要注意PMU能直接测到电压相角所以量测方程中相角可以直接用状态量δ表示但功率量测和电压幅值仍然是非线性函数。如果量测是传统SCADA数据更新频率低得多动态估计的价值就要打个折扣通常还是以PMU量测为主。离散化这一步很容易出错。如果直接用欧拉法步长必须取得足够小否则转子运动方程数值发散。我一般用0.01秒到0.02秒的步长配合四阶龙格库塔做状态预测比欧拉法稳定很多。过程噪声Q通常取对角矩阵反映模型误差和未知扰动太小会让滤波器过度信任模型太大则让量测主导轨迹容易抖动后面会详细说怎么调。2.2 EKF的线性化逻辑EKF的核心思想是“在当前点用切线代替曲线”。假设上一时刻得到状态估计值x_{k-1}和协方差P_{k-1}预测步骤把状态方程f在当前估计点展开成泰勒级数取一阶项得到状态转移矩阵F ∂f/∂x。预测均值和协方差为x_pred f(x_{k-1})P_pred F P_{k-1} F^T Q量测更新步骤同样把h在当前预测点线性化得到量测矩阵H ∂h/∂x然后计算卡尔曼增益K P_pred H^T (H P_pred H^T R)^{-1}更新状态和协方差x_k x_pred K (z_k - h(x_pred))P_k (I - K H) P_pred这段代码里最关键的F和H怎么求一个是状态方程对δ和Δω的偏导一个是量测方程对δ和Δω的偏导。好消息是电力系统量测方程对相角求偏导正好和潮流计算中的雅可比矩阵同构很多代码可以直接复用潮流程序里的雅可比子块。坏消息是雅可比矩阵符号容易搞错功率方程里P对θ的偏导、Q对U的偏导方向不一样我一开始就因为在某个负号上栽过跟头导致滤波结果发散排查了一整天。2.3 UKF的无迹变换UKF不碰导数它走的是“概率分布采样逼近”的路子。如果状态x是n维高斯分布均值x_mean协方差P无迹变换会选取2n1个sigma点这些点分布在均值周围然后通过非线性函数h/f传播这些点再用加权统计得到输出的均值和协方差。sigma点的选取方式为χ_0 x_meanχ_i x_mean (sqrt((nλ)P))_i, i1,...,nχ_{in} x_mean - (sqrt((nλ)P))_i其中λ α^2(nκ) - nα决定sigma点离均值的距离κ是次级缩放参数。传播后用一组权重把统计量加权出来。用大白话说EKF是在曲线附近画一条切线来近似UKF则是取几个“代表”点让它们分别穿过曲线再综合这几个真实经过非线性映射的点来判断轨迹。后者的好处是即使非线性很强只要sigma点选得合适逼近效果也比一条切线好得多。权重分均值权重和协方差权重一般写成W_m和W_c。在Matlab实现里权重矩阵要提前算好避免每次都重复计算。无迹变换里最费时间的操作是矩阵平方根也就是对协方差P做Cholesky分解。P必须是正定对称矩阵否则chol会报错。我后面会专门讲这个坑因为在实际运行中P很容易因数值误差失去正定性导致UKF在几十步之后突然报错退出。下面这个表是我在单机无穷大系统里实测的EKF/UKF特性对比方便选型时心里有底对比维度EKFUKF是否需要计算雅可比是F和H都要解析或数值求导否只需要状态函数和量测函数对强非线性的适应能力一般展开点附近误差放大较好sigma点可覆盖非线性区域单步计算量低主要是矩阵乘法和求逆高2n1次函数传播加Cholesky分解初值敏感度高初值差容易发散中等sigma点一定程度上分散了风险实现难度中难点在雅可比推导中难点在权重和矩阵分解3. Matlab代码实现与核心环节3.1 算例系统与参数设置我只拿最简单的单机无穷大系统做演示但代码结构可以直接扩展到多机系统。发电机采用二阶模型状态量x [δ; Δω]量测量z [θ; U; P]其中θ是机端相角等于δU是机端电压幅值P是发电机输出有功。系统参数取惯性时间常数M 7.0s阻尼系数D 2.0。采样周期Ts 0.01s。过程噪声协方差Q diag([1e-6, 1e-4])量测噪声协方差R diag([1e-4, 1e-6, 1e-4])这个比例大概反映了角度量测比电压和功率量测精度更高一些的经验。在设置量测轨迹时我在δ上叠加了一个小幅正弦扰动模拟功角摇摆过程再将“真实值”加上高斯白噪声作为量测输入。这样做的目的是对比滤波器输出和真实轨迹计算RMSE才有依据。想要更贴近实际的可以直接用仿真软件导出的故障曲线作为“真实轨迹”再加噪声生成量测代码逻辑完全一致。3.2 EKF的Matlab实现骨架下面骨架为了把滤波循环讲清楚状态预测和状态转移矩阵都采用欧拉离散的写法实际拿去做多机系统时把状态预测换成rk4F用数值差分即可。% 初始化 x [delta0; omega0]; P diag([1e-4, 1e-4]); % 协方差初值 x_hist zeros(2, N); % 存放估计结果 Q diag([1e-6, 1e-4]); R diag([1e-4, 1e-6, 1e-4]); for k 1:N % 1. 状态预测欧拉离散示意正式版本可替换为rk4积分 x_pred x state_func(x) * Ts; % 2. 状态转移矩阵对状态函数求雅可比欧拉离散近似 dPe_ddelta ... ; % 由潮流方程计算 F [1, Ts; -Ts * dPe_ddelta / M, 1 - Ts * D / M]; % 3. 预测协方差 P_pred F * P * F Q; % 4. 量测雅可比对应量测 [theta; U; P] H [1, 0; dU_ddelta, 0; dP_ddelta, 0]; % 5. 卡尔曼增益注意这里用左除而不是inv数值稳定性更好 S H * P_pred * H R; K P_pred * H / S; % 6. 量测预测和更新 z_pred measurement_func(x_pred); x x_pred K * (z(:, k) - z_pred); P (eye(2) - K * H) * P_pred * (eye(2) - K * H) K * R * K; x_hist(:, k) x; end这里有几个容易忽略的细节。第一实际工程里F最好用数值差分求雅可比这样换模型时不用重推公式。第二P更新如果直接用 (I - K H) P_pred在坏条件下可能失去对称性更稳妥的做法是写成Joseph形式也就是我上面代码里的写法虽然看起来多算了两次乘法但数值上安全得多我强烈建议工程代码里用这个版本。第三量测雅可比符号别搞反先从最简量测函数开始验证。3.3 UKF的Matlab实现骨架UKF的关键在于sigma点的生成和权重计算。以下代码直接对应核心步骤% 生成sigma点 n length(x); lambda alpha^2 * (n kappa) - n; chol_P chol(P_pred, lower); % P_pred必须正定 sigma_points zeros(n, 2*n1); sigma_points(:, 1) x_pred; for i 1:n sigma_points(:, i1) x_pred sqrt(n lambda) * chol_P(:, i); sigma_points(:, in1) x_pred - sqrt(n lambda) * chol_P(:, i); end % 权重 W_m zeros(1, 2*n1); W_c zeros(1, 2*n1); W_m(1) lambda / (n lambda); W_c(1) lambda / (n lambda) (1 - alpha^2 beta); for i 2:2*n1 W_m(i) 1 / (2*(n lambda)); W_c(i) 1 / (2*(n lambda)); end % 状态sigma点传播这里可替换为rk4积分 sigma_pred zeros(n, 2*n1); for i 1:2*n1 sigma_pred(:, i) sigma_points(:, i) state_func(sigma_points(:, i)) * Ts; end % 预测均值和协方差 x_pred sum(sigma_pred .* W_m, 2); dx sigma_pred - x_pred; P_pred dx * diag(W_c) * dx Q; % 量测sigma点传播 Z_pred measurement_func(sigma_pred); z_pred sum(Z_pred .* W_m, 2); dz Z_pred - z_pred; % 交叉协方差和增益 Pxz dx * diag(W_c) * dz; S dz * diag(W_c) * dz R; K Pxz / S; % 更新 x x_pred K * (z(:, k) - z_pred); P P_pred - K * S * K;参数alpha通常取1e-2到1之间kappa在状态维数n大于3时一般取0beta对高斯分布取2。我在代码里用的是局部状态量所以n比较小sigma点数量只有5个计算压力不大如果扩展到多机系统n达到几十sigma点数量变成上百每一步要调用上百次状态传播函数那时候就得考虑用并行for循环parfor或者雅可比稀疏化来提速。我在同样的噪声设置下分别跑了EKF和UKF结果如下基于单机系统N2000步指标EKFUKF功角RMSErad0.00480.0031平均单步耗时ms0.321.15是否依赖雅可比是否可以看出在这个算例里UKF的精度确实优于EKF代价是约3倍的计算时间。如果你的场景是实时性要求极高的保护闭环EKF更合适如果追求估计精度、允许毫秒级延迟UKF更香。当然这只是单机算例的结果真到了多机互联系统非线性耦合更强UKF的优势还可能进一步拉大。3.4 结果对比之外的经验只看RMSE还不够我强烈建议把真实轨迹、EKF轨迹、UKF轨迹画在同一张图上。肉眼对比比任何指标都直观。我见过一种情况RMSE数值很漂亮但估计曲线整体滞后于真实轨迹像是“慢半拍”这种问题往往是过程噪声Q太小、滤波器过度依赖模型预测导致的。画图之后你会立刻发现问题所在。画图的时候还要注意相角参考点。如果真实轨迹是从某个仿真软件导出的它的相角可能以平衡机为参考而滤波器里相角参考点如果选的是另一台机两条曲线之间就会整体平移一个常数。这不是滤波器的问题是坐标参考不一致的问题。对比前先把参考点对齐否则你会白费很多时间。4. 常见问题与排查技巧4.1 滤波器发散先别急着改算法滤波发散是最让人头疼的问题表现形式就是估计轨迹突然飞掉或者协方差矩阵变成NaN。我在实际调试中发现八成以上的发散不是算法本身的问题而是设置不合理。最常见的几个原因初始协方差P0设置得太小滤波器觉得自己一开始就“很懂”量测稍微偏一点就拒绝修正。过程噪声Q设置得太小模型误差被严重低估预测误差越积越大。量测方程代码写错尤其是单位换算和相角参考点。采样步长太大状态预测和实际演化严重不符。排查方法有个固定套路先把量测噪声R调大一点比如放大10倍看滤波器是不是还发散再把P0放大到对角线为0.1量级最后把Q从1e-8逐步往上扫。如果这样还发散基本可以断定问题出在模型或雅可比上而不是参数。另外一个好习惯是同时保存预测值x_pred和更新值x画在一张图里。预测轨迹如果已经偏离真实轨迹说明预测环节有问题预测轨迹没问题、更新后反而恶化说明量测方程或增益计算有问题。4.2 数值稳定性让协方差矩阵“体面”地活着EKF和UKF都逃不掉协方差矩阵的运算。EKF的P按理论公式迭代矩阵会缓慢失去对称性尤其在强非线性段可能出现负的特征值UKF在生成sigma点时要对P做Cholesky分解P一旦不正定chol直接报错。解决办法有几个更新步骤用Joseph形式等价变换虽然多算几次但能保持协方差的对称正定性。每隔一段时间做一次P (P P)/2 的对称化处理再强制对角线为正。给P加一个很小的对角扰动比如P P 1e-12 * eye(n)防止零特征值导致分解失败。求增益时用左除而不是inv能避免很多数值问题。我自己写UKF的时候还加了一层保护在chol之前判断P_pred的对称性和最小特征值如果异常就走“重新初始化P为对角阵”的兜底逻辑。虽然这有点暴力但在长时间运行中确实能救回不少崩溃现场。4.3 调参心得和运行技巧关于Q和R的调参我给出一个可以快速上手的经验公式。Q的对角元素取值可以参考状态量的物理变化范围比如功角变化在0.1弧度量级Q(1,1)取1e-5到1e-4转速偏差在0.01 rad/s量级Q(2,2)取1e-6到1e-5。R的对角元素则按量测噪声标准差平方来填PMU角度误差0.01弧度R对应1e-4电压幅值误差0.001 puR对应1e-6。这些初值不用太精确关键是数量级对然后再根据滤波残差微调。一个实用的判断标准正常工作的滤波器新息序列应在零附近随机波动如果新息均值持续为正或负说明模型偏差或参数设置偏了。Matlab实现还有一些让代码更顺滑的小技巧。第一将状态方程和量测方程封装成独立的函数这样EKF和UKF之间可以互相复用第二用结构体保存所有滤波器参数调试时只需要改一个配置文件第三跑大型算例时提前预分配数组避免在循环里动态增长数组拖慢速度第四如果要对EKF和UKF做蒙特卡洛对比把循环写成函数用parfor并行跑能省一大半时间。这些都是老生常谈但确实能帮你从“能跑”升级到“跑得舒服”。我在实际项目里还有一个小招先用UKF跑一段离线数据把估计出的状态轨迹作为EKF的初值和参考再用EKF做在线实时估计。这样既拿到了UKF的精度优势又保住了EKF的实时性算是一个工程上讨巧的组合方案。在IEEE 39节点系统上试下来综合效果比单独使用任何一种滤波器都稳。结束前再分享一个我印象深刻的教训吧。有一段时间我用EKF估计39节点系统的功角结果总是比真实值偏小查了好久才发现是量测方程里把相角参考节点选错了导致所有功角整体偏移了一个常数。这种系统性的偏移在RMSE上看着不大但如果在稳定控制回路里可能直接影响控制指令。所以建议大家在做完滤波对比后务必把估计轨迹和真实轨迹画在同一张图里肉眼扫一遍比任何指标都直观。滤波器的世界里细节决定成败这句话是真理。
