电力系统动态状态估计：EKF与UKF算法对比及Matlab实现

电力系统动态状态估计这几年因为PMU大规模铺开逐渐从实验室走向工程应用。用卡尔曼滤波的思路去跟踪发电机的功角和转速是其中最常见的技术路线。EKF和UKF作为两种主流非线性滤波算法到底怎么选、怎么在Matlab里落地我把自己跑过的代码和踩过的坑整理成这篇文章给做动态状态估计的同行一个直接可参考的方案。这篇文章不是教科书式原理堆砌而是以实操为导向适合刚接触动态状态估计的研究生也适合要在Matlab里快速验证算法、对比EKF和UKF性能的工程师。你会得到的核心内容包括两种算法的适用边界、单机无穷大系统下的建模与离散化、可直接改用的Matlab代码骨架以及我在调参时遇到的大量现场问题和排查思路。读完你可以直接动手跑通一个完整的动态状态估计流程。1. 为什么电力系统动态状态估计需要卡尔曼滤波1.1 动态状态估计解决什么实际问题传统的电力系统状态估计尤其是基于SCADA的量测本质上是一个静态估计问题。数据采样周期是秒级甚至分钟级默认系统运行在一个相对平稳的工况然后用加权最小二乘去解一个非线性代数方程组。这个思路在稳态调度、EMS分析里没有问题可一旦系统出现负荷快速波动、低频振荡、暂态过程SCADA的量测时间分辨率根本不够用静态估计解出来的结果只是一个模糊的“平均状态”无法反映发电机转子运动和电压随时间变化的真实轨迹。动态状态估计就是要解决这个时间维度上的缺口。它的输入来自两方面一个是系统动态模型比如发电机的转子运动方程告诉我们状态是怎么随时间变化的另一个是高时间分辨率的量测比如PMU测得的相量、频率、功率刷新率可达每秒几十帧甚至上百帧。卡尔曼滤波这类递推算法天然适合把模型预测和在线量测融合起来每一帧都给出状态的后验估计并附带不确定性信息也就是协方差矩阵。这个特性让动态状态估计在动态安全评估、振荡监测、模型参数校核等场景里有直接价值。实际工程中我们关心的状态量通常是发电机的功角δ和转速ω可能再加上暂态电势等。量测则来自PMU可能是电压相角、频率也可能通过功率变换得到。这里的关键点在于状态方程和量测方程都不是简单的线性关系特别是功率与功角之间存在正弦函数这就引出了非线性滤波问题。1.2 从线性卡尔曼到非线性滤波标准卡尔曼滤波是最小方差意义下的最优线性滤波器但它有两个硬性前提状态方程和量测方程都是线性的噪声都服从高斯分布。如果把发电机摆动方程离散化状态更新里会出现正弦项如果把功率作为量测量测方程又会是电压、相角的非线性函数。强行套用标准卡尔曼滤波会失效因为协方差传播时无法保留非线性函数的真实分布特性。针对这个问题工程上有两条主流近似路线。一条是EKF把非线性函数在当前估计点做一阶泰勒展开用雅可比矩阵替代原来的线性转移矩阵另一条是UKF通过一组精心选择的Sigma点去传播非线性函数用采样点的方式近似状态的后验均值和协方差。EKF思路上更接近“局部线性化”UKF则是“无需算导数直接传点”。两种方法我都在Matlab里实现过各有取舍下面我会从原理和代码两个层面拆开细说。2. EKF与UKF的核心原理对比2.1 EKF线性化近似EKF的核心思想不难理解既然系统是非线性的那就把非线性函数在当前估计值附近展开成线性函数剩下的照搬标准卡尔曼滤波流程。每一步递推中预测阶段用系统方程直接算先验状态同时用状态方程的雅可比矩阵F去传播协方差更新阶段用量测方程算预测量测用雅可比矩阵H去计算卡尔曼增益。这里最容易被忽略的是雅可比矩阵的“时效性”。EKF的线性化是在每一步的标称点上做的标称点本身随估计值移动所以每次都需要重新计算F和H。如果系统强非线性标称点处的线性近似误差会很大甚至让滤波发散。我在多机系统上试过当扰动幅度增大到一定程度EKF的估计误差会出现明显突变这时候不是代码写错了而是线性化假设撑不住了。计算量方面EKF在低维状态空间里确实便宜。雅可比矩阵如果解析求导麻烦可以用数值差分代替Matlab里也支持符号求导再转成函数句柄。我在实际项目里更推荐先用符号工具箱把雅可比推导出来生成matlabFunction然后在滤波循环里调用这样既避免了手算公式出错又比数值差分稳定。2.2 UKF无迹变换UKF走的是完全不同的路线。它不追求线性化而是用无迹变换也就是选取一组权重确定的Sigma点让这些点在均值附近以特定方式分布经过非线性函数传播后用它们的统计特性来近似真实的后验分布。对我来说UKF最直观的好处就是不需要推导雅可比矩阵对强非线性问题有更强的鲁棒性。具体来说对于一个n维状态向量UKF生成2n1个Sigma点。这些点按照状态的均值和协方差展开其中中心点就在均值上其余的沿协方差矩阵的主轴方向按比例偏移。偏移量由参数alpha、beta、kappa控制。每个Sigma点都有对应权重当它们全部通过非线性函数之后加权平均和加权协方差就重构出了状态的后验统计数据。我在电力系统动态状态估计里用UKF之后最明显的感受是当量测方程里三角函数带来的非线性很强时UKF估计出的功角曲线更平稳没有EKF那种毛刺感。代价是单个滤波步里需要多次调用系统模型函数相当于2n1次函数评估计算量比EKF高出数倍。不过对于单机或中小规模系统这个代价完全可接受很多研究者也习惯在实验室先跑通UKF再考虑切换到更复杂的Sigma点改进算法。2.3 精度与计算代价的取舍选择EKF还是UKF不能只盯着滤波精度要连同系统模型复杂度和实时性约束一起判断。我列一个自己的经验对照方便你做初步选型对比维度EKFUKF实现原理一阶泰勒展开需要雅可比矩阵无迹变换需要Sigma点生成非线性处理能力弱非线性下足够强非线性易发散强非线性下仍能保持较好的鲁棒性计算量较小适合机群规模较大时偏大约是EKF的3倍左右调试难度雅可比矩阵写错很隐蔽Sigma点参数需要调节但调试相对直观适用场景弱非线性、要求实时性高、模型相对可信强非线性、模型未知部分多、可以采用偏重精度这个表是我在实践中最常用来向别人解释取舍的框架。你如果是做在线实时估计状态维数又高EKF已经是工程上稳妥的起步选择如果项目允许更大的计算开销或者你发现EKF始终无法稳定收敛那就果断换UKF。3. Matlab代码实现从建模到滤波3.1 电力系统动态模型与量测模型我在验证算法时最常用的是单机无穷大系统。这个模型虽然简单但它包含了发电机转子动态、负荷变化和线路传输的非线性关系足够用来对比EKF和UKF的算法行为。连续时间下的发电机摆动方程可以写成delta_dot omega - omega_s omega_dot (Pm - Pe) / M - D * (omega - omega_s) / M其中delta是发电机功角omega是电角速度omega_s是同步角速度Pm是机械功率Pe是电磁功率M是惯性时间常数D是阻尼系数。电磁功率Pe在经典二阶模型中通常表达为Pe (E * V / X) * sin(delta)E是暂态电势V是无穷大母线电压X是传输电抗。这样状态方程本身就是一个非线性常微分方程组。离散化时我倾向于用欧拉法做粗步长验证但如果你想更接近真实动态建议用二阶龙格库塔代码复杂度也不高。量测方面最方便的做法是假设PMU可以直接给出功角和转速的带噪声测量这时量测方程是线性的。但为了把EKF和UKF的非线性特性充分体现出来我更习惯用电磁功率作为量测因为Pe与delta之间是正弦关系非线性和缓但不平凡。加上实际PMU也能给出功率量测这种设置更贴近真实工程。整个动态状态估计的过程就是给定系统模型函数f(x)和量测函数h(x)初始化状态与协方差然后在每一时刻递推执行预测和更新两个步骤。3.2 EKF实现步骤与代码骨架EKF在Matlab里实现起来的代码结构非常清晰。我先把核心循环骨架贴出来注意这个版本为了可读性把雅可比矩阵用数值差分处理方便跑通后再换成解析式。% 系统参数 M 10; D 2; Pm 0.8; E 1.2; V 1.0; X 0.3; omega_s 1.0; dt 0.005; % 状态: x1 delta, x2 omega f (x) [x(2) - omega_s; (Pm - (E*V/X)*sin(x(1)) - D*(x(2) - omega_s)) / M]; h (x) (E*V/X)*sin(x(1)); % 量测: 电磁功率 % 数值雅可比 Jf (x) numericalJacobian(f, x, dt); Jh (x) numericalJacobian(h, x, dt); % 初始化 x_hat [0.3; 1.0]; P diag([0.1, 0.01]); Q diag([1e-4, 1e-4]); R 1e-3; % 滤波循环 for k 1:N % 预测 x_pred x_hat dt * f(x_hat); F Jf(x_hat); P_pred F * P * F Q; % 更新 z_pred h(x_pred); H Jh(x_pred); K P_pred * H / (H * P_pred * H R); x_hat x_pred K * (z_meas(k) - z_pred); P (eye(2) - K * H) * P_pred; end有几个细节我想特别提醒。第一上面的f是连续时间导数预测时我用了一个简单的欧拉离散x_pred x_hat dt * f(x_hat)。如果你把Jacobian也按连续时间直接算再离散容易把时间步长因子弄丢导致协方差传播不匹配。我习惯是先离散状态方程再对离散方程求雅可比。第二量测方程h只返回一个标量功率H在这时是1x2的雅可比向量代码里用数值差分会自然得到正确的维度。如果要扩展到多个量测只需把h改成列向量输出H变成mxn矩阵即可。第三数值雅可比函数需要小心选择扰动步长。我通常用中心差分步长取sqrt(eps)相对量级例如sqrt(eps)*max(1,abs(x))比一阶前向差分稳得多。工程调试中数值雅可比造成的发散很常见解析式能避免这个问题。3.3 UKF实现步骤与代码骨架UKF代码比EKF多一个Sigma点生成环节。我给出一个标准的Matlab实现骨架关键步骤都做了注释。n 2; % 状态维数 alpha 1e-3; kappa 0; beta 2; lambda alpha^2 * (n kappa) - n; % 权重 Wm zeros(2*n1,1); Wc zeros(2*n1,1); Wm(1) lambda / (n lambda); Wc(1) lambda / (n lambda) (1 - alpha^2 beta); for i 2:2*n1 Wm(i) 1 / (2*(n lambda)); Wc(i) 1 / (2*(n lambda)); end % 初始化 x_hat [0.3; 1.0]; P diag([0.1, 0.01]); Q diag([1e-4, 1e-4]); R 1e-3; % 滤波循环 for k 1:N % Sigma点生成 S chol((n lambda) * P, lower); X zeros(n, 2*n1); X(:,1) x_hat; for i 1:n X(:,i1) x_hat S(:,i); X(:,in1) x_hat - S(:,i); end % 预测 X_pred zeros(n, 2*n1); for i 1:2*n1 X_pred(:,i) X(:,i) dt * f(X(:,i)); end x_pred zeros(n,1); for i 1:2*n1 x_pred x_pred Wm(i) * X_pred(:,i); end P_pred Q; for i 1:2*n1 diff X_pred(:,i) - x_pred; P_pred P_pred Wc(i) * (diff * diff); end % 更新 Z_pred zeros(1, 2*n1); % 若量测为标量功率 for i 1:2*n1 Z_pred(i) h(X_pred(:,i)); end z_mean sum(Wm .* Z_pred, 2); Pzz R; for i 1:2*n1 diff_z Z_pred(i) - z_mean; Pzz Pzz Wc(i) * (diff_z * diff_z); end Pxz zeros(n,1); for i 1:2*n1 diff_x X_pred(:,i) - x_pred; diff_z Z_pred(i) - z_mean; Pxz Pxz Wc(i) * (diff_x * diff_z); end K Pxz / Pzz; x_hat x_pred K * (z_meas(k) - z_mean); P P_pred - K * Pzz * K; end这段代码里最容易出错的地方是权重初值的计算。beta参数在高斯分布下取2是常规经验值如果量测噪声不是高斯分布这个值就未必最优。alpha控制Sigma点偏离均值的程度太小会让协方差矩阵出现数值病态太大会损失高阶信息。我一般在alpha1e-3附近调试。另外一个细节是量测更新中Pxz的维度。如果量测是标量Pxz是n维列向量K是n维列向量代码里用矩阵乘法“/”是可行的但如果你把量测扩成m维向量Pxz会是n×m矩阵建议改用Pxz * inv(Pzz)或Pxz / Pzz的Matlab标准写法可以保持维度一致。我在运行的对比实验中同样的单机系统UKF的功角估计曲线在负荷突变后的动态过程中比EKF更平滑尤其在线性化误差被放大时表现更明显。当然UKF每一步需要多跑2n次函数调用循环体变大但状态维数只有2时耗时几乎感觉不到差异。4. 实操中的参数调优与关键细节4.1 协方差矩阵Q和R的设置经验做卡尔曼滤波Q和R几乎是所有调试问题的源头。我的基本思路是先固定R再调Q。R代表量测噪声方差如果PMU量测的精度已知直接用厂家给出的误差标准差平方即可不需要拍脑袋。Q则是模型可信度的衡量它和系统建模误差有关比如忽略调压器动态、负荷模型不准确都需要适当增大Q来“吸收”这些模型误差。初始值方面我常用的是Q对角线先设一个小量比如1e-4然后跑一次滤波观察新息序列也就是量测值与预测量测的差值。如果新息里存在明显相关性和均值偏移说明Q偏小或模型失配如果新息白噪声特性良好但状态估计抖得厉害说明Q偏大。这个调整过程很像在PID里看稳态误差和动态响应理论不多但手感很重要。另外状态变量之间的Q往往不是独立的。例如功角和转速本身存在物理耦合如果两个Q都取很小状态协方差就容易收缩得过快后续量测无法纠正状态导致滤波发散。我在多机系统上遇到过类似问题最后把转速对应的Q调到功角的5到10倍才恢复了跟踪能力。你也可以用对角线外的非零项去描述状态噪声相关性但调试难度会显著上升新手建议先保持对角结构。4.2 状态初始化与收敛性动态状态估计是递推算法初始状态给得不对结果通常很难看。我的做法是先用潮流计算或静态状态估计得到系统的初始运行点把发电机功角和转速的初值对应赋上。转速初值通常等于同步角速度这个偏差不会太大功角初值却可能差出几十度这时候P0不能给太小否则算法会自以为很确定后续量测很难把状态拉回真实轨迹。P0本质上是对初值不确定性的量化。我第一次跑EKF时把P0取得特别小结果滤波出来的功角一直保持初值附近明显跟不上仿真里的振荡就是因为初始协方差太小导致增益被压死。后来我把P0设成比如diag([0.1, 0.01])让滤波在前几十个采样点内快速收敛效果就正常了。还有一个工程技巧如果初始偏差很大可以先用一组真实量测做开环预测不去更新状态等方式运行几步后再启动滤波。这相当于用数据先把模型预测校正到合理区间。我见过有文献直接把这种方法叫“预热”在强非线性场景下非常实用尤其对EKF来说能显著降低线性化误差导致的早期发散风险。4.3 数值稳定性Cholesky分解与矩阵奇异问题UKF里的核心数值操作是协方差矩阵P的Cholesky分解因为Sigma点生成需要把P分解成下三角矩阵。实际运行中P经过多步递推可能失去正定性出现特征值接近0甚至为负的情况。这通常发生在状态可观测性不强或量测噪声非常小的时候。我常用的一个“保底”做法是在做Cholesky分解之前给P加一个很小的单位阵修正项比如P_reg P 1e-6 * eye(n)。这个修正会稍微增大Sigma点分布的协方差但能防止chol函数报错。还有一个更主流的方法是采用平方根UKF它直接传播P的Cholesky因子避免反复分解造成的数值舍入实现复杂度会高一些但值得在长期运行项目中推广。EKF那边虽然没有Cholesky但在更新步需要计算逆矩阵或者解线性方程组。如果量测噪声R取得特别小加上P_pred又比较接近奇异增益计算就会出现数值振荡。我在代码里习惯用K P_pred * H / (H * P_pred * H R)而不是显式求逆Matlab的“/”会走线性求解器比inv稳定得多推荐你也这样写。5. 常见问题与排查技巧实录5.1 滤波发散滤波发散是我碰到最多的现象表现很直观估计值早早偏出合理范围或者协方差矩阵迅速变成一个很小且不再变化的数值。这时你先别急着调代码按下面顺序排查。第一个嫌疑是雅可比矩阵错了。EKF对雅可比的要求很高哪怕某个元素的符号反了整个滤波都会乱。我建议在调试阶段用数值差分和解析表达式做一次对比打印几个采样点的Jacobian差异。第二种常见原因是Q设太小滤波器完全信任模型而模型本身有误差一旦状态被模型带到错误方向量测也拉不回来。这时候把Q调大几个数量级通常能立刻看到改善。第三个原因是初始P0给太少导致早期的卡尔曼增益过低。如果你发现滤波初期就发散把P0放大到diag([1, 0.1])之类的量级让滤波器有足够的“学习空间”。如果以上都不行就要怀疑模型本身有问题比如单机无穷大系统中的电磁功率公式符号反了或者时间步长过大导致离散化误差太大这就得回到建模层面检查。5.2 结果滞后或振荡滤波曲线跟真实轨迹相比总是慢半拍看起来像有相位滞后多半是Q取太小滤波器过度依赖模型预测量测更新对状态的修正作用被削弱。解决办法与发散问题相反适当增大Q对应状态的元素让滤波器更相信量测。如果你调了Q还滞后另一个方向是看R是否被设得过大量测被当成噪声滤掉了压低了更新增益。振荡则通常是另外一个极端Q过大导致每个量测都被高度重视估计值跟着量测噪声来回跳。这种情况下滤波曲线会有高频毛刺看起来像锯齿。你需要在Q和R之间找一个中间值。我的调试技巧是记录新息协方差的理论值与实际统计值理论上Pzz是预测的协方差实际统计可以从多步新息序列里计算出来如果实际远大于理论说明Q小如果实际远小于理论说明Q大。这样做等比用肉眼看好操作得多。对于单机模型最佳Q往往不是唯一的只要不引起发散状态跟踪效果就能接受。但到了多机系统各节点Q/R的比例关系会直接影响估计的空间分布那就需要更多调参经验。我的建议是先在单机上找到规律再迁移到多机会省掉很多摸索时间。5.3 EKF与UKF的计算量对比计算量方面我做过一个简单的Matlab计时实验状态维数取2量测取单个功率值跑10万步滤波。EKF单步耗时大约在零点几毫秒UKF因为每次要生成5个Sigma点并传播耗时大约是EKF的3倍左右。这个差距在单机模型上几乎可以忽略但如果你做的是几十台发电机的动态状态估计UKF的计算开销就会成为一个需要认真评估的因素。除了滤波器本身量测维数也会影响耗时。PMU量测一旦包含多个通道比如同时用功率、电压幅值和电压相角EKF里计算雅可比矩阵的开销会明显上升而UKF只是增加Sigma点传播后的量测聚合计算相对没那么敏感。这也是UKF在“多量测、多状态”场景下反而更划算的原因。我自己的实践感受是做学术型对比研究EKF和UKF都要实现两者都值得当成基线做工程落地先评估系统规模和非线性强度再决定用哪一套。如果你有实时性压力不妨先用EKF顶着数值上如果出现无法避免的发散再换UKF。最后分享一个小技巧我在Matlab里跑这套东西时习惯用Profiler分析耗时瓶颈。很多情况下真正拖时间的不是滤波算法本身而是重复计算的非线性函数里的三角函数。如果你能把sin、cos这些计算做向量化或者预计算某些不变参数性能提升比纠结选EKF还是UKF更明显。我个人在实际操作中的体会是先把单机EKF跑通再换UKF对比波形最后再研究多机扩展这条路走下来最稳。希望这些代码骨架和调参经验能帮你少走一些弯路。