做电力系统动态状态估计这些年EKF和UKF是我最常用的两个非线性滤波工具。今天这篇就来记录一下我用Matlab从零实现这两种滤波器并在IEEE标准节点系统上跑通全过程的思路、代码和踩坑记录。内容偏实操我会把模型怎么建、雅可比怎么求、Sigma点怎么采、参数怎么调一次性讲清楚适合刚接触动态状态估计的研究生也想给正在做PMU数据融合的工程师一点参考。先说结论EKF胜在计算量小、实现直观但强非线性场景容易因线性化误差发散UKF不用求导精度和鲁棒性都要好一截代价是多算2n1个Sigma点。真要在Matlab里跑两者代码量其实差不多真正花时间的是模型搭建和参数整定。1. 为什么要做电力系统动态状态估计1.1 静态估计的局限性和动态估计的必然性传统调度中心用的状态估计本质上是加权最小二乘(WLS)对一段时间的断面做拟合依赖SCADA系统慢速周期的数据通常几分钟才出一个结果。这个时间尺度应对系统稳定运行可以但遇到负荷快速波动、新能源出力随机变化或者故障后的暂态过程几分钟前的断面代表不了当前状态。更关键的是紧急控制决策需要毫秒级到秒级的状态感知SCADA数据根本喂不过来。所以有了基于PMU的动态状态估计。PMU以20到100帧每秒的速率同步采集电压电流相量配合发电机动态方程就能在线递推地跟踪功角、转速这些状态量。这个“动态”指的不是简单的时变参数而是把状态估计嵌入到状态方程里用上一时刻的估计值预测下一时刻再用实时量测修正形成闭环递推。这种思路其实就是卡尔曼滤波的经典框架。电力系统模型是非线性的所以标准卡尔曼滤波不能用得用非线性版本。EKF和UKF就是我在实际项目里用到最多的两个非线性卡尔曼变体。搞清楚这两个就能解决电力系统中绝大多数动态状态估计问题从单机无穷大系统到多机系统算法框架是通用的。1.2 EKF和UKF非线性卡尔曼滤波的两个主力EKF的思路很直接非线性函数在估计点附近做一阶泰勒展开用雅可比矩阵代替原本的线性状态矩阵。好处是只要会求导就能用计算效率高。坏处也明显线性化误差大尤其当系统强非线性或者初始误差大的时候协方差传播不准容易出现滤波发散。另外求雅可比矩阵容易出错尤其是多机电力系统里相互耦合的功率方程。UKF走的是另一条路用一组确定的Sigma点去近似状态分布这些点经过非线性函数传播后再按权重合成新的均值和协方差。整个过程完全不涉及求导非线性函数就是原样用的。所以UKF的处理精度在二阶以上比EKF高一阶对强非线性系统的适应能力也强。我在做IEEE 9节点系统时两个都实现了。要说哪个更好用我的体会是UKF更省心——不需要手推雅可比参数设置就alpha、beta、kappa三个数整定好就不太容易发散。但EKF也有存在价值状态维数比较高的时候UKF的Sigma点数量会膨胀计算量增长快而EKF计算量相对可控。所以实际选型往往是在精度和实时性之间做权衡不是无脑选UKF。2. 从原理到代码EKF和UKF的核心推导2.1 状态空间模型与离散化在做动态状态估计之前必须先确定被控对象的状态方程和量测方程。我用的模型是发电机二阶经典模型虽然简单但足以说清楚整个算法流程适合做教学和基准对比。状态量取功角δ和角速度ω标幺值形式输入是机械功率P_m输出是电磁功率P_e。连续时间模型如下dδ/dt ω - ω_s dω/dt (P_m - P_e - D*(ω - ω_s)) / (2H)其中D是阻尼系数H是惯性时间常数ω_s是同步转速。P_e是节点注入功率和系统其他节点电压、相角非线性相关。对于单机系统可以简化成P_e EV/X*sin(δ)E是暂态电动势X是电抗。这个式子是非线性的正好用来验证EKF和UKF。离散化我用的是一阶欧拉法步长dt取0.01秒对应PMU采样频率100Hz。离散后的状态方程delta(k1) delta(k) (omega(k) - omega_s)*dt omega(k1) omega(k) (P_m - P_e(k) - D*(omega(k)-omega_s))/(2H)*dt因为dt较小欧拉法精度够用。如果仿真步长更大可能要改用梯形法或Runge-Kutta。量测方程就直接取P_e(k)和V(k)如果可观测电压幅值。这里要特别提醒状态量和量测量的单位要统一用标幺值就别混入有名值不然卡尔曼增益会算得一团糟。2.2 EKF的线性化雅可比矩阵从哪来EKF的核心在于两个雅可比矩阵状态转移矩阵F和量测矩阵H。对于上面的离散方程F是状态方程对状态向量的偏导数。因为状态方程里omega(k1)含有P_e(k)而P_e是δ的函数所以F不是简单的常数矩阵需要逐项求偏导。以单机系统为例设状态x[delta; omega]量测zP_e雅可比矩阵H就是H [dP_e/ddelta, 0]其中dP_e/ddelta EV/X*cos(delta)。如果量测还包含电压幅值V那V对delta、omega也有偏导这里不展开。实际写Matlab代码时我建议先用符号运算做验证再手写解析式。因为手写偏导容易漏项尤其多机系统里功率对相角的偏导直接用B矩阵和G矩阵表示时更容易写错。我踩过的坑就是在推导H矩阵时把B和G的位置搞反导致滤波一会就发散。EKF主循环的逻辑很简单先根据当前状态和协方差做预测再用量测更新。预测部分用状态方程递推协方差用F、Q传播更新部分用H算卡尔曼增益然后修正状态和协方差。代码结构后面会给出。2.3 UKF的Sigma点采样不求导也能传播不确定性UKF的思想是无迹变换。假设状态x的均值为x_mean协方差为P那么构造2n1个Sigma点其中n是状态维度常见生成方式X_0 x_mean X_i x_mean (sqrt((nlambda)*P))_i , i1..n X_i x_mean - (sqrt((nlambda)*P))_i , in1..2n这里的lambda alpha^2*(nkappa) - nalpha控制Sigma点的散布程度kappa是次级缩放参数。sqrt((nlambda)*P)表示对矩阵做Cholesky分解得到的下三角矩阵的每一列。这组Sigma点通过非线性函数传播后按权重计算新的均值和协方差。权重定义W_m_0 lambda/(nlambda) W_c_0 lambda/(nlambda) (1 - alpha^2 beta) W_m_i W_c_i 1/(2*(nlambda)), i1..2nbeta是用来包含状态分布先验信息的参数对高斯分布取2是最优。UKF代码比EKF长一些但胜在不需要求雅可比直接把你写好的状态方程和量测方程当黑盒用就行。这也意味着如果以后换了更复杂的发电机模型只要替换函数句柄滤波框架完全不用动。3. Matlab实现的关键环节与实操步骤3.1 仿真数据生成用MATPOWER搭IEEE 9节点系统真试验证算法我推荐用MATPOWER导入标准节点系统。IEEE 9节点系统是最经典的三机九节点数据量不大但足够体现多机交互的动态。MATPOWER里提供case9.m文件直接加载就能得到节点、线路、发电机的稳态参数。先做一次潮流计算得到系统稳态工作点作为动态模型的初值。接下来生成“真值轨迹”我在稳态工作点基础上设置一个扰动比如在某个节点加一个持续0.1秒的三相短路故障然后用数值积分ode15s或自己写的R-K跑一段时间得到功角、转速、节点功率的演化轨迹。这个轨迹就是真实状态用来生成量测数据和控制滤波效果评估。生成量测时我分别在功率量测和功角量测上叠加高斯白噪声。PMU的幅值测量误差可以按0.01到0.02标幺值算相角误差按0.01弧度算转换成协方差就能得到量测噪声矩阵R。不能用一个凭空的R最好结合PMU的出厂指标和实际测试结果来定。3.2 EKF代码实现与Q/R参数整定EKF代码不难写但有几个细节决定成败。先说初始化状态初值在真实初值附近加一个随机偏差模拟实际情况下只知道大概状态。协方差P0设为一个对角阵对角线可以取状态误差的平方。比如初值偏差0.05弧度P0对角元素就取0.0025。Q矩阵表示过程噪声协方差代表模型误差和扰动的不确定性一般也设成对角阵。我常用的Q整定经验是先用小量级比如1e-6到1e-4跑看滤波轨迹是否平滑如果估计结果过于依赖量测波动大就适当增大Q如果估计结果跟不上真实轨迹变化说明Q太小模型预测过于自信。这是一种“黑盒整定”的笨办法但在没有精确噪声统计的时候很实用。R则根据传感器精度来如果PMU标称误差是1%R对角元素就取(0.01)^21e-4。EKF主循环的Matlab伪代码如下for k 2:N % 预测 x_pred f(x_est(:,k-1), u, dt); F compute_F(x_pred); P_pred F * P_est(:,:,k-1) * F Q; % 更新 H compute_H(x_pred); K P_pred * H / (H * P_pred * H R); x_est(:,k) x_pred K * (z(:,k) - h(x_pred)); P_est(:,:,k) (eye(n) - K * H) * P_pred; endcompute_F和compute_H是手写的雅可比函数h是量测方程。注意每一步H要在预测点x_pred处计算不是在上一时刻的估计点算这是新手容易忽略的。3.3 UKF代码实现与参数调优UKF的实现可以抽象成几步给定当前均值和协方差生成Sigma点然后分别通过状态方程和量测方程传播再计算增益。Matlab里可以用函数句柄传入状态方程不用像EKF那样维护雅可比代码复用性更好。我贴上核心部分的伪代码for k 2:N % 生成Sigma点 [X_sig, Wm, Wc] generateSigmaPoints(x_est(:,k-1), P_est(:,:,k-1), alpha, beta, kappa); % 状态传播 X_sig_pred f(X_sig, u, dt); x_pred X_sig_pred * Wm; P_pred (X_sig_pred - x_pred) * diag(Wc) * (X_sig_pred - x_pred) Q; % 量测传播 Z_sig h(X_sig_pred); z_pred Z_sig * Wm; Pzz (Z_sig - z_pred) * diag(Wc) * (Z_sig - z_pred) R; Pxz (X_sig_pred - x_pred) * diag(Wc) * (Z_sig - z_pred); % 更新 K Pxz / Pzz; x_est(:,k) x_pred K * (z(:,k) - z_pred); P_est(:,:,k) P_pred - K * Pzz * K; endgenerateSigmaPoints里要对P做Cholesky分解所以P必须保持正定。前面提到alpha取1e-3时Sigma点离均值很近对于状态变化剧烈的场景反而容易低估不确定性我实际跑下来的经验是alpha取0.01到0.1之间比较稳kappa取0beta取2。如果Cholesky分解报错说明P矩阵不是正定可以在分解前加一个很小的对角阵比如1e-12*eye(n)手动保证数值稳定。4. 结果对比与性能分析4.1 精度对比RMSE指标怎么算评价滤波精度我一般用均方根误差(RMSE)来对比分别统计状态量功角和转速以及输出功率的估计误差。公式不复杂RMSE sqrt(mean((x_est - x_true).^2))。注意要排除前几个点的暂态我通常让滤波器跑0.5秒后再开始统计给滤波器一个收敛过程。用IEEE 9节点系统的三台发电机状态做对比在相同初值偏差、噪声和故障扰动下UKF的功角RMSE通常比EKF小30%到50%。原因还是EKF的线性化误差尤其在系统经历大的扰动、功角摆动幅度比较大的时段EKF的协方差传播精度不足修正效果打折扣。下面是我一次典型仿真中得到的数值具体数值因扰动设置略有差异但趋势一致状态量EKF RMSEUKF RMSE提升幅度功角δ弧度0.0210.011约47%转速ω标幺0.00380.0021约45%电磁功率Pe标幺0.0260.017约35%UKF在非线性更强的时段优势会更明显系统平稳时两者差别不大。4.2 计算效率和实时性对比精度之外计算效率也是选型的重要考量。我在同一台机器、相同步数和数据长度下分别跑了100次记录单步平均耗时。UKF生成Sigma点、传播每个点计算量大约是EKF的1.5到2倍。对于9节点系统状态维度n6三台发电机每台2个状态UKF需要生成13个Sigma点还不算太夸张。但如果把模型换成四阶甚至六阶状态维度上到20以上UKF的计算负担会明显上升。如果你要做实时在线估计且系统规模不大UKF完全可以跑在实时仿真器上。如果状态维度很大且模型有比较平滑的非线性特性EKF的计算优势就体现出来了。实际项目中也可以考虑“混合策略”系统运行平稳时用EKF省资源检测到大扰动时切到UKF提高精度这个思路后续可以扩展。4.3 噪声协方差Q/R的敏感性分析Q和R的匹配程度直接影响滤波质量。我把Q设成固定对角阵发现当Q过小时滤波器对新量测的信任度很低估计轨迹会“滞后”于真值当Q过大时状态被噪声污染估计曲线毛刺明显。R的设定错误也类似R过大滤波器过度相信模型无法跟踪突变R过小量测噪声会被放大估计值波动剧烈。建议在调试阶段先用一组名义Q/R跑出基准然后固定其中一个参数把另一个乘以系数0.1到10扫描一遍看RMSE变化。这个过程可以用脚本批量跑画出RMSE随参数变化的曲线帮助理解系统的灵敏方向。实测下来EKF对Q/R比UKF更敏感因为EKF的协方差传播不准如果R设小了卡尔曼增益会异常偏高极易发散UKF在这方面稳定得多这也是我推荐新手先用UKF起步的原因。5. 常见问题与排查技巧5.1 滤波器发散先查这几个地方我在调试EKF时遇到过最典型的问题就是滤波发散估计值飞掉曲线直接冲出坐标系。排查顺序很有讲究先看状态方程对不对再看量测方程对不对最后才怀疑参数问题。具体检查点包括初值是否落在可行域内δ不可能超出一堆电抗决定的稳定边界P0是否反映了真实的不确定性如果P0设得过大会让早期增益过大Q是否过小导致模型预测误差被低估R是否过小导致过度相信量测雅可比是否计算正确这个我专门用符号工具箱验证过。如果一时查不出问题可以在滤波前几步打印新息序列z - h(x_pred)的均值和协方差。新息均值不该持续偏大新息协方差应该和理论值大致匹配。如果新息序列呈强相关性或偏差很大多半是模型与量测不匹配。5.2 雅可比矩阵写错用符号工具自查EKF的雅可比矩阵是手写最容易出错的地方。多机系统的状态方程里P_e对每个发电机功角都有偏导交叉项特别容易漏。我的做法是先用Matlab Symbolic Toolbox把f和h的解析式写出来用jacobian函数求偏导然后对比手写的函数数值结果。比如随机撒10个状态点比较手写雅可比和符号雅可比的计算结果看到所有点误差小于1e-8再放心用。还有一个细节雅可比是离散化后的状态方程对状态的偏导不是连续微分方程的偏导。如果你直接用连续模型雅可比套到离散递推里步长一大会有明显的截断误差。实际调试中我发现步长0.01s时连续雅可比还能用但步长大于0.05s时误差就比较明显了最好自己对离散方程重新推导雅可比或者直接用数值逼近的雅可比但那会增加计算量。5.3 UKF三个参数怎么设最省心UKF参数里alpha、beta、kappa网上有各种说法。我结合自己的测试给出一个推荐起点alpha0.01beta2kappa0。这个组合在大多数电力系统动态估计场景下都是个不错的起点。alpha越小Sigma点越集中对于接近线性的区域精度高但对于强非线性可能低估不确定性alpha在0.01到0.1之间是一个折中区间。beta2是所有高斯分布场景下的理论最优没什么好调的。kappa取0保证半正定为3-n可能是另一个常用选择但很多文献实际测试对单峰高斯影响很小。如果Cholesky分解报错我的经验是把alpha调大一点或者给P加上一个很小的对角扰动不要轻易调kappa到负值因为负的kappa可能会导致权重为负进而协方差失去正定性。5.4 步长与采样周期不匹配动态状态估计对时间离散化很敏感。我的仿真里用dt0.01s做状态递推量测假设每0.02s来一个对应50Hz采样所以滤波器每两步更新一次。这样做的原因是模拟PMU实际采样间隔可能大于积分步长。如果直接把观测更新步长设成和状态递推步长一样结果看起来没问题但和真实场景有差距。另一种情况是PMU采样率很高量测每0.01s来一次但状态递推步长也需要0.01s这样倒简单。最怕的是把状态递推步长设得很大比如0.1s那么EKF和UKF的离散化误差都会显著增大估计性能变差。我的经验是状态递推步长至少要比系统主要动态时间常数小一个数量级电力系统机电振荡周期大约在0.5到2s之间所以0.01s到0.02s的步长是合适的。如果模型包含更快的电磁暂态那就要用更小的步长或者把状态增广这已经超出常规动态状态估计的范围了。如果读者做的是这部分可以留言交流我对电磁暂态下卡尔曼滤波也做过一些尝试但那是另一个话题了。最后再分享一个实用习惯。我写Matlab代码时会把EKF和UKF封装成两个通用函数输入是状态方程句柄、量测方程句柄、初始状态、协方差以及Q/R输出是状态估计序列。这样后续想换模型、换节点系统只需要改函数句柄和初始参数滤波器主体完全不用动。配合MATPOWER批量生成场景可以快速做多工况对比。这一套下来不仅省了重复写代码的时间也让我更专注于模型本身的问题。