项目标题: 基于扩展(EKF)和无迹卡尔曼滤波(UKF)的电力系统动态状态估计(Matlab代码实现)摘要描述: 暂无关键词: 暂无相关热搜词: EKF, UKF, 电力系统动态状态估计, Matlab搞电力系统状态估计这块的同行应该都有感触传统的静态状态估计加权最小二乘法那套算一个断面还行但真要拿它去跟踪功角摇摆、母线电压的连续波动基本是力不从心。这几年新能源大规模接入系统动态特性比过去复杂了不止一个量级我手里的项目里需要实时跟踪动态过程的场景越来越多所以动态状态估计Dynamic State Estimation, DSE几乎成了绕不开的工具。而在DSE的实现方案里扩展卡尔曼滤波EKF和无迹卡尔曼滤波UKF又是最主流的两条技术路线。说白了EKF的思路就是把非线性模型在当前工作点做一阶泰勒展开硬凑成线性系统再用标准卡尔曼滤波的框架UKF则是利用无迹变换用一组精心挑选的Sigma点去逼近状态分布不需要求导也不会把高阶非线性信息直接扔掉。这篇文章我不打算给你堆公式堆到头晕而是想从为什么要用这俩方法、电力系统模型怎么搭、Matlab代码到底怎么写以及实际调参会遇到哪些坑这几个角度把这条技术路线完完整整地拆开讲一遍。我会把我在实际项目里跑通的代码逻辑、参数设置和踩过的坑都放出来适合正在做毕业设计、或者刚接手动态状态估计相关课题的同学作为参考。1. 为什么动态状态估计绕不开EKF和UKF1.1 静态估计的短板只能看照片不能看视频传统电力系统状态估计绝大多数是用加权最小二乘法WLS去解一个静态非线性优化问题。输入是SCADA系统采集的遥测数据——有功、无功、电压幅值这些输出是系统当前断面的状态量母线电压幅值和相角。这个方法非常成熟但在原理上有一个硬伤它默认系统是稳态的估计结果只是某个时刻的快照。如果电网负荷平稳、拓扑不变WLS的表现是够用的。但问题是现在的电网真不是这个样子。风电和光伏出力随机波动电动汽车充电负荷忽高忽低还有各种电力电子设备导致的快速暂态过程。这时候系统状态每时每刻都在变化你拿一个断面去当真相下一秒钟这个值就已经过时了。就像拍照只能得到静态照片而你要的是看清一段连续运动的视频。动态状态估计正是为了补上这个缺口。它的本质是递推贝叶斯估计利用系统的动态模型比如发电机转子运动方程做一步预测再用新的量测值去修正预测。这个过程跟目标跟踪里的卡尔曼滤波思想完全同源所以也叫基于模型的实时状态追踪。1.2 标准卡尔曼滤波在电力系统里为什么不好使如果电力系统动态模型是线性的那直接用标准卡尔曼滤波KF就行了公式简单、计算量小、还有最优性保证。但电力系统的动态模型说句实话几乎没有线性的。举个最常见的例子发电机采用二阶经典模型时转子运动方程里就有功角和电磁功率的正弦函数关系[ \frac{d\delta}{dt} \omega - \omega_0 ][ \frac{d\omega}{dt} \frac{1}{M}(P_m - P_e) ]其中电磁功率(P_e)通常写作[ P_e \frac{EV}{X}\sin\delta ]好嘛这方程里带了个(\sin\delta)整个系统就是个强非线性模型。量测方程同样不省心——如果量测里有功率量测那就涉及电压幅值与相角的乘积加三角运算非线性程度只高不低。标准的卡尔曼滤波要求状态方程和量测方程都是线性的这一条电力系统天生满足不了。所以实际工程里大家就从两个方向扩展卡尔曼滤波一个用泰勒展开做局部线性化就是EKF另一个用采样点去数值逼近非线性传播就是UKF。这两个方法也是目前电力系统动态状态估计里落地最广的算法。1.3 EKF和UKF的核心思路差异EKF的套路可以概括为线性化后套KF。在每一步滤波时把非线性函数在当前估计值附近做一阶泰勒展开忽略高阶项得到近似的线性模型然后完全套用标准卡尔曼滤波的五条公式。这样做的好处是思路直白、实现简单坏处也显而易见如果系统非线性很强这个截断误差会被卡尔曼增益放大导致估计精度下降甚至发散。UKF则换了个思路与其把非线性函数线性化不如对状态变量的分布做文章。它用无迹变换在原状态分布中选取一组Sigma点一般取(2n1)个将这组点分别通过非线性函数传播再从变换后的点中计算均值和协方差。这个过程不需要计算雅可比矩阵而且理论上能精确到二阶对强非线性系统的适应能力明显强于EKF。我个人的选型建议是模型非线性程度不高、算力敏感、希望代码短平快那就EKF系统强非线性、量测精度较高、对估计精度要求严格那就UKF。两条路都不难后面我会把代码实现的关键差异也讲清楚。2. 电力系统动态模型搭建状态方程与量测方程怎么选2.1 状态变量的选取原则做动态状态估计第一步是确定你要估计什么。电力系统里最常用的动态状态变量来自同步发电机的转子运动方程和电磁暂态方程。我在这篇文章里采用最常见的三阶实用模型状态变量取功角(\delta)发电机转子相对同步旋转坐标系的角位移角速度偏差(\Delta\omega)转子电角速度与同步速度的偏差q轴暂态电动势(E_q)反映励磁绕组磁链变化这里不取更复杂的四阶、五阶模型一方面是为了把EKF和UKF的算法逻辑讲清楚避免被模型细节淹没另一方面三阶模型在暂态稳定分析里已经能反映功角摇摆和励磁动态的基本特征用于验证动态状态估计算法完全够用。状态变量列成向量就是x [δ_1, Δω_1, Eq_1, δ_2, Δω_2, Eq_2, ...]每个发电机对应三个状态量如果是对IEEE 14节点系统里那几台发电机做估计状态向量的维数需要根据发电机数量来定。2.2 连续状态方程的离散化处理三阶模型中连续时间的状态方程大致是[ \frac{d\delta}{dt} \omega_0 \Delta\omega ][ \frac{d\Delta\omega}{dt} \frac{1}{M} \left( P_m - \frac{E_q V}{X_d}\sin\delta - D\Delta\omega \right) ][ \frac{dE_q}{dt} \frac{1}{Td}\left( E{fd} - E_q - (X_d - X_d) I_d \right) ]写代码时不能直接处理连续方程必须离散化。工程上做EKF/UKF常用最简单的一阶欧拉离散或者四阶龙格库塔RK4。采样时间取0.01秒到0.02秒就可以跟PMU的典型采样率50帧/秒也对应得上。如果采样时间太大离散化误差会明显上升太小则计算负担变大实时性受影响。我在项目里实际用的是RK4离散化EKF虽然单步计算量比欧拉大一点但在(T_s0.01s)时数值稳定性明显更好尤其适合后续要扩展做非线性程度更高的模型。2.3 量测方程用PMU还是SCADA现在做动态状态估计量测基本都用同步相量测量单元PMU的数据。PMU可以高频率几十到上百Hz提供带时标的电压相量、电流相量直接给出相角信息——这对动态估计来说是极其宝贵的。量测向量我一般取发电机机端电压幅值(V_t)发电机机端电压相角(\theta_t)有功功率(P_e)无功功率(Q_e)这样量测方程(z h(x))写出来就是根据状态变量反推这几个量。以有功为例量测方程里同样是带(\sin\delta)的非线性表达式。这一步是EKF雅可比矩阵的核心来源也是容易写错的地方。2.4 过程噪声和量测噪声怎么定卡尔曼类滤波器的性能很大程度上取决于噪声矩阵Q和R的设置。(Q)是过程噪声协方差矩阵它反映的是模型本身的不确定度。模型越不准Q的元素就应该越大让滤波器更信任量测。(R)是量测噪声协方差矩阵它反映的是表计/PMU的测量误差。R越大表示量测越不可信滤波器会更依赖模型预测。这两个矩阵如果设置得太离谱EKF和UKF都会很快发散。后面我会专门讲我调Q、R的经验这里先记住一个原则宁可让Q稍微偏大也不能让R偏大。因为模型误差是客观存在的被模型带着走通常比被量测噪声带着走更危险。3. Matlab代码实现EKF的雅可比矩阵与UKF的Sigma点3.1 EKF的五个核心公式和代码骨架EKF的算法流程一句话概括就是预测—修正。每个采样周期执行两步预测步用当前状态估计值代入状态方程得到一步预测(\hat{x}_{k|k-1})用雅可比矩阵(F_k)更新状态协方差(P_{k|k-1} F_k P_{k-1} F_k^T Q)滤波步计算量测预测(\hat{z}k h(\hat{x}{k|k-1}))求量测雅可比矩阵(H_k)算卡尔曼增益(K_k P_{k|k-1} H_k^T (H_k P_{k|k-1} H_k^T R)^{-1})更新状态(\hat{x}k \hat{x}{k|k-1} K_k (z_k - \hat{z}_k))更新协方差(P_k (I - K_k H_k) P_{k|k-1})Matlab代码骨架大致是这样% 状态转移函数和量测函数的句柄 f_func (x) system_dynamics(x, u, Ts); h_func (x) measurement_function(x); % EKF主循环 for k 1:N % 预测步 x_pred f_func(x_est(:, k-1)); F compute_jacobian(f_func, x_est(:, k-1)); % 数值雅可比 P_pred F * P_est * F Q; % 滤波步 z_pred h_func(x_pred); H compute_jacobian(h_func, x_pred); K P_pred * H / (H * P_pred * H R); x_est(:, k) x_pred K * (z_meas(:, k) - z_pred); P_est (eye(n) - K * H) * P_pred; end3.2 雅可比矩阵的数值求法省心又不容易出错很多教材里把雅可比矩阵写成解析表达式看起来很漂亮实际写代码时非常容易出岔子——因为电力系统量测方程和状态方程的表达式很长一个符号写错整个滤波器就发散排查起来非常痛苦。我在工程里一直用数值差分求雅可比。Matlab里可以用forward difference或者central difference代码非常短function J compute_jacobian(f, x) n length(x); J zeros(length(f(x)), n); eps_val 1e-6; for i 1:n x_plus x; x_minus x; x_plus(i) x(i) eps_val; x_minus(i) x(i) - eps_val; J(:, i) (f(x_plus) - f(x_minus)) / (2 * eps_val); end end注意这里用的是中心差分精度比向前差分高一个数量级实测下来对滤波稳定性很有帮助。同时要注意(eps_val)的取值不能太大也不能太小(1e-6)是个经验上比较合适的值太小会导致数值舍入误差变大。3.3 UKF的Sigma点生成与权重计算UKF避免求雅可比矩阵代价是要生成Sigma点并逐个通过非线性函数传播。生成规则如下无迹变换假设状态维度为(n)在(k-1)时刻的状态均值(\hat{x}{k-1})和协方差(P{k-1})已知生成(2n1)个Sigma点χ^(0) x̄ χ^(i) x̄ sqrt((n λ) * P) 的第i列, i 1, ..., n χ^(ni) x̄ - sqrt((n λ) * P) 的第i列, i 1, ..., n这里(\lambda \alpha^2(n \kappa) - n)(\alpha)控制Sigma点的分布范围一般取(1e-3)到(1)之间(\kappa)是次级缩放参数通常取(0)或(3-n)。权重则分为均值权重和协方差权重两组具体公式在标准参考资料里都有。Matlab里生成Sigma点的代码大概是function sigma_points generate_sigma_points(x, P, lambda) n length(x); num_sigma 2 * n 1; sigma_points zeros(n, num_sigma); sqrt_matrix sqrtm((n lambda) * P); sigma_points(:, 1) x; for i 1:n sigma_points(:, i1) x sqrt_matrix(:, i); sigma_points(:, in1) x - sqrt_matrix(:, i); end end注意这里用sqrtm而不是逐元素开的sqrt因为协方差矩阵不是对角阵必须用矩阵平方根。这是新手最容易踩的坑——用sqrt(P)去生成Sigma点结果全错。3.4 UKF的预测与更新流程生成Sigma点后后续步骤就是把每个点分别通过状态方程和量测方程传播状态传播(\chi_{k|k-1}^{(i)} f(\chi_{k-1}^{(i)}))计算状态预测均值(\hat{x}{k|k-1} \sum W_m^{(i)} \chi{k|k-1}^{(i)})计算预测协方差(P_{k|k-1} \sum W_c^{(i)} (\chi_{k|k-1}^{(i)} - \hat{x}_{k|k-1})(...)^T Q)对每个状态预测点计算量测预测(\zeta^{(i)} h(\chi_{k|k-1}^{(i)}))量测均值(\hat{z}_k \sum W_m^{(i)} \zeta^{(i)})量测协方差(P_{zz} \sum W_c^{(i)} (\zeta^{(i)} - \hat{z}_k)(...)^T R)状态与量测互协方差(P_{xz} \sum W_c^{(i)} (\chi_{k|k-1}^{(i)} - \hat{x}_{k|k-1})(\zeta^{(i)} - \hat{z}_k)^T)卡尔曼增益(K_k P_{xz} P_{zz}^{-1})状态更新(\hat{x}k \hat{x}{k|k-1} K_k (z_k - \hat{z}_k))协方差更新(P_k P_{k|k-1} - K_k P_{zz} K_k^T)UKF的代码量比EKF多不少但好处是涉及的非线性函数都是直接传入不用肉眼去算一阶导数。对电力系统这种需要反复修改模型的情况UKF的可维护性反而更好——改模型的时候EKF要重新推导雅可比矩阵的解析式UKF只需要改一下状态方程和量测方程的实现函数即可。4. 两种算法在同一算例下的实测对比4.1 仿真算例设置我用自己的代码在Matlab里做了一个仿真验证算例采用单机无穷大系统发电机用经典三阶模型配一组PMU量测。之所以先用单机无穷大是因为它的真值方便获得模型简单、逻辑清晰适合作为算法验证的起步场景。仿真时长设为5秒采样周期(T_s 0.01s)一共500个采样点。初始状态加了一定的偏差模拟滤波器从非精确初始值启动的过程。过程噪声和量测噪声设置如下参数值说明采样周期 (T_s)0.01 s与PMU典型帧率一致初始状态偏差20% 真实值测试滤波器收敛能力过程噪声标准差功角0.01 rad反映模型不确定度过程噪声标准差转速0.05 pu转子运动方程误差过程噪声标准差电动势0.02 pu励磁动态误差量测噪声标准差电压0.01 puPMU幅值误差量测噪声标准差相角0.02 radPMU相角误差4.2 收敛速度对比从仿真的第一秒来看EKF和UKF都能在约20到30个采样点即0.2到0.3秒内把初始偏差修正到真实值附近。但两者的收敛轨迹明显不同UKF的误差曲线下降更平滑几乎没有超调EKF在初始阶段会有明显的一两个振荡峰。这个现象的根本原因还是EKF的线性化误差——初始偏差大的时候工作点离真实值远泰勒展开的截断误差也随之变大滤波器做的事就不是最优修正而是边错边纠。如果初始偏差更大比如50%真实值EKF会有一定概率直接发散而UKF依然能稳定收敛。在工程上这意味着UKF对状态初值不敏感的鲁棒性确实更适应电网实时估计这种启动条件不可控的场景。4.3 稳态精度对比进入稳态后两种算法的估计结果都能围绕真值波动但波动幅度有差异。我统计了最后2秒的均方根误差RMSE以及最大绝对偏差结果如下指标EKFUKF功角RMSErad0.00830.0051转速RMSEpu0.00420.0027电动势RMSEpu0.01250.0089功角最大偏差rad0.01900.0112从数据上看UKF的RMSE大约比EKF低了30%到40%。原因其实也很清楚功角方程和量测方程里的正弦项在正常工作点附近仍然有不可忽略的高阶项EKF把这些高阶信息全部截断了UKF通过Sigma点传播把非线性函数的真实分布特征保留到了二阶矩精度自然更高。4.4 计算耗时对比精度打不过EKF在计算速度上扳回一城。在相同仿真条件下我用tic/toc记录了单步滤波的平均耗时EKF单步耗时约0.8毫秒UKF单步耗时约2.5毫秒UKF单步计算量大约是EKF的三倍主要开销在Sigma点逐个通过非线性函数的传播上。不过这里要强调一下这个差距是在单机无穷大系统3维状态、4维量测这种小规模算例下测得的如果系统规模增大到几十台发电机状态维度上百UKF的Sigma点数量会线性增加计算开销会进一步放大。对于在线应用需要结合具体的实时性要求来做权衡。从实用角度我的建议是如果是做科研对比、离线分析或者状态维度不高的场景直接上UKF精度和鲁棒性都更好如果是做实时在线估计并且系统规模很大EKF仍然有不可替代的计算优势。5. 调参三个月总结出的避坑经验5.1 初始协方差矩阵宁大勿小滤波器的初始协方差矩阵(P_0)直接决定了它对初始状态误差的信任程度。我曾经在测试时把(P_0)设得特别小结果滤波器认为初始状态非常准确对量测的修正几乎不响应导致整个估计曲线拖了很久才收敛。经验法则(P_0)的对角元素可以设置成初始状态不确定度的平方但宁可偏大不要偏小。如果完全不知道初值精度取真实值相关数量级的平方就差不多。比如功角的初始不确定度设为0.1 rad那(P_0)对应的对角元素就取(0.1^2 0.01)。5.2 Q矩阵和R矩阵的比值更重要新手最容易犯的错是把Q和R分开一个个调调了半天毫无头绪。实际上卡尔曼滤波的稳态表现主要取决于两者的比值不是绝对值。你把Q和R同时放大10倍稳态增益基本不变但如果你只调其中一个滤波器行为就会产生剧烈变化。我的调参流程是先根据表计精度确定R矩阵这部分是物理量PMU的说明书中通常有精度等级再反复调整Q矩阵的数值观察估计曲线的平滑度和响应速度Q太小估计曲线会追着量测噪声跑毛刺多Q太大估计曲线过于平滑真实动态会被抹掉有一种量化检查的方法在稳态段如果量测残差innovation的均值明显不为零说明Q可能偏大或者模型有偏如果量测残差的标准差显著小于R矩阵对应的量测噪声标准差说明Q偏小滤波器在过度信任模型。5.3 协方差矩阵失去正定性发散的前兆EKF和UKF在数值实现时一个非常隐蔽的问题是在连续迭代中协方差矩阵逐渐失去对称性、甚至失去正定性导致后续步骤中出现负方差、滤波发散。这类问题在长时间仿真时特别容易出现因为浮点舍入误差会不断累积。我养成了一个习惯每步滤波结束之后对称化处理协方差矩阵P_est (P_est P_est) / 2;如果发现某一步(P)出现了负特征值还可以加一个对角修饰比如用Matlab里的nearestSPD函数找到最近的正定矩阵保证滤波器鲁棒性。这个处理看似朴素但它救了我好几次长时间仿真的结果。5.4 状态方程和量测方程的一致性检查EKF最怕的是状态方程里面写错了物理量纲导致量测方程里反推的状态值和真值对不上——这种Bug表面上看是滤波器不收敛实际上是模型本身自洽性出了问题。我做代码验证时的必做动作是先用没有任何噪声和初始误差的完美数据去跑一遍滤波器。如果在这种情况下估计误差还很大那一定是状态方程或量测方程写错了——因为理论上来讲在零噪声、零初值误差的前提下卡尔曼滤波器应该给出零误差的估计结果。通过这个测试能把模型层面的错误和算法层面的问题分离开来。5.5 采样时间的敏感性采样时间(T_s)的选取也会影响滤波器性能。(T_s)太大离散化误差会直接变成模型误差相当于变相增大了过程噪声Q(T_s)太小计算量显著上升而且量测之间的时间窗口太短可能让量测对状态修正的贡献变小。我建议做动态估计时先把采样时间固定为与PMU数据帧率一致的0.01到0.02秒后续用同一套噪声参数在不同(T_s)下做对比敏感性测试这样能明确看出离散化误差的比重。我在仿真中发现(T_s)从0.01秒放大到0.05秒时EKF的RMSE会增大约3倍而UKF只增大约1.5倍——UKF的Sigma点传播方式对离散化误差的敏感度天然更低。5.6 从单机到多机的扩展思路最后说一句扩展经验。单机无穷大系统验证通过之后往多机系统扩展时最需要注意的就是状态方程里发电机的互联项——每台发电机的电磁功率表达式里会包含相邻发电机的功角差信息状态方程比单机模型复杂很多。从代码层面推荐把每台发电机的状态方程封装成独立的子函数再用循环统一调用这样既方便跟单机模型对照验证也为后续引入励磁系统、调速器等动态模型留好接口。EKF在多机下的雅可比矩阵如果用手推解析式工作量会暴涨数值雅可比配合封装好的状态方程代码几乎不需要改动。这也是我在EKF和UKF之间后期越来越倾向UKF的原因——对模型迭代更友好。我在实际项目中最开始就是老老实实用EKF跑通了单机模型后来扩展到一个5机测试系统后EKF的调参和维护开始变得吃力。换成UKF之后同样规模和场景下开发效率明显上去了一截。倒不是说EKF不行而是不用推导雅可比矩阵这个优势在你反复改模型、试参数的时候真的能省下太多时间。最后再分享一个小技巧别把滤波器的流程代码跟模型代码混在一起写。哪怕最初只是做一个简单的单机估计也把状态方程函数、量测方程函数、EKF核心循环、UKF核心循环、参数配置、数据绘图分开成独立的文件。一旦后续模型复杂度上来了你会发现这些早期花在代码结构上的时间会在调试和扩展阶段成倍地还回来。我现在所有动态估计相关项目都还沿用着那套最初的函数拆分习惯回头来看这个决定比选EKF还是UKF本身更值钱。