在电力系统的在线监测与运行控制里动态状态估计一直是个绕不开的核心话题。用扩展卡尔曼滤波EKF和无迹卡尔曼滤波UKF去跟踪发电机功角、转速这些动态状态是目前学术研究和工程尝试里最主流的做法。这套方案落地成Matlab代码以后既能跑通仿真验证算法效果又能为后续接入PMU实测数据打底。这篇博文围绕这套代码的实现思路、算法选型、建模细节和排坑经验展开凡是正在做电力系统动态状态估计、或者准备在Matlab里实现卡尔曼系列算法的朋友都能从中找到可以直接参考的路径。1. 动态状态估计为何是电力系统的实时眼1.1 静态状态估计的短板数据与模型的滞后传统电力系统状态估计以加权最小二乘为代表它处理的是某一断面下的静态快照利用SCADA系统采集的遥测、遥信数据解算出母线电压幅值和相角。这套体系运行了几十年在稳态工况下表现尚可但面对现在的电网形态短板越来越明显。SCADA的数据刷新率通常在几秒到几十秒一次量测到达时间不严格同步。当系统发生负荷快速波动、新能源出力突变、短路故障后的暂态过程SCADA拿到的数据往往已经过时了。再加上传统状态估计模型以代数方程为主不包含发电机转子运动这类动态方程它给出的结果本质上是一个准稳态解无法反映功角、转速在动态过程中的演化轨迹。换句话说它在时间尺度上就达不到动态追踪的要求。有人可能会说既然有PMU同步相量测量单元了不是可以直接测功角、转速吗问题是PMU只能测电气量电压、电流相量功角和转速属于机械量PMU并不能直接完全测量。虽然可以间接推导但实际中量测噪声、线路模型误差、角度基准漂移都会让直接推导的结果不可靠。所以需要动态状态估计把系统模型和PMU高精度量测融合起来通过滤波算法实时估计内部动态状态。1.2 动态状态估计的核心任务在线追踪发电机的运动轨迹动态状态估计的目标很明确基于发电机转子运动方程、励磁绕组动态方程等状态方程结合PMU提供的快速量测在每一个时间步递归估计出发电机的动态状态。最具代表性的状态量包括转子功角δ、转速偏差ω、暂态电动势的d轴和q轴分量对应发电机四阶模型或更详细的状态。这样做的工程价值很直接。调度和稳控系统需要知道当前发电机是否处于稳定运行范围功角摆开有多大转速偏差是否收敛当系统受扰动后整个过程是否向稳定方向演化。如果这些状态能被实时、平滑、去噪地估计出来后续的紧急控制、低频减载、失步预测就有了一组可靠的输入信号。也正因为如此动态状态估计被普遍看作电力系统在线动态安全分析的数据前置环节。要完成这个任务滤波器选型就非常关键。卡尔曼滤波框架天然适合状态方程预测量测更新的递归模式而系统中的状态方程和量测方程都是非线性的这就引出EKF和UKF两个经典选项。2. EKF和UKF的原理与选型逻辑2.1 EKF对非线性做一阶线性化代价是精度和稳定性扩展卡尔曼滤波的思路非常直接在每一个时间步把非线性状态方程和量测方程在当前估计点附近做一阶泰勒展开用雅可比矩阵代替线性卡尔曼滤波中的状态转移矩阵和量测矩阵然后套用标准的预测-更新流程。预测阶段的核心计算是状态预测x_pred f(x_est)其中f是离散化的机电暂态状态方程协方差预测P_pred A * P_est * A QA是状态方程对状态的雅可比矩阵卡尔曼增益K P_pred * H * (H * P_pred * H R)^(-1)H是量测方程对状态的雅可比矩阵状态更新x_est x_pred K * (z - h(x_pred))协方差更新P_est (I - K * H) * P_pred。EKF看似简单实战中却有几处暗坑。电力系统状态量的数量级跨层很大功角以弧度计通常在零点几到几之间转速偏差以标幺值或rad/s计范围也在零点几量级但不同状态对应的雅可比矩阵元素可能相差几十倍到上百倍。这会造成协方差矩阵的病态影响数值稳定性。此外EKF只是用一阶泰勒近似逼近非线性函数系统非线性强比如故障后的剧烈摆动阶段时线性化误差会被放大滤波结果可能明显偏离真值甚至发散。2.2 UKF用一组sigma点绕开求导逼近非线性传播无迹卡尔曼滤波的核心思想是无迹变换。它不再对非线性函数做泰勒展开而是选取一组带权值的sigma点让这些点经过非线性函数传播后再用加权统计方法重建均值与协方差。具体说如果状态维度是n那么选取2n1个sigma点用尺度参数调整点在均值附近的散布程度其中常见的参数配置包括α决定sigma点的散布通常取1e-3~1之间、β用于融入先验分布信息高斯分布取2、κ通常取0或3-n。这些sigma点通过状态方程和量测方程传播后用它们计算预测均值和预测协方差再相同框架下计算卡尔曼增益并更新状态。UKF最直接的好处是不需要手推雅可比矩阵。对于电力系统这种状态方程和量测方程形式复杂、含三角函数和代数变量消除过程的场景避免求导是节省人力和降低出错率的巨大优势。在非线性程度较高的暂态过程中UKF一般能比EKF获得更高的估计精度因为它对非线性函数的传播精度能达到三阶矩水平对高斯分布而言而不是像EKF那样只保留一阶。代价也很明确计算量大约是EKF的2n1倍因为每个时间步需要多次求值非线性函数。对n4~9的发电机动态状态估计来说这个计算量在Matlab仿真和现代硬件上完全可接受所以UKF的性价比相当高。2.3 两个放一起比较什么时候选谁下表是我在实际项目里总结的选型对照对比维度EKFUKF是否需要雅可比矩阵需要手推或符号计算不需要非线性逼近精度一阶强非线性时偏大近似三阶强非线性表现更好计算量低约2n1倍函数求值实现难度中难在求导和调雅可比低只需正确设置sigma点参数协方差数值稳定性可能因雅可比病态恶化需要保证协方差正定性模型更换的灵活性模型一变导数重推直接改状态函数即可我的建议是如果你刚开始做动态状态估计先用电台模型把UKF实现跑通因为代码链路短、不容易在求导环节卡住如果论文或工程方案中同时需要对比方法再把EKF补上两者对比也能更直观地凸显算法差异。这套Matlab代码把两种滤波器都实现了正好可以同台对比。3. 系统建模发电机动态模型与量测模型3.1 用几阶发电机模型如何离散化动态状态估计的模型选择精度和复杂度要平衡。经典二阶模型摇摆方程只描述转子运动状态变量是功角δ和转速偏差ω结构最简单适合快速验证滤波器代码是否正确收敛。更贴合实际的是四阶模型在δ、ω基础上增加d轴暂态电动势Ed和q轴暂态电动势Eq描述励磁绕组和阻尼绕组的动态能反映暂态过程中的电压变化。以四阶模型为例连续时间状态方程可写成dδ/dt ω * ω_b - ω_s不同文献量纲写法有差异另一种常用写法是dδ/dt ω - ω_sdω/dt (Pm - Pe - D*(ω - ω_s)) / (2H)dEq/dt (-Eq Ef (Xd - Xd)*Id) / TdodEd/dt (-Ed - (Xq - Xq)*Iq) / Tqo需要注意这里电气量Id、Iq、Pe并非状态变量的直接函数还要通过代数网络方程消去。要将这些连续方程变成离散状态空间形式最简单的做法是采用一阶欧拉离散它直观易实现但步长小时精度才够。实际中我建议采用四阶Runge-Kutta离散这样在PMU量测步长通常20ms到100ms下仍能保持足够的精度。3.2 量测方程怎么构造量测方程建模是这套代码里最容易出问题的地方。以机端PMU量测为例假设量测量为机端电压幅值Vt、机端电压相角θ或者还可以加发电机输出有功功率Pe和无功功率Qe。量测方程的形式并不唯一取决于你选择哪个坐标系和母线方程。如果机端电压相角θ与功角δ之间的角度差作为变量Y θ - δ那么量测方程就可以写成关于状态量和端电流的非线性方程Vt sqrt(Vd^2 Vq^2)θ δ atan(Vq / Vd)具体符号需根据参考轴定义调整这样量测方程天然是非线性的而且量测矩阵HEKF需要推导起来比状态方程更繁琐因为要经过潮流网络方程间接求导。这也是很多人在EKF实现中被卡住的环节后面排查章节我会给出验证手段。3.3 量测噪声、过程噪声和初值怎么定卡尔曼滤波里Q和R矩阵的取值直接影响滤波性能。R矩阵可以从PMU的技术参数中估计例如电压幅值误差0.1%~0.2%相角误差0.01~0.02弧度把它对角化设置为量测协方差矩阵即可。Q矩阵则复杂得多它代表过程噪声本质是对模型误差的容忍度。模型省略了励磁调节器、调速器动态这些未知输入都会以过程噪声的形式体现。初值的设置基于稳态潮流解先跑一次潮流计算得到发电机端电压和功角的初始稳态值并以该稳态点计算其他状态初值。协方差P_est初始值通常设置为一个对角矩阵数值不宜太小若设置过小会让滤波器过于自信导致初期的量测更新跟不上真值变化太大又会让前几步估计波动明显。一般取各状态初值平方的一定比例即可。4. Matlab代码实现从框架到核心函数4.1 工程化的模块划分写Matlab代码不要一个脚本塞到底我建议按模块划分主脚本用于设置参数、生成仿真数据、调用滤波器并绘图系统模型函数负责状态方程和量测方程滤波器函数独立封装EKF和UKF。这样不同算例之间切换只需要修改参数代码复用性高。主脚本中需要定义开关量控制滤波器类型例如用methodekf或methodukf选择调用哪个函数最终统一输出估计状态序列方便对比。4.2 状态方程和量测方程的Matlab封装状态方程函数建议写成统一形式[x_next] dst_func(x, u, dt, para)其中u表示输入如机械功率Pm、励磁电压Efpara结构体放发电机和时间常数等参数。量测方程函数则写成[z] meas_func(x, para)。两个函数需要保持输入输出接口固定因为UKF的sigma点传播和EKF的预测阶段都要反复调用它们。四阶Runge-Kutta离散实现的关键代码大致是function xn f_rk4(x, u, dt, para) k1 f_cont(x, u, para); k2 f_cont(x 0.5*dt*k1, u, para); k3 f_cont(x 0.5*dt*k2, u, para); k4 f_cont(x dt*k3, u, para); xn x (dt/6)*(k1 2*k2 2*k3 k4); end其中f_cont是连续时间状态方程负责根据当前状态和输入计算导数。4.3 EKF核心代码逻辑EKF的关键点在于雅可比矩阵的计算。这里给出两种方式一是解析法用符号工具箱求出A和H二是数值差分法用有限差分近似。后者实现简单但要注意步长选择太大会引入截断误差太小会淹没在浮点误差中。核心更新代码如下function [x_est, P_est] ekf_update(x_est, P_est, z, u, dt, para) % 预测 x_pred f_rk4(x_est, u, dt, para); A jacobian_state(x_est, u, dt, para); P_pred A * P_est * A para.Q; % 量测预测 z_pred meas_func(x_pred, para); H jacobian_meas(x_pred, para); S H * P_pred * H para.R; K P_pred * H / S; % 更新 x_est x_pred K * (z - z_pred); P_est (eye(length(x_est)) - K * H) * P_pred; end这里我特意加了数值差分辅助函数当你对解析导数没有十足把握时先用有限差分结果对比验证function A jacobian_state(x, u, dt, para) n length(x); A zeros(n, n); f0 f_rk4(x, u, dt, para); h 1e-6; for i 1:n x_pert x; x_pert(i) x_pert(i) h; f_pert f_rk4(x_pert, u, dt, para); A(:, i) (f_pert - f0) / h; end end4.4 UKF核心代码逻辑sigma点生成与权重计算UKF的实现关键在无迹变换和权重计算。按常用的比例修正形式权重计算如下function [chi, Wm, Wc, c] ut_transform(x, P, alpha, beta, kappa) n numel(x); lambda alpha^2 * (n kappa) - n; c n lambda; % 计算协方差平方根 sqrtP chol((n lambda) * P, lower); chi zeros(n, 2*n 1); chi(:, 1) x; for i 1:n chi(:, i1) x sqrtP(:, i); chi(:, ni1) x - sqrtP(:, i); end Wm zeros(2*n 1, 1); Wc zeros(2*n 1, 1); Wm(1) lambda / (n lambda); Wc(1) lambda / (n lambda) (1 - alpha^2 beta); for i 2:2*n 1 Wm(i) 1 / (2*(n lambda)); Wc(i) Wm(i); end end注意这里用了chol函数求平方根后续如果遇到矩阵非正定导致chol失败的问题参见第5章排查办法。UKF预测和更新阶段的代码逻辑是把每个sigma点分别通过状态方程和量测方程传播然后加权统计均值和协方差。4.5 仿真数据生成让真值有处可查做滤波器最怕没有参照物。这套代码里我把真实系统用精细的仿真模型生成状态方程中不加过程噪声、量测方程输出后叠加高斯白噪声模拟PMU量测。这样既能得到noisy量测序列z又能保留干净的状态真值x_true用于计算估计误差。生成方式类似for k 1:N x_true(:, k1) f_rk4(x_true(:, k), u(:, k), dt, para); z(:, k) meas_func(x_true(:, k1), para) mvnrnd(zeros(1, m), R); end然后以z作为滤波器的输入比较x_est与x_true。常用的性能指标包括RMSE均方根误差和最大绝对误差还可以记录单次滤波的总耗时。4.6 参数整定与初始化示例下面是我在单机无穷大算例中的一组典型参数可供参考发电机惯性常数H 3.5s阻尼系数D 2.0同步电抗Xd 1.8暂态电抗Xd 0.3Xq 1.7Xq 0.45励磁绕组时间常数Tdo 7.5s阻尼绕组时间常数Tqo 0.45s量测噪声电压幅值标准差0.005相角标准差0.01弧度Q矩阵对角元取1e-4到1e-3量级R矩阵按上述噪声标准差平方设置仿真时长5s步长0.01sPMU量测每0.02s采一次50Hz。这套参数跑下来我见过UKF的功角RMSE在0.001~0.01弧度量级转速RMSE在1e-4量级具体数值取决于故障场景强度。5. 实测中遇到的问题与排查技巧实录5.1 滤波器发散最常见的翻车现场动态状态估计里最让人头疼的就是滤波发散前几步看起来正常忽然某一步估计值突变之后误差越来越大直接把曲线跑飞。根据我调试的经历原因通常归为几类。一是Q矩阵设置不合理。Q太小意味着过程模型被过度信任一旦模型与实际系统存在偏差滤波器会不断把误差归结到量测上输出出现震荡Q太大会导致滤波过于依赖量测失去平滑去噪能力。排查办法是先固定R用一组对比仿真扫描Q的数量级观察RMSE曲线变化。二是初值严重偏离真值。在仿真中可以把初值从真值偏移10%~20%测试滤波器随时间收敛的能力。若发散优先检查状态方程离散化是否稳定尤其二阶模型里ω的量纲很容易写错。攻角对时间的导数到底是ω还是ω-1标幺值基准下的转速偏差必须要一致否则状态演化方向都是错的。三是协方差阵数值病态。电力系统状态量量纲跨度过大时P矩阵条件数可能达到10^8甚至更高EKF中尤甚。解决办法是归一化处理功角用弧度转速直接用rad/s的量纲不要混用标幺值和其他量纲系统或者对P做定期对称化和特征值裁剪把极小负特征值修剪到0。5.2 Jacobian矩阵算错EKF的隐形炸弹EKF对雅可比矩阵的准确性极其敏感。H矩阵算错一个符号结果可能不是数值偏差而是直接发散。最典型的错误集中在量测方程对功角的导数上因为量测方程不仅含有状态变量的三角函数项还要通过电流Id/Iq间接依赖状态这层复合求导极易漏项。我建议的验证方法是有限差分对照法先写解析H矩阵再用数值差分生成参考矩阵两者对比误差在1e-4以内基本可以认定正确。这个验证可以做成一次性辅助脚本不放进主循环避免影响性能。还有个小技巧当EKF和UKF在同样的条件下UKF正常而EKF发散优先怀疑雅可比矩阵而不是滤波器结构因为UKF不依赖雅可比。5.3 UKF的协方差非正定问题UKF虽然避开了雅可比求导但引入了另一个典型问题协方差矩阵P在递推中可能失去正定性导致chol分解失败报错信息通常是Matrix must be positive definite。这类问题多出现在初始P矩阵设置不当、量测噪声过小导致S矩阵接近奇异、或者步骤中浮点误差累积。应对方案在每次更新后强制对称化P (P P) / 2诊断矩阵最小特征值若接近0或为负加一个很小的正对角阵比如1e-8 * eye(n)作为正则项改用平方根UKFSR-UKF它直接递推协方差的Cholesky因子从根上避免非正定问题代价是实现更复杂。检查参数kapppa当n较大时取kappa3-n会导致lambda为负传播中可能出现负权值这也是数值不稳定的来源之一。建议使用GTFGaussian-to-First参数组合即alpha1beta0kappa3-n或者经典alpha1e-3beta2kappa0的组合并确认lambda满足必要的条件。5.4 离散化步长和采样周期怎么选PMU上送速率多为10~50帧每秒而发电机暂态过程的时间常数从毫秒级次暂态分量到秒级功角摆动都有。离散化步长要能覆盖快变过程否则系统微分方程的数值解本身就发散。我的经验是仿真/状态预测的积分步长取0.001~0.01s量测更新步长取0.02~0.1s两者之间可以不等距——积分步长小时量测更新之间做多个积分小步等更新时刻到了再用当时的量测做校正。这样滤波器结构更贴近实际也利于对比不同PMU速率对估计精度的影响。感兴趣的话还可以做一个量测丢失率测试模拟PMU丢帧场景看看滤波器在缺少某些时刻量测时是否能靠模型预测继续维持基本可用。这也是动态状态估计走向实际应用时不得不面对的问题。6. 算例对比EKF与UKF的实测表现6.1 同一场景下的性能对比结果以某单机无穷大系统为算例在1s时设置一个三相短路故障、1.15s切除的暂态场景分别用EKF和UKF做动态状态估计我测试得到的结果具有明显的代表性。EKF在功角摆动幅度较大、系统非线性较强的时段出现了明显的估计偏差最大功角估计误差约0.02弧度转速误差约2e-3 rad/s。而UKF在整个过程中保持了更平滑的估计轨迹功角最大误差约0.005弧度转速误差约5e-4 rad/s精度高了一个量级左右。计算耗时方面n4时UKF每步需要9次函数求值耗时约为EKF的3~5倍因为EKF每次迭代还需要计算雅可比所以倍数低于理论预期的2n1。当然在系统非线性很弱、接近稳态运行时EKF与UKF的估计结果几乎重合这时候UEKF在计算速度上的优势就体现出来了。这就印证了选型阶段的结论强非线性暂态过程选UKF轻载稳定场景EKF性价比更高。6.2 面对量测噪声差异时的鲁棒性我还做过一组量测噪声敏感度测试把量测噪声水平从0.005逐步加大到0.05观察两种滤波器的RMSE变化。结果是UKF的性能下降更缓慢在量测噪声较大时依然能维持可用的估计质量EKF则对噪声的容忍度明显更低噪声大时滤波轨迹出现较明显的跟随滞后。这背后的原理是UKF在更新步骤中用的协方差传播更准确对量测异常值的反应更平缓。6.3 从单机走向多机系统时要注意什么当从单机系统扩展到多机系统比如IEEE 9节点、39节点状态维度上升状态方程和量测方程之间的耦合变复杂量测方程需要包含互联母线的电压相量信息。此时有两个建议一是状态可分割为每台机一个局部滤波器分散式动态状态估计局部量测覆盖本机及相邻母线降低全局维度适合并行计算二是如果坚持集中式全状态估计务必关注P矩阵维度增大后的数值稳定性Q矩阵也要相应调整因为不同发电机之间过程噪声相关性需要建模否则多机系统背景下滤波器容易误发散。7. 扩展方向这套代码还能往哪走做滤波算法和仿真代码最忌讳的就是跑通一次就封存。这套EKF和UKF代码框架其实可以平滑迁移到更多场景。叶片尺度上可以替换发电机模型例如加入励磁调节器AVR和调速器Governor的简单动态模型状态维度扩展到8~10维观察动态状态估计在闭环控制影响下的表现。再进一步把光伏、储能逆变器的动态模型纳入估计范围研究新型电力系统下多时间尺度的状态估计问题这在当前碳中和背景下很有研究价值。算法层面可以对比更先进的滤波方案中心差分卡尔曼滤波器CDKF、平方根UKF、容积卡尔曼滤波CKF还有粒子滤波和H∞滤波这些方法在处理更强的非线性和非高斯噪声时各有擅长。动态状态估计结合深度学习也是热门比如用神经网络自动调整Q、R矩阵或者用RNN学习模型误差残差作为过程噪声补偿都能在现有Matlab框架上扩展。把PMU实测数据接入这套代码也不难只需把量测数据来源从仿真生成换成数据文件读取注意时间对齐和标幺值换算即可。到这一步代码就不再只是算法验证工具而是具备工程原型潜质的应用底座了。回头再看整个动态状态估计项目最核心的心得其实只有两条第一建模一致性比算法选择更重要状态方程和量测方程一旦在量纲、参考系、离散化上出错什么滤波器都救不回来第二Q和R的整定不是一次性的配合不同场景反复调参才能真正体会每种滤波器的性格。这套Matlab代码作为起点足够你在EKF和UKF之间来回切换、设置各种故障场景并评估性能了。踩过那些发散的坑、调过参数之后你对动态状态估计的理解会扎实很多。