我最初接触电力系统动态状态估计DSE时第一反应是拿EKF和UKF直接各跑一遍Matlab代码比一比谁的误差更小。但真到动手写代码才发现关键瓶颈根本不是滤波器本身而是怎么把发电机动态模型和量测方程组织成可计算的形式。静态状态估计只需要解一个非线性最小二乘问题动态状态估计却要同时处理状态预测和量测修正两条链路模型稍微给错EKF和UKF都会在几十个采样点内飘飘飘地发散掉。这篇文章就围绕我用Matlab实现EKF和UKF做电力系统动态状态估计的完整过程展开。适合两类人一类是刚接触动态状态估计、想要一份能跑通的代码框架的初学者另一类是已经在用静态WLS加权最小二乘状态估计做断面分析想往预测型估计器转的工程师。你不需要提前懂很多但最好对卡尔曼滤波有一点点印象否则建议先把线性卡尔曼的五条公式过一遍再来读。1. 为什么静态状态估计满足不了动态分析的需求1.1 静态估计的冻结断面假设错在哪里传统电力系统状态估计是典型的静态问题给定某一时刻的全部量测求电压幅值和相角的最优解。SCADA系统大概每2到5秒推送一次遥测数据调度端拿到这批数据后做一次WLS估计输出一个断面。这个逻辑在运行方式变化缓慢的年代没有问题但新能源出力波动、负荷快速变化、直流输电功率调整这些场景出现后断面的间隔期里系统照样在演化量测到了下一帧却还没传上来调度员看到的状态就已经滞后了。更麻烦的是SCADA的刷新率并不均匀有的厂站遥测晚到几百毫秒甚至几秒状态估计器只能等全部数据齐了再算。而动态状态估计的思路完全相反用上一时刻的状态和系统模型去预测当前时刻的状态再用最新量测修正预测值量测不需要等齐来一个修正一次天然适合处理数据延迟和丢失。1.2 动态模型从哪里来做动态状态估计必须先有一个状态方程。电力系统里最经典的动态模型是发电机的摇摆方程它描述了转子转速与电磁功率之间的平衡关系。如果手头有发电机的详细参数状态量可以取功角、转速、暂态电动势如果只做母线级的估计文献里也有一种更务实的做法——把状态限定为节点电压幅值与相角状态转移用随机游走或一阶惯性方程近似。我在Matlab实现里采用的是后者状态向量直接取除参考母线外的所有节点电压幅值和相角。这样做有两个原因一是通用性强换一个系统拓扑不需要重新推导发电机模型二是量测方程直接复用潮流计算的映射关系省去了从发电机内电势到网侧状态量的复杂转换。代价是状态转移模型的物理意义弱一些Q矩阵需要认真调这部分留到后面的调参章节细讲。1.3 量测方程与状态方程的分工动态状态估计能成立靠的是状态方程和量测方程各管一段。状态方程负责描述系统的连续演化规律量测方程负责描述量测设备和状态之间的映射关系。在母线级模型里量测方程就是节点注入功率和线路潮流的非线性表达式和静态估计完全一样。这两条链路的配合很有意思状态方程预测不准量测修正还能往回拉量测数不够状态预测也能顶一阵子。冗余度低的可观测系统在静态估计里很难收敛动态估计反而能利用模型信息维持可观性这一点在实际工程里很有价值。2. EKF方案在非线性系统里做卡尔曼的代价2.1 从线性卡尔曼到扩展卡尔曼卡尔曼滤波本身是线性系统的产物五条公式里最核心的一步是误差协方差的传播P(k1|k)等于A乘以P乘以A转置再加Q。一旦状态方程或量测方程变成非线性这套传播公式就不能直接用。EKF的做法很简单也很暴力在工作点附近对非线性函数做一阶泰勒展开求出雅可比矩阵然后把雅可比当成线性卡尔曼里的系统矩阵A和量测矩阵H。电力系统的量测方程是潮流方程非线性程度中等偏强EKF的性能很大程度上取决于线性化工作点选得准不准。状态预测一步完成之后量测更新之前的预测状态就是线性化点预测误差越小雅可比越准确滤波效果越好。这正是EKF的软肋遇到强突变工况预测误差大线性化点偏离真实状态远雅可比矩阵失真然后在下一拍继续劣化。2.2 电力系统里的雅可比矩阵构造以母线注入有功为例极坐标下第i条母线的有功注入可以写成P_i V_i * Σ_j V_j (G_ij * cos(θ_i - θ_j) B_ij * sin(θ_i - θ_j))对状态向量里的每个电压幅值和相角求偏导才会得到量测矩阵H。我在代码里用符号工具箱先推导解析表达式再转成函数句柄。一开始偷懒试过数值差分求雅可比跑IEEE 39节点系统时误差不大但换到重负荷算例后协方差矩阵经常出现非正定后来全部改成解析雅可比才稳定下来。EKF的量测更新公式没有变化仍然是卡尔曼增益K、状态修正和协方差更新三条主线只不过H矩阵变成了当前工作点下的雅可比。需要注意的坑是雅可比矩阵每一列都要对应状态向量里的一个元素如果状态里混入了常数项或者参考母线的相角维度对不上Matlab里矩阵乘法的报错会排山倒海一样来。2.3 EKF的已知短板EKF最大的问题是一阶线性化在大扰动场景下不够用。振荡、甩负荷这些工况会让状态量跑出线性化区间雅可比矩阵的局部近似失效。还有一个工程上很痛的点雅可比矩阵的解析表达式和系统拓扑强相关改一次接线方式所有偏导公式都要连带修改维护成本极高。但我依然建议初学者先跑通EKF再上UKF。原因很简单EKF的每一行代码都能和卡尔曼五条公式对应上出错时容易排查。UKF虽然数学上更优雅调试起来反而更黑盒协方差出现问题不好定位。3. UKF方案绕开雅可比矩阵的确定性采样3.1 sigma点传播的核心思想UKF的思路和EKF完全不同它不去求导而是用一组精心挑选的sigma点来逼近非线性函数的统计特性。假设状态维度是nUKF会生成2n1个采样点这些点经过非线性函数传播之后用加权均值和加权协方差来近似真实的后验分布。在线性卡尔曼里协方差传播严格精确在EKF里是一阶近似在UKF里则是无迹变换的近似精度至少达到泰勒展开的二阶项对强非线性系统的适应性明显更好。sigma点的生成逻辑是在状态均值周围按照协方差矩阵的Cholesky分解结果做偏移偏移量由尺度参数决定。Matlab里直接调chol函数就能得到下三角矩阵但要注意如果P矩阵因为数值误差变得非正定chol会报错所以必须在每次协方差更新后加一个对称化和正定性修正。3.2 权重设计与参数选择sigma点的均值权重和协方差权重由几个参数控制通常取alpha1e-3beta2kappa0或3-n。alpha决定sigma点离均值的距离取太小可能让数值稳定性变差beta在高斯分布假设下取2是最优的kappa用于保证协方差矩阵的半正定性。这些参数不是金银细软随便调的alpha取多少直接关系到高阶项误差的放大程度我个人的经验是alpha落在1e-4到1e-2之间比较稳妥再小就要注意数值精度问题。UKF的不变性是一个容易被忽略的优点由于不需要求导换量测函数、换网络拓扑都只需要改函数本身滤波器的骨架完全不用动。对于EKF必须重新推导雅可比矩阵的场景UKF省掉的工程量非常大。3.3 为什么UKF在电力系统里更稳电力系统的量测函数是三角函数和乘积项的叠加在重负荷或者电压偏低的时候非线性程度接近临界EKF的一阶线性化会产生明显截断误差。UKF的sigma点传播等效于用多个采样点去拟合真实映射即便工作点移动精度也不会瞬间恶化。我在算例里观察到一个典型现象电压幅值从0.95pu往下掉的过程中EKF的误差开始周期性振荡UKF还能平稳跟踪。需要泼一点冷水UKF不是银弹。sigma点数量是2n1状态维度一旦上百每一拍的计算量会明显上涨。另外UKF对协方差的数值品质更敏感尤其在角度类状态上很小的舍入误差可能让权重出现负值最后导致协方差非正定。所以做UKF一定要在更新后强制检查P矩阵。4. Matlab代码实现状态推进、量测更新与整体架构4.1 数据组织方式我把整套代码按模块拆分最顶层是一个主脚本run_dse.m负责加载系统参数、生成量测、循环调用滤波器的predict和update方法、最后画图存结果。系统参数统一放在一个结构体里包括母线导纳矩阵Ybus、节点编号、量测位置索引这样换算例只改数据文件不需要动滤波逻辑。滤波器用Matlab的类来封装基类定义predict和update两个抽象接口EKF类和UKF类分别继承实现。这个设计让我后来加平方根UKF和CKF变得非常轻松只需要多写一个类文件主循环一行不需要改。对初学者来说面向对象不是必需品但如果你想反复换滤波器做对比研究这个设计能省很多重复劳动。4.2 量测数据生成动态状态估计仿真需要先有真值才能评价滤波器好坏。我先用潮流计算在若干时间断面上生成一组基准状态再把基准状态代入量测方程得到无噪声量测最后叠加上指定方差的高斯噪声。这里有一个关键技术细节状态真值的时间演化必须由状态方程生成而不是随便离散几条曲线。做法是先用状态方程从初始状态递推出真实状态序列再把每拍状态映射到量测空间。否则滤波器里的状态方程和仿真里的真实动态不一致评估结果毫无意义。量测类型我采用PMU风格的电压相量加上少量注入功率PMU采样率设定为每周期50帧即采样间隔0.02秒。这个时间尺度能体现动态估计的价值SCADA风格的秒级数据在这种算法演示里看不出太大区别。4.3 滤波主循环与EKF/UKF的分叉点每拍滤波循环中的逻辑可以浓缩成下面这段伪代码结构% 状态预测 x_pred f(x_est, u); if strcmp(filterType, ekf) [A, ~] stateJacobian(x_est, u); P_pred A * P_est * A Q; else P_pred sigmaPropagate(x_est, P_est, f); end % 量测更新 if strcmp(filterType, ekf) H measJacobian(x_pred); K P_pred * H / (H * P_pred * H R); else K ukfUpdate(x_pred, P_pred, z, h); end x_est x_pred K * (z - h(x_pred)); P_est P_pred - K * (H * P_pred); % EKF P_est (P_est P_est) / 2; % 强制对称EKF的分叉点主要在雅可比的计算和常规卡尔曼更新UKF的分叉点在sigma点的生成和传播。核心量测更新公式的形态一致这让比较两者时的变量控制很干净除了滤波器机理不同其他环节完全相同。4.4 误差评估与绘图我习惯同时保存每拍的估计误差、协方差迹以及卡尔曼增益的范数。单看误差曲线会误导人因为某一次随机噪声的波动可能很大一定要配合协方差迹观察滤波器是否自信得过火。绘图时用三张图电压幅值跟踪曲线、相角跟踪曲线、RMSE随时间变化曲线。相角曲线里要注意角度归一化不然从179度变到-179度的那一拍会被误判成巨大误差。5. 算例测试IEEE 39节点系统上的EKF与UKF对比5.1 测试场景设置我在IEEE 39节点系统上做了测试状态量取除参考母线外的38个相角加上39个电压幅值共77维。量测配置为一部分母线配置PMU电压相量再加若干注入有功和无功功率保证全网可观测。状态方程的动态模型采用一阶惯性加随机游走的混合形式时间常数按系统惯量大致整定。噪声协方差R按PMU精度设定电压幅值标准差0.001pu相角标准差0.001弧度。过程噪声协方差Q我花了不少时间调最后电压幅值对应元素取1e-6量级相角对应元素取1e-5量级。Q太大会让状态预测完全不信任模型滤波退化成逐拍静态估计Q太小又会让滤波器过于相信模型量测失灵时协方差收缩过快后面专门讲。5.2 实测曲线与误差行为以小扰动场景为例系统稳定运行到第200拍时某条母线负荷阶跃增加5%。两个滤波器都跟住了状态变化但细节差异明显。EKF在扰动发生后的第1至第3拍出现一个误差尖峰随后收敛回来UKF从第一拍起误差就控制在较小范围几乎看不到尖峰。这符合理论预期负荷阶跃瞬间状态轨迹快速弯曲一阶线性化的局部假设短暂失效。我还跑了一个更容易发散的场景把某台机组出力按正弦曲线大幅波动EKF产生了几次幅度较大的协方差收缩P矩阵一度接近奇异UKF则一直保持平滑。需要注意的是这两个都是单次仿真结果不同随机种子下的表现会有波动应该用Monte Carlo多跑几十遍再下结论。5.3 数值指标与时间开销单次仿真结束以后我统计了RMSE和单拍平均耗时。典型结果EKF的电压幅值RMSE在0.0031pu左右UKF在0.0022pu左右相角RMSE两个滤波器差距更明显EKF约0.0028radUKF约0.0017rad。时间上EKF的单拍耗时要低一些因为一次量测更新只需要计算一个雅可比矩阵UKF要传播155个sigma点。但这不代表EKF永远更快解析雅可比推导需要人工成本如果算上推导时间UKF反而更划算。我在项目里已经把解析雅可比提前算好存成函数文件所以单拍耗时才压得比较低。5.4 什么时候EKF够用如果系统运行平稳、量测冗余度高、没有频繁的大扰动EKF的精度和UKF差距并不大有时甚至因为数值更稳定反而表现更好。UKF的优势在强非线性或扰动频繁的场景下才能体现出来。所以我不建议盲目追求UKF先把系统的动态特性想清楚再选。6. 调参心得与滤波器发散避坑6.1 Q矩阵和R矩阵怎么定R矩阵相对好办按量测设备的精度指标来就可以。Q矩阵是动态估计里最让人头疼的参数。它没有明确的物理标定方法只能从模型误差的估计出发你的状态方程离真实动态有多远Q就设多大。模型越粗糙Q应该越大给量测修正留出余地。我在代码里写了一个辅助函数可以在一定范围内扫描Q的乘子系数自动输出不同系数下的估计RMSE。这样调参不是盲目的而是先看RMSE谷值在哪一段再在附近做细化。对于77维状态全矩阵扫描不可能我把它拆成幅值子块和相角子块分别给出统一的乘子。经验上电压幅值子块和相角子块的Q系数可以差一个数量级混在一起调很难收敛。提示如果滤波曲线出现量测修正过头的表现——估计值高频抖动、误差方差忽大忽小优先减小Q如果出现量测拉不动的表现——真值突变后估计值缓慢爬行优先增大Q。这两个症状是调Q的可靠风向标。6.2 协方差矩阵的非正定问题P矩阵非正定在EKF和UKF里都出现过是动态估计调试期最常遇到的坑。症状是Matlab报错说chol输入矩阵必须正定或者RMSE曲线突然爆表。我做了三件加固第一每次P更新后强制对称化也就是(PP)/2第二加一个小的对角修正项类似于Levenberg-Marquardt的做法让最小特征值保持在某个阈值之上第三定期检查P的特征值如果出现显著负值说明模型发散了靠数值修补救不回来必须回头调Q。UKF对P矩阵的正定性要求比EKF更高因为sigma点生成依赖Cholesky分解。我最开始在39节点系统上跑UKF发散十次里有七八次都是这里出的问题。修正之后稳定了很多但也要记住UKF发散的路径往往比EKF更剧烈调参时要边跑边看协方差迹曲线。6.3 采样周期与离散化的陷阱动态状态估计的采样周期直接决定状态转移模型的形式。PMU的50帧每秒采样对应20毫秒间隔在这个尺度上很多连续动态可以用简单的一阶离散公式近似但如果数据源是SCADA的秒级刷新状态方程里的时间常数就必须重新整定否则一拍之内状态变化太大预测根本不准。还有角度处理的细节相角是周期量预测值和量测值之间做差一定要归一到[-pi, pi]区间否则遇到180度附近的跳变会产生假的大误差进而拉歪卡尔曼增益。这个坑在Matlab里尤其隐蔽因为Matlab的sin和cos函数不会提醒你角度差异异常。6.4 我做这套代码后的几条总结滤波发散时先别动滤波器结构先把Q、R两个矩阵和初始协方差P0过一遍大部分问题都出在这三个参数上。P0设得太小会让滤波器一开始就过度自信量测稍微不稳就发散P0设太大又会让前几十拍的估计噪声偏大。我一般取对角元素0.1到1之间的对角阵让滤波器先通过前几十拍把协方差降下来。如果你也是从静态估计切到动态估计的我的建议是先别急着上UKF。老老实实把EKF在简单系统上跑通体会状态预测和量测修正各占多少分量再把系统换到39节点最后才换成UKF做对比。这样每一步出了偏差你都清楚该去查哪一段代码。这套Matlab框架本身可以直接拿去做二次开发换数据文件、换量测配置、换滤波器类型都不需要重写主循环。我后来加扩展卡尔曼的强跟踪版本也只用了一天时间架构上省下来的功夫相当可观。