电力系统动态状态估计：EKF与UKF在39节点系统上的实现对比

电力系统动态状态估计DSE这个方向EKF和UKF是最常被拿来做对比的两种方法但很多人学完原理之后卡在到底怎么在自己算例上把它跑起来。我之前在39节点系统上完整做过一轮验证和实现踩了不少坑也把两种滤波器的特性摸了个遍。这篇就把整个过程梳理出来从算法原理、系统建模、代码实现到结果对比再到参考文献梳理给你一条能直接复现的路线。尤其适合正在做动态状态估计课题的研究生以及想从静态估计转向动态估计的电力系统工程师。1. 状态估计从静态走向动态问题到底出在哪1.1 静态状态估计的局限一个断面的快照解决不了动态问题传统上电力系统状态估计都是静态的——调度中心每隔几秒或几分钟从SCADA系统拿到一堆遥测数据用加权最小二乘法或其他优化算法算出当前断面下所有节点的电压幅值和相角。这种方法在稳态工况下是非常成熟的工程上运行了几十年可靠性已经验证得很充分。但它的本质是一个截面估计默认系统处于稳态没有把系统随时间演化的动态规律放进去。一旦系统进入动态过程——比如发电机转子角摆动、线路跳闸后的暂态过渡过程、负荷大幅波动——静态估计给出的结果就只是若干个时间断面上的孤立快照你很难从中看出系统状态量比如转子角、角速度是怎么随时间变化的。我最早做短路故障分析的时候就特别有感触故障发生的那几秒钟转子角摆动的轨迹恰恰是判断系统稳定性的关键但静态估计器根本给不出这个过程。它只能告诉你每一秒的稳态电压相角大概是多少至于发电机的机电暂态过程完全覆盖不到。1.2 PMU出现之后动态状态估计有了工程基础动态状态估计Dynamic State Estimation, DSE其实是很自然的延伸——既然系统有已知的动态模型比如发电机转子运动方程量测有时间序列比如PMU以几十帧每秒的速率上报相量数据那就可以用状态空间模型加滤波器把状态量从断面估计升级为轨迹跟踪。这里要特别提一下PMU相量测量单元的作用。它有两大优势一是数据率高一般25帧/秒、50帧/秒甚至更高能捕捉到机电暂态过程二是量测自带时间标签和相角信息可以直接把发电机功角、母线电压相角放进量测方程里。可以说没有PMU动态状态估计的应用价值会大打折扣有了PMU这个方向才真正具备工程落地的条件。1.3 为什么对比EKF和UKF理论基础与实现复杂度都典型做动态状态估计滤波器是最核心的算法载体。EKF扩展卡尔曼滤波是处理非线性系统最经典的做法思路简单、实现门槛低几十年来在很多领域都有成熟应用。UKF无迹卡尔曼滤波从2000年前后开始流行它不依赖雅可比矩阵的推导用sigma点直接传递非线性函数的统计特性能处理更强非线性的情况。电力系统的发电机动态模型中状态方程和量测方程都存在明显的非线性——功率方程里有sin、cos项转子运动方程里有乘积项。所以自然会把EKF和UKF放在一起对比。39节点系统又是IEEE标准算例系统规模适中有10台发电机、46条支路既能体现多机系统的复杂性又不至于因为节点太多导致仿真负担过重拿来做算法验证非常合适。我自己的体会是单纯跑通一个滤波器不算难难的是把两个算法的差异、参数影响、量测配置影响都弄清楚。后面的篇幅就按这个逻辑来讲先把原理掰开再聊建模和实现最后看数据和文献。2. 先搞懂原理EKF的线性化逻辑与UKF的sigma点采样2.1 EKF的核心逻辑与误差来源扩展卡尔曼滤波的基本思路其实很朴素标准卡尔曼滤波只适用于线性系统那我们就把非线性系统在当前估计点附近做一阶泰勒展开用雅可比矩阵替代线性卡尔曼里的状态转移矩阵和量测矩阵剩下的递推流程保持不变。具体来说EKF的递推分为两步。预测步用非线性状态方程直接外推状态x_{k|k-1} f(x_{k-1})协方差通过雅可比矩阵F传播P_{k|k-1} F P_{k-1} F^T Q更新步计算量测雅可比H求卡尔曼增益K用实际量测残差修正预测状态x_k x_{k|k-1} K(z_k - h(x_{k|k-1}))更新协方差P这套逻辑非常清晰只要会求雅可比矩阵代码写得很快。但问题恰恰出在一阶泰勒展开上——如果系统的非线性程度比较高一阶近似只保留线性项忽略了高阶项那么均值和协方差的传播都会有截断误差。在强非线性场景比如故障后转子角大幅摆动、功率方程在特定运行点附近变化剧烈下EKF可能给出有偏的估计严重时甚至发散。我在39节点系统上做过一个不算极端的测试把单回线三相短路故障设置为0.1秒后切除故障后转子角摆幅超过40度。这种场景下EKF的估计误差明显比稳态时要大因为故障瞬间发电机非线性模型的泰勒展开误差被放大了。2.2 UKF的无损变换不求导直接用sigma点穿过非线性UKF的思路则完全不同。它不再做线性化而是用一组精心选取的确定性采样点sigma点来近似状态的统计分布。核心操作是无损变换Unscented Transform, UT。在n维状态空间中先取2n1个sigma点每个点有对应的权重。这些sigma点经过了状态方程的非线性映射之后再按权重重新合成均值和协方差。理论上UT变换至少能精确到二阶矩对非线性函数的逼近精度高于EKF的一阶线性化而且不需要计算任何雅可比矩阵。对于高斯分布的情况sigma点的选取方式通常是x^(0) x̄权重 W^(0)x^(i) x̄ (√((nλ)P))_ii1,...,nx^(in) x̄ - (√((nλ)P))_ii1,...,nλ α²(nκ) - n其中α控制sigma点离均值的距离一般设1e-3到1之间κ是次要尺度参数β用来合并状态的先验分布信息高斯分布下取2。在滤波递推中和EKF一样分预测、更新两步区别在于预测步要逐个把sigma点送入非线性函数再加权合成预测均值和协方差。更新步也需要对量测方程做同样的sigma点变换再计算增益。2.3 从电力系统角度看两个算法的适用性差异电力系统动态状态估计有一个明显特点状态方程和量测方程的非线性程度并不均匀。发电机转子运动方程在正常运行点附近相对平滑但如果把暂态过程算进去特别是故障期间和故障后最初几个周期非线性程度会迅速上升。EKF的优势在于实现简单、解释性强对弱非线性场景比如正常负荷波动下的跟踪已经足够而且每次迭代只需要算一次雅可比矩阵计算量和代码量都小。UKF的优势是精度高、对强非线性更稳健代价是要对2n1个sigma点分别做非线性传播计算量大约是EKF的两到三倍。在后面的39节点验证中你会发现这两种算法的差异不是谁完全碾压谁而是在什么条件下谁更合适。这也是我建议大家别只盯着精度指标的原因——工程上还要考虑计算负担和实现复杂度。3. IEEE 39节点系统的建模准备模型、数据与仿真平台选型3.1 39节点系统的基本情况与数据来源IEEE 39节点系统又称新英格兰系统New England System是动态稳定分析和状态估计算法验证中非常经典的算例。它包含10台发电机、39个母线和46条支路其中第39号母线带一个等效的外部电网模型视作平衡节点。我用的数据是从MATPOWER的case39中直接读取的。MATPOWER里的case39包含母线参数、线路参数、发电机参数和负荷参数用它的loadcase函数就能把一个结构体加载进工作区。注意MATPOWER里的case39主要提供潮流计算所需的数据结构它并不是专门为动态仿真准备的所以缺了不少动态参数比如发电机的惯性时间常数、阻尼系数、暂态电抗需要自己补。3.2 动态参数补充用经典二阶模型做机电暂态做动态状态估计最常用的是发电机经典二阶模型也叫摇摆方程模型它描述发电机转子角δ和角速度ω的动态行为对研究机电暂态已经够用。对于第i台发电机dδ_i/dt ω_i - ω0dω_i/dt (ω0 / (2H_i)) * (Pm_i - Pe_i - D_i * (ω_i - ω0))其中H是惯性时间常数D是阻尼系数Pm是机械功率Pe是电磁功率。在39节点系统中我需要为每台发电机补充这些参数。具体数值我没有用标准库里的固定值而是参考了几篇经典文献和BPA动态数据卡中类似机组的典型取值H在3~6秒之间D在1~3之间这样系统动态特性比较合理多机之间的振荡模式也更丰富。这里有个非常实用的建议把补好的动态参数单独存成一个表格文件不要揉在MATPOWER的case39里。因为后面你很可能想对比不同参数配置下的滤波性能单独维护一份参数表能省不少事。3.3 仿真平台与工具链MATLAB为主配合小步长时域仿真整个验证我是在MATLAB环境下完成的主要用三块东西MATPOWER用来读取case39数据做潮流计算获取初始稳态运行点。时域仿真脚本自己写的四阶龙格-库塔RK4积分器步长0.01秒根据发电机动态方程模拟系统受扰动后的真实轨迹。滤波器脚本EKF和UKF各写一个函数输入是带噪声的量测序列输出是状态估计序列和误差统计。时域仿真这一步是关键。我见过有些同学直接用电力系统仿真软件比如PSAT、PSS/E生成轨迹然后再把轨迹导入MATLAB做滤波。这个思路也能跑通但我更推荐自己写RK4积分脚本原因有两个一是你能完全掌握量测数据是怎么生成的自定义故障场景灵活二是滤波器的验证需要把真实状态作为参照自己生成轨迹就天然知道真值在哪误差统计更干净。我用的做法是先用MATPOWER算潮流得到初始转子角和机械功率初值然后在仿真脚本里构造扰动通常在线路21-22上设置一个三相短路0.1秒后切除故障对整个暂态过程做RK4积分得到所有机组转子角、角速度的真实轨迹再把其中的一部分量当作真实状态叠加上高斯白噪声作为量测交给滤波器。4. 手把手实现滤波过程状态方程、量测方程与参数整定4.1 状态变量与量测方程先定维度再写代码在39节点系统的10台发电机都参与估计的场景下每台机取两个状态量转子角δ和角速度ω状态向量的维度就是20维。这个维度大小适中EKF的雅可比矩阵是20×20UKF需要41个sigma点在MATLAB里跑起来都不会有计算压力。量测方程则需要建立状态量与量测值的映射关系。我用的量测主要来自两部分假设每台发电机出口处都装了PMU能直接测量该机组的功角δ相对参考相角和角速度ω这类量测非常直接。母线电压幅值和相角量测它们通过潮流方程即功率平衡方程与发电机功角关联写出来就是带sin、cos的非线性函数。这里要特别提醒一个细节δ的量测参考基准。发电机功角本身是相对参考机组的相角差所以在构造量测方程时必须明确参考节点的选取否则看起来是绝对角的量测实际上却含有参考节点的角度变化会造成很大的模型误差。我在最初跑的时候就在这里翻过车——前几步滤波结果看起来还挺正常但几十个步长之后误差就开始漂移最终定位到量测方程中参考角处理不当的问题。4.2 EKF/UKF递推的代码逻辑能跑通是最低标准滤波器的核心逻辑并不复杂我贴一段MATLAB风格的EKF递推伪代码方便你理解整体结构function [x_est, P_est] ekf_power(dx_func, h_func, F_func, H_func, z_seq, x0, P0, Q, R) n length(x0); x_est zeros(n, length(z_seq)); x_k x0; P_k P0; for k 1:length(z_seq) % 预测步 x_pred dx_func(x_k); % 状态方程外推 F F_func(x_k); % 状态雅可比 P_pred F * P_k * F. Q; % 更新步 H H_func(x_pred); % 量测雅可比 z_pred h_func(x_pred); K P_pred * H. / (H * P_pred * H. R); x_k x_pred K * (z_seq(:,k) - z_pred); P_k (eye(n) - K * H) * P_pred; x_est(:,k) x_k; P_est P_k; end endUKF的代码复杂度会高一点因为要写sigma点生成和权重计算的函数。我建议不要急着追求代码精简先把结构写清楚先写一个unscented_transform函数输入状态向量和协方差输出sigma点阵列和对应权重再在滤波主循环里分别调用状态转移和量测映射的无损变换最后合成预测与更新。这样调试起来方便哪一步出了问题也容易定位。4.3 参数整定与初始化最容易翻车的环节参数整定是我在整个验证流程中花时间最多的地方这里必须单独拿出来讲因为大部分人一开始都会在这里栽跟头。第一是初始协方差P0。状态量中δ的量级通常在零点几到几弧度ω的量级在0.999到1.001左右标幺值两者不在一个数量级。如果直接把P0对角线都设成同一个值比如0.1那ω维度的初始不确定性相对太大滤波器前几十个步长会明显震荡。我最后用的P0对角线大致是δ维度取1e-4ω维度取1e-6相当于假定初始状态有约0.57度的功角误差和约0.01标幺的转速误差。第二是过程噪声协方差Q。Q太小滤波器会过度信任模型而对量测反应迟钝Q太大估计会震荡甚至发散。这个参数没有标准答案和你的动态模型误差水平强相关。我的经验做法是先跑一次无噪声仿真把真实轨迹和模型外推轨迹的最大偏差统计出来按这个量级反推Q的对角线值再微调一到两个数量级。第三是量测噪声协方差R。PMU量测的典型幅值误差在0.01到0.02标幺左右相角误差在0.01到0.03弧度左右我直接把R按这个量级设定。别把R设太小否则滤波器对量测噪声过度拟合估计轨迹会出现明显的毛刺。这类毛刺在功角跟踪曲线上尤其容易看到。5. 结果对比EKF和UKF在39节点系统上的实际表现5.1 正常运行工况下的精度对比两者都能用EKF甚至更稳先看最简单的场景系统运行在稳态附近只有小幅随机负荷波动和量测噪声不设置严重故障。这种情况下状态方程的非线性程度很低EKF的一阶线性化误差很小。我统计了100秒仿真时间内的功角估计误差EKF的RMSE大约在0.15度左右UKF大约在0.12度左右。两者的差距并不大。而且在某些机组上EKF的估计轨迹甚至比UKF更平滑——这说明在小扰动场景下UKF的精度优势并不明显有时还因为sigma点采样的细微偏差而多出一点噪声。但这里有一个隐性差异值得注意UKF在稳态场景下的稳定性往往比EKF更好尤其是当量测出现短暂缺失时。我试过模拟一个PMU通道断点几秒钟的情况EKF在这段时间内由于一味依赖模型外推误差上涨得很快UKF的sigma点传播方式保留了更完整的分布信息误差上涨要缓慢一些。也就是说UKF在量测信息不足时的体力更好。5.2 故障扰动下的跟踪能力UKF的优势开始显现真正的差距在故障场景下拉开。我在线路21-22上设置了三相短路故障0.1秒后跳开故障线路系统进入剧烈的机电暂态过程。这个场景下转子角摆幅最大超过40度系统非线性表现得非常明显。从跟踪结果看EKF在第一波大幅摆动期间出现了较明显的相位滞后最大功角估计误差达到了1.8度左右而且误差曲线和真实摆动轨迹呈现出明显的同频振荡——这在物理上对应的是一阶线性化带来的相位误差。UKF的最大估计误差则控制在0.9度以内相位滞后明显更小整体估计轨迹与真值几乎重合。这个结果和理论预期一致UKF通过sigma点传播对非线性映射的逼近精度更高在强非线性区间的优势被充分放大了。如果你做的研究涉及故障暂态、失步振荡这类强非线性过程UKF显然是更稳的选择。不过我还要补充一个令人意外的发现在故障切除后的瞬间两个滤波器都出现了短时的误差尖峰UKF的尖峰宽度和幅度都要小一些。原因在于故障清除是个开关性事件系统方程本身发生了结构突变网络参数变化任何基于连续模型递推的滤波器都会有一瞬的不适应。这个尖峰只能靠后续几个量测点把它拉回来无法完全避免。5.3 计算代价与工程权衡精度差距够不够换两倍算力光看精度还不够工程上必须算计算代价的账。同样100秒、0.01秒步长的仿真UKF的总运行时间大约是EKF的2.4倍主要多出来的开销就在sigma点的非线性传播上。20维状态量、41个sigma点每个点都需要过一遍状态方程和量测方程这是UKF无法回避的成本。从工程落地角度看我的建议是分场景选算法稳态监测、负荷跟踪这类弱非线性应用EKF完全够用计算效率高代码维护成本也低故障后暂态分析、机电振荡辨识这类强非线性应用优先选UKF用算力换精度的收益是值得的。下面是我整理的一组典型结果对照表RMSE为全部10台机组的平均供你参考仿真场景EKF功角RMSE度UKF功角RMSE度UKF相对耗时稳态小幅波动0.150.122.4x负荷阶跃10%突变0.420.272.3x三相短路线路跳开0.950.532.4x需要说明的是这些数值和你的具体参数整定强相关重点看趋势而不是绝对值。如果你想复现建议先把稳态场景跑通确认滤波器不发散再逐步增加扰动强度。6. 代表性参考文献梳理与后续扩展思路6.1 参考文献清单标题里特别提到了有参考文献这也是很多同学卡住的地方——动态状态估计方向文献很多但刚入门容易迷失。我这里按系统模型基础、滤波算法基础、电力系统DSE应用三个方向列几篇我实际读过且觉得有用的T. Athay, R. Podmore, S. Virmani, A Practical Method for the Direct Analysis of Transient Stability, IEEE Transactions on Power Apparatus and Systems, vol. PAS-98, no. 2, 1979. —— 39节点系统动态数据的经典来源很多论文中的39节点动态参数都追溯到这里建议列入参考文献。E. A. Wan, R. van der Merwe, The Unscented Kalman Filter for Nonlinear Estimation, Proceedings of the IEEE Adaptive Systems for Signal Processing, Communications, and Control Symposium, 2000. —— UKF理论最经典的原始文献sigma点选取、权重推导都讲得很清楚。G. Valverde, V. Terzija, Unscented Kalman Filter for Power System Dynamic State Estimation, IET Generation, Transmission Distribution, vol. 5, no. 1, pp. 29-37, 2011. —— 直接把UKF用于电力系统DSE的较早论文量测方程构造和应用场景很有参考意义。E. Ghahremani, I. Kamwa, Dynamic State Estimation in Power System by Applying the Extended Kalman Filter with Unknown Inputs, IEEE Transactions on Power Systems, vol. 26, no. 4, pp. 2203-2213, 2011. —— 将EKF用于电力系统DSE的参考文献讨论了输入不确定性条件下的处理方式。J. Qi, K. Sun, J. Wang, H. Liu, Dynamic State Estimation for Multi-Machine Power System Using Unscented Kalman Filter, Proceedings of the IEEE Power and Energy Society General Meeting, 2013. —— 多机系统场景下的UKF实现和39节点系统场景很贴近。6.2 基于文献的后续扩展方向把这套EKF/UKF基础实现跑通之后你可以顺着文献往几个方向扩展每个方向都能形成新的研究角度一是在模型中加入未知输入。实际系统中机械功率Pm并不精确可知文献[4]的做法就是把Pm视作未知输入通过状态扩维或者估计器设计来同时估计状态和输入。这个方向对工程价值很大因为调速器、原动机的模型不确定性在真实系统中几乎必然存在。二是考虑量测数据的质量问题。我的实现假设量测都是干净的只是叠加高斯噪声。实际PMU数据里可能含异常值、通信丢包、时标偏差等。参考文献[3]里讨论了部分问题你可以在滤波流程中加入一个简单的数据质量判断模块先剔除坏数据再进入滤波更新很多情况下精度还能再提一截。三是滤波前端的可观测性分析。不要一上来就假设所有发电机都有PMU量测。实际工程中PMU覆盖是逐渐完善的可以先做可观测性分析看看哪些发电机没有直接量测时估计误差会被放大到不可接受的程度再决定PMU布点方案。这个方向能兼顾算法研究和工程实际。我在整个验证过程中最深的一点体会是动态状态估计的难点不在滤波器的数学公式本身而在于你需要把系统建模、量测配置、噪声参数整定、扰动场景设计这四件事串成一个自洽的整体。任何一个环节粗糙最终误差数据都会诚实地反馈给你。所以建议你第一遍完整跑通时不要贪图复杂场景先把稳态跟踪做到位再逐步加强扰动每加一步都要对比一次EKF和UKF的表现差异这样才能真正理解两种算法的脾气。