去年做电力系统暂态分析课题的时候我需要一套能实时跟踪发电机动角轨迹的估计程序。一开始只想跑通扩展卡尔曼滤波EKF交差后来发现论文里大家都在用无迹卡尔曼滤波UKF做对比干脆两个版本都用Matlab写了一遍顺便把模型、代码骨架、调参过程和踩坑记录全部整理出来了。这篇东西就是那套代码的完整复盘适合做电力系统课程设计、毕业设计或者刚开始接触动态状态估计、想从公式落到Matlab代码的研究生和工程师参考。先说明一个观念问题很多人一提“电力系统状态估计”就想到加权最小二乘WLS那套方法是给SCADA稳态断面用的几秒甚至几十秒才能出一帧根本追不上暂态过程。电力系统动态状态估计用的是PMU提供的高频量测结合发电机转子运动方程把状态估计从“一张断面的快照”升级成“一条连续轨迹的追踪”。这个转变是所有后续代码逻辑的出发点。1. 为什么电力系统需要动态状态估计从静态潮流到实时轨迹1.1 静态状态估计的局限在哪静态状态估计处理的是某一时刻的冗余量测本质上在做潮流断面的拟合输出的是某一瞬间的电压幅值、相角、注入功率。它不关心系统从一个断面到下一个断面之间发生了什么。问题在于现代电网的动态行为越来越快新能源出力波动、电力电子变流器响应、连锁故障这些过程的特征时间只有几十到几百毫秒SCADA的扫描周期根本跟不上。PMU同步相量量测单元把采样率拉到了几十甚至上百赫兹量测数据里包含了机电暂态的完整过程。但要利用这些高频数据就必须引入系统动态模型也就是发电机的运动方程让滤波器沿着“模型预测量测修正”的递推路径去追踪状态变化。这就是动态状态估计的核心价值它能实时给出功角、转速这些无法直接量测的“内状态”为动态安全评估、广域保护、参数辨识和紧急控制提供输入。1.2 为什么只保留发电机经典二阶模型全阶电力系统模型是典型的微分代数方程组发电机有六阶甚至更精细的转子绕组动态还要考虑励磁、调速器、负荷模型。直接拿来做滤波器的状态方程状态维数会爆炸而且很多参数根本整定不了。我在这套代码里采用了发电机经典二阶模型也就是摇摆方程。这种做法的理由很工程动态状态估计关心的是机电暂态的时间尺度这个尺度下暂态电动势E′可以近似看成恒定值机械功率Pm在短时间窗内视为常数或缓慢变化负荷用恒阻抗等效后并入网络导纳矩阵。系统状态就被压缩为每台发电机的功角δ和转速ωn台发电机对应2n维状态。状态维数降下来之后EKF和UKF的矩阵运算量、参数整定难度都在可接受范围内。具体状态方程写成连续时间形式就是dδ_i/dt ω_i - ω_sM_i * dω_i/dt P_mi - P_ei - D_i * (ω_i - ω_s)其中M_i 2H_i / ω_s是惯性时间常数D_i是阻尼系数ω_s是同步角速度。这里最关键的电磁功率P_ei是状态的非线性函数P_ei E_i′² * G_ii Σ_{j≠i} E_i′ * E_j′ * [ G_ij * cos(δ_i - δ_j) B_ij * sin(δ_i - δ_j) ]G_ij和B_ij来自网络化简后的导纳矩阵。量测方程我取的是母线电压幅值和注入有功同样是关于功角的非线性表达式。这样整套系统的状态方程和量测方程都带有明显非线性正好用来考察EKF线性化近似与UKF采样逼近之间的性能差异。2. EKF在电力系统中的Matlab实现把雅可比矩阵当门槛2.1 从连续方程到离散化状态转移Matlab里做滤波第一步是把连续状态方程离散化。考虑到机电暂态的典型时间尺度我选了固定步长Δt 0.01s用最简单的欧拉法近似δ_k δ_{k-1} Δt * (ω_{k-1} - ω_s)ω_k ω_{k-1} Δt / M_i * (P_mi - P_ei(δ_{k-1}) - D_i * (ω_{k-1} - ω_s))这一步看似朴实无华但步长选择有讲究。Δt太长欧拉法对摇摆方程的截断误差会污染预测均值Δt太短PMU量测往往达不到那么高的采样率量测更新频率跟不上。0.01s是兼顾模型精度和工程可行性的常用选择。很多入门代码把状态转移写成线性矩阵那就失去了动态估计的意义。2.2 EKF预测-更新循环的Matlab骨架EKF的思想是在每个工作点把非线性函数做泰勒展开只保留一阶项得到一个近似的线性系统然后套用标准卡尔曼滤波框架。Matlab骨架大致长这样function [xHat, PPred, P] ekfPredictUpdate(xLast, PLast, z, sys) % 预测步 xPred state_fun(xLast, sys); [F, ~] jacobians(xLast, sys); PPred F * PLast * F sys.Q; % 更新步 [~, H] jacobians(xPred, sys); zPred h_fun(xPred, sys); S H * PPred * H sys.R; K PPred * H / S; xHat xPred K * (z - zPred); P (eye(length(xHat)) - K * H) * PPred; endstate_fun封装状态转移方程h_fun封装量测方程jacobians返回状态转移雅可比F和量测雅可比H。整套循环结构并不复杂真正的复杂度都集中在雅可比推导上。状态方程里电磁功率P_e对δ求导会得到一组三角函数线性组合的表达式量测方程对δ求导也类似。如果你用手推公式建议先在草稿纸上把n台机的功率方程展开再逐项对每个δ_j求偏导否则排错非常痛苦。2.3 用数值差分验证解析雅可比解析雅可比最大的问题是写错了不一定马上报错而是表现为滤波收敛慢、发散、或者在某些扰动场景下误差突然跳变。等你在一个大模型里定位错误时排查成本极高。我的做法是写一个check_jacobian.m用中心差分把F和H的数值解算出来和解析结果逐元素对比eps0 1e-6; for j 1:n xp x; xp(j) xp(j) eps0; xm x; xm(j) xm(j) - eps0; F_num(:, j) (state_fun(xp, sys) - state_fun(xm, sys)) / (2*eps0); H_num(:, j) (h_fun(xp, sys) - h_fun(xm, sys)) / (2*eps0); end disp(max(abs(F_num - F_ana), [], all)); disp(max(abs(H_num - H_ana), [], all));正常情况解析和数值差应在1e-6量级。我在第一版代码里把cos和sin的符号搞反了正是靠这个脚本半分钟就定位到了问题。这个习惯我后来一直保留就连只跑UKF的项目也会跑一遍这个检查——因为UKF虽然不需要解析雅可比但状态转移函数和量测函数本身的正确性仍然可以通过这个方式快速验证。EKF还有一层要意识到它只保留了一阶泰勒项本质上是在当前估计点用切线代替非线性函数。当故障后功角摆开到几十度、系统非线性被充分激发时一阶截断误差会直接体现为滤波滞后和偏差。这正好引出UKF。3. UKF的sigma点机制不求导也能捕捉非线性3.1 sigma点如何编码均值和协方差UKF的核心是无迹变换Unscented Transform。它的出发点很直观既然直接对非线性函数做线性化困难且精度有限不如找一组特殊的采样点——sigma点——让这组点的样本均值和协方差与原状态分布完全一致然后把每个点都丢进非线性函数里跑一遍再从跑完的结果中重新加权计算输出均值和协方差。这里可以打个比方EKF的做法是在弯曲的盘山公路上用当前点的切线方向代替整条路UKF的做法是路上取几个有代表性的位置每辆车实际开过去以后再统计整队车的位置分布。后者没有在函数曲线上做任何近似代价是多了几次函数计算。对于n维高斯分布状态选取2n1个sigma点λ α² * (n κ) - nχ_0 x̄χ_i x̄ (sqrt((n λ) * P))_i 第i列χ_{in} x̄ - (sqrt((n λ) * P))_i对应的权重为W_0^m λ / (n λ)W_0^c λ / (n λ) (1 - α² β)W_i^m W_i^c 1 / (2 * (n λ))这里sqrt表示矩阵平方根Matlab里通过chol((n λ) * P)实现。所有sigma点经状态函数传播后加权得到预测均值经量测函数传播后加权得到量测预测均值和协方差再算出状态与量测的互协方差就能像标准卡尔曼滤波一样计算增益并更新。3.2 UKF滤波循环的Matlab骨架UKF代码骨架和EKF高度相似只是把雅可比计算换成了sigma点生成与传播function [xHat, PPred, P] ukfPredictUpdate(xLast, PLast, z, sys) n length(xLast); lambda sys.alpha^2 * (n sys.kappa) - n; % 生成sigma点 sqrtP chol((n lambda) * PLast, lower); chi zeros(n, 2*n1); chi(:, 1) xLast; for i 1:n chi(:, i1) xLast sqrtP(:, i); chi(:, in1) xLast - sqrtP(:, i); end % 状态传播 chiPred zeros(n, 2*n1); for i 1:2*n1 chiPred(:, i) state_fun(chi(:, i), sys); end [xPred, PPred, Wm, Wc] weightedStats(chiPred, lambda, sys.alpha, sys.beta); % 量测传播 gammaPred zeros(length(z), 2*n1); for i 1:2*n1 gammaPred(:, i) h_fun(chiPred(:, i), sys); end zPred gammaPred * Wm; Pzz (gammaPred - zPred) * diag(Wc) * (gammaPred - zPred) sys.R; Pxz (chiPred - xPred) * diag(Wc) * (gammaPred - zPred); K Pxz / Pzz; xHat xPred K * (z - zPred); P PPred - K * Pzz * K; end核心开销在于2n1次状态函数和量测函数计算。n6时每步需要13次函数求值比EKF单次求值多一个量级但公式复杂度比手推雅可比低得多而且不用维护解析导数。3.3 alpha、beta、kappa三个参数怎么配UKF比较“劝退”新手的地方是参数整定。我的经验配置是α取1e-3到1e-2β2κ0。这套组合在大多数高斯假设下都能用也是很多文献的默认值。α控制sigma点离均值的远近。α太大sigma点离均值远低阶非线性下权重容易产生负方差α太小sigma点过于集中在均值附近对非线性传播的捕捉能力下降。β和先验分布有关。高斯分布下取2是理论最优非高斯时需要重新扫描。κ的经典选择是0也有人取3-n来消除高阶矩误差。对于6维状态3-n-3会让λ出现负值导致W_0^c为负虽然理论可行但数值上更容易碰到非正定问题。我一般就用κ0。建议把三组参数扫一个小网格对比RMSE变化趋势。你很快会发现α从1e-3调到1e-1估算误差变化远没有Q矩阵一个数量级的变化影响大。所以别在UKF参数上钻牛角尖后面那些坑才是真正的大头。4. 同一套算例下的实测对比EKF与UKF谁更值得4.1 算例与量测噪声的设定方式为了对比两个滤波器我在三机九节点系统上搭了一套测试基准。三台发电机对应6维状态状态量为[δ1, δ2, δ3, ω1, ω2, ω3]。仿真步长0.01s时长10s。故障设置为母线三相短路1.0s发生故障1.1s切除制造一段强非线性暂态过程。真实轨迹由同一套模型在不加量测噪声的情况下积分得到。给量测叠加的高斯白噪声设定为电压幅值噪声标准差0.01pu注入有功噪声标准差0.02pu相角噪声标准差0.005rad。滤波器初始状态取真实轨迹附近的一个小偏移初始协方差P0取0.01*eye(6)。这个设定很关键如果初始状态偏差过大EKF的第一轮线性化就会严重偏离直接发散。4.2 精度与计算开销的量化对比跑完10s仿真统计三段状态量稳态段、故障段、故障清除后暂态段的估计误差。典型结果如下指标EKFUKF功角RMSErad0.0210.008转速RMSEpu0.00420.0019最大功角估计偏差rad0.0780.031单步平均耗时ms0.320.51故障清除后功角摆开幅度较大系统非线性非常明显EKF切线和真实轨迹的偏差被放大功角估计在摆开初期明显滞后于真实值最大偏差接近0.08rad。UKF虽然也有滞后但sigma点传播更好地保留了非线性信息最大偏差只有EKF的不到一半。计算开销方面UKF单步耗时大约是EKF的1.6倍对于6维状态来说完全可接受。如果把扰动幅度减小只做负荷阶跃或者正常运行工况两者的RMSE差距会迅速缩小到同一数量级。这说明EKF在小扰动场景下足够用UKF的精度优势主要体现在大扰动、强非线性阶段。4.3 什么时候直接用EKF别折腾UKF我给自己定的选择标准很简单如果只是做正常运行状态追踪或者量测更新频率很高每次修正幅度很小EKF够用且代码量少。但如果算例里有短路、切机、大负荷突变这类强扰动或者你不想手工推导高阶雅可比UKF是更稳妥的选择相当于用多一点计算量去换鲁棒性。另一个实用建议先跑通EKF再把循环替换成UKF。EKF的调试链路短、报错逻辑清晰一旦模型函数和量测函数验证无误切换到UKF通常只需要改滤波循环那几十行代码出问题的概率小得多。5. Matlab工程实现中五处容易翻车的细节5.1 参数单位与基准值明明模型对却全错标幺制是电力系统仿真的默认语言但混用有名值和标幺值的情况非常常见。比如转速ω在状态方程里用有名值rad/s和用标幺值pu对应的M 2H / ω_s差别巨大。我曾见过有人把阻尼系数D从某个数据源直接搬过来那个数据源用的是有名值而代码里是标幺制滤波结果直接发散。我的建议是模型内部全部用标幺制所有外部数据在接口处统一转换。ω_s在标幺制下就是1但M依然是2H / 314.16如果H以秒为单位而有名同步角速度是314.16rad/s。把这条规则写进代码注释里能帮你省下大量排查时间。5.2 初值X0与P0滤波器的起跑姿势动态估计的初值一般来自静态状态估计或潮流计算结果X0不需要精确但方向要大致对。P0代表初始不确定度它的量级比数值本身更敏感。P0太小滤波器开局“过度自信”前几步量测更新幅度极低真实轨迹跑出去了还追不回来P0太大初始阶段量测噪声被完全吸收功角曲线前期剧烈抖动。我的经验是从P0 1e-2 * eye(2n)起步观察前50步的估计曲线。如果前期抖动大下调到1e-3如果收敛太慢上调到0.1。这个扫描过程十分钟就能完成不要嫌麻烦。5.3 Q和R矩阵整定哪些能调哪些不能动Q表示模型误差R表示量测误差。这两个矩阵本质上是滤波器的“信任分配器”Q给得大说明你觉得模型不可靠更多相信量测R给得大说明你觉得量测噪声大更多相信模型。电力系统动态估计里R往往可以由PMU的精度指标直接定下来属于“不能随便动”的参数。Q则没有清晰标定途径基本靠调。我的整定套路是固定R不变让Q q * eye(2n)q从1e-6到1e-2按对数扫描比较各组的功角RMSE。通常RMSE随q变化的曲线呈U形中间总有一个稳定平台区在平台区取值最稳妥。注意Q太小比Q太大更危险Q过小时滤波器对量测修正不敏感一旦真实动态偏离模型假设误差会一直累积到发散。5.4 Cholesky分解失败与数值稳定性UKF每步都要做chol((n λ) * P)P必须是正定对称矩阵。浮点累积会让P慢慢失去对称性甚至出现小负特征值Matlab直接报错。这种问题在仿真时间较长或量测更新很频繁时容易出现。处理办法有两个一是每步强制对称化P (P P) / 2二是给P的对角线加一个极小量eye(n) * 1e-9。我自己在代码里两个措施都加了也只加在UKF分支不影响EKF。这个保护性修改对结果精度的影响可以忽略但能让你免遭跑两小时程序然后在最后一步崩溃的体验。5.5 快速定位发散的三板斧滤波器发散是最让人崩溃的调试场景。我总结了一套固定的排查顺序按这个顺序来基本半小时内定位问题第一板斧打出新息序列。新息是量测值减去量测预测值即z - zPred。如果它均值明显偏离零说明系统模型偏差或量测偏差而不是滤波器自身的数值问题。第二板斧同时画出预测状态和更新状态。如果更新状态相对预测状态出现剧烈跳变说明量测噪声过大或R太小如果预测状态本身已经偏离真实轨迹问题出在状态方程或Q上。第三板斧把量测噪声设为0R趋近0跑一遍理想情况。如果这时候还是发散说明雅可比、状态方程或量测函数有硬错误和滤波器参数无关。这个操作能快速把“参数问题”和“模型问题”分开。这三板斧本质上是把滤波器当成一个黑盒系统逐步缩小可疑范围。我每次调试新系统都会按这个流程走一遍效率比对着矩阵数值瞎猜高得多。这套EKF/UKF代码跑通之后我自己最大的体会是滤波器本身不难难的是把电力系统模型准确翻译成Matlab函数。你只要把单位制、初始协方差、Q/R整定这三件事理顺EKF和UKF切换只是换几十行循环的事。后续如果你想把代码扩展成迭代扩展卡尔曼IEKF、平方根UKF或者处理PMU量测时延、做分布式多机估计底座都不用动太多重点扩展量测模型和滤波循环即可。