1. 项目概述为什么要做电力系统动态状态估计电力系统动态状态估计这个话题在学术圈和工程圈里一直是热点但也是不少人的拦路虎。很多同学拿到题目第一反应是卡尔曼滤波我懂EKF和UKF我也看过原理可一到电力系统里那些矩阵、非线性函数、量测方程到底怎么对应一下子就懵了。先说清楚这个项目到底是干什么的。电力系统运行过程中发电机转速、功角、母线电压幅值和相角这些状态量一直在动态变化。我们不可能在每个节点都装量测装置去实时读取全部状态而且PMU同步相量测量装置和SCADA数据采集与监控系统的采样周期不一样、噪声也难以避免。动态状态估计要解决的问题就是利用已有的量测数据结合系统模型把那些测不到或者测不准的实时状态尽可能准确地估计出来。这里为什么非要动态传统的静态状态估计加权最小二乘法那一套只能给出某个断面下的稳态状态没法跟踪系统的动态演化过程。当系统发生扰动或者负荷快速变化时静态估计就力不从心了。动态状态估计在静态的基础上增加了一步用上一时刻的状态去预测下一时刻的状态再用量测数据去修正预测值。这一步“预测-修正”正是卡尔曼滤波家族的核心思想。之所以选EKF和UKF做对比而不是直接用经典的线性卡尔曼滤波KF是因为电力系统本质上是一个强非线性系统。同步发电机的转子运动方程、电磁功率方程、潮流方程全都是非线性的。线性KF假设系统是线性的、噪声是高斯的拿到电力系统里直接套用必然出现模型失配。EKF和UKF是在非线性条件下处理状态估计的两个主流路线EKF靠泰勒展开做线性化UKF靠一组Sigma点直接传播概率分布。两者各有优缺点放在同一个框架下对比既能验证各自算法的有效性也能在实际应用中根据工况选型。这套代码适合谁如果你是电力系统方向的研究生或者做控制理论、状态估计相关工作的工程师这代码可以直接作为基线来改。代码里实现了从系统建模、EKF/UKF算法到误差分析的全流程拿到手跑通DEMO再结合自己课题的模型替换进去就能快速产出对比结果。2. 核心原理EKF和UKF在电力系统里的对应关系2.1 发电机动态模型与状态空间表达要理解状态估计得先搞清楚被估计的对象是谁。电力系统动态状态估计最经典的被估计对象是同步发电机。在单机无穷大系统或者多机系统中发电机的三阶实用模型即考虑励磁动态的模型是常用的折中方案。以单机无穷大系统为例状态向量通常是x [δ, ω, E′q]^T其中δ是发电机功角转子相对于同步旋转参考轴的角度ω是转子角速度偏差E′q是q轴暂态电动势。微分方程如下dδ/dt ω - ωs dω/dt (Pm - Pe - D(ω - ωs)) / (2H) dE′q/dt (Efd - E′q - (Xd - XdI)Id) / Tdo这里的Pm是机械功率Pe是电磁功率D是阻尼系数H是惯性时间常数Efd是励磁电压Tdo是励磁绕组时间常数。而Pe的计算又涉及网络方程它和E′q、功角δ之间存在非线性关系在单机无穷大系统里可以写成Pe E′q·V∞·sinδ / XdΣ其中XdΣ是发电机暂态电抗和线路电抗的等效和V∞是无穷大母线电压。这样一来状态方程f(x)就是一个非线性向量函数。这组方程是EKF和UKF的公共底座两种滤波都是围绕这个非线性状态方程展开的。再说量测方程。在实际工程中我们能直接观测到的是发电机的输出有功功率Pe、机端电压幅值Vt、功角δ通过PMU设备可以测到等。量测方程h(x)也呈现非线性特征比如机端电压幅值与E′q、δ的关系就不是简单的线性叠加。量测噪声通常假设为零均值高斯白噪声协方差矩阵R由量测设备的精度决定。PMU的量测误差一般在0.5%左右这在仿真中可以折算成对应的标准差。把状态方程和量测方程合并起来就得到标准的非线性离散系统x(k1) f(x(k)) w(k) z(k) h(x(k)) v(k)其中w(k)是过程噪声代表模型本身的误差v(k)是量测噪声。这两个噪声的统计特性直接决定了滤波器的收敛速度和估计精度后面仿真部分有实际验证。2.2 扩展卡尔曼滤波EKF线性化的艺术与代价EKF的思路非常直观既然系统是非线性的那就在当前估计值附近做一阶泰勒展开把非线性函数线性化然后套用标准KF的框架。具体到实现EKF在每个时刻需要计算两个雅可比矩阵状态转移矩阵F(k) ∂f/∂x在x̂(k|k)处求导和量测矩阵H(k) ∂h/∂x在x̂(k1|k)处求导。这两个矩阵的计算是整个EKF的难点也是容易出错的地方。电力系统的发电机模型里f对状态变量求导会得到很多含三角函数的项一个不小心符号错了整个滤波就发散。我在代码里用符号工具箱辅助推导再手工核对最终确定了解析表达式这样精度更高、运行时也不怕数值微分带来的额外误差。EKF的预测-修正流程跟标准KF一致预测步x̂(k1|k) f(x̂(k|k))P(k1|k) F(k)P(k|k)F(k)^T Q修正步K(k1) P(k1|k)H(k1)^T [H(k1)P(k1|k)H(k1)^T R]^(-1)x̂(k1|k1) x̂(k1|k) K(k1)[z(k1) - h(x̂(k1|k))]P(k1|k1) [I - K(k1)H(k1)]P(k1|k)EKF的优势是计算量小三个状态量的情况下雅可比矩阵只有3×3矩阵运算负担很轻适合实时性要求高的场合。但它的劣势也很明显一阶线性化在大非线性强度下有截断误差可能导致估计偏差甚至滤波器发散。在电力系统重负荷、大扰动工况下状态量变化剧烈时这种线性化误差会被放大。2.3 无迹卡尔曼滤波UKF不线性化用Sigma点直接传播UKF的出发点跟EKF完全不同它不试图去线性化非线性函数而是近似状态的概率分布。一个高斯分布可以用一组确定性的采样点称为Sigma点来表征这些点经过非线性函数传播后再用加权统计的方式得到均值和协方差。这种做法在数学上能保证至少二阶精度而EKF只有一阶精度。UKF的Sigma点生成方式是核心。对于n维状态向量选取2n1个Sigma点χ0 x̄ χi x̄ (√((nλ)P))i, i 1,...,n χi x̄ - (√((nλ)P))i, i n1,...,2n其中λ α²(nκ) - nα决定Sigma点在均值周围的扩散程度通常取1e-3到1κ是次级缩放参数通常取0或3-nβ用于合并先验分布的高阶信息高斯分布取2最优。这些参数看着繁琐实际仿真中对滤波性能的影响非常大后面有专门的参数调参心得。每个Sigma点经过非线性状态方程f传播后得到一组预测点再通过加权求和得到预测均值和协方差χi(k1|k) f(χi(k|k))x̂(k1|k) Σ Wi^m χi(k1|k)P(k1|k) Σ Wi^c [χi(k1|k) - x̂(k1|k)][...]^T Q量测更新的类似也是把Sigma点通过量测方程h传播后做加权统计。UKF全程不需要计算雅可比矩阵这在系统模型复杂、导数推导困难时是极大的优势。代价是计算量比EKF大需要反复传播2n1个点但在现代计算机上这点开销完全不是问题。在电力系统中等规模模型比如10机39节点系统中UKF的高精度优势会体现得越来越明显。不过UKF对非高斯分布的适应性也有一个隐含假设它依然假设后验分布是高斯的只是对抗非线性传播的能力比EKF强。如果真的遇到强非高斯噪声UKF和EKF都会受影响那就得考虑粒子滤波PF了。3. 仿真实现Matlab代码框架与关键细节3.1 仿真场景设置与参数初始化首先明确仿真场景。这里我搭建的是单机无穷大系统发电机采用三阶实用模型。仿真参数如下发电机额定容量100 MVA惯性时间常数H4 s阻尼系数D2 p.u.暂态电抗Xd0.3 p.u.系统等效电抗XdΣ0.9 p.u.无穷大母线电压V∞1.0 p.u.励磁时间常数Tdo6 s机械功率Pm0.8 p.u.采样周期Δt0.01 s仿真总时长10 s过程噪声协方差矩阵Q设为diag([1e-6, 1e-5, 1e-6])对应三个状态量的不确定性程度。量测噪声协方差R由设定的信噪比决定这里取量测误差标准差约为量测真值的0.5%。初始状态x0取[0.5 rad, 1.0 p.u., 1.1 p.u.]初始协方差P0取diag([0.1^2, 0.01^2, 0.05^2])。初始协方差的选取值得多说一句。P0不能取太小如果P0过小滤波器会过度相信初始状态导致收敛缓慢P0过大又容易让滤波器在初始阶段出现较大波动。经验做法是根据系统初值的可信度来确定数量级一般比状态真值的方差略大一个量级就是安全的。我试过P0取diag([1, 1, 1])时滤波器仍然能收敛只是前几步的估计轨迹更粗糙一些。3.2 EKF核心代码逐行解读完整代码在文末提供这里先贴出EKF核心循环部分逐段说明。% 状态转移函数 f f_func (x) [ x(1) dt * (x(2) - omega_s); x(2) dt * (Pm - Pe_calc(x(3), x(1)) - D*(x(2) - omega_s)) / (2*H); x(3) dt * (Efd - x(3) - (x(3) - Xd_prime * Id_calc(x(3), x(1))) ) / Tdo_prime; ];注意Pe_calc和Id_calc是从当前状态x(3)和x(1)算出发电机电磁功率和d轴电流的非线性函数它们是整个系统非线性的来源。在代码中我单独写成子函数了方便单独调试。雅可比矩阵F的计算这里用符号工具箱推导后直接写解析解% 状态转移矩阵 F ∂f/∂x 的解析表达式 F [1, dt, 0; dt * dPe_ddelta, 1 - dt*D/(2*H), dt * dPe_dEq / (2*H); 0, 0, 1 - dt*(1 Xd_prime * dId_dEq) / Tdo_prime];由Eq对功角δ的偏导、Eq对E′q的偏导都涉及三角函数和电压方程的链式法则推导过程相当繁琐但一旦写对就是固定的代码块后期维护成本极低。建议大家在写自己的模型时也优先用解析雅可比矩阵而不是用数值微分比如有限差分法。数值微分每次滤波都要额外多次计算f函数精度还受步长选择的制约。量测矩阵H是输出有功Pe和机端电压Vt对状态量的偏导这里ds量测选取的是[Pe; Vt; δ]三组量测所以H是3×3矩阵。对于Pe的偏导项dPe/dx用链式法则展开后填入即可。修正步代码% 计算卡尔曼增益 S H * P_pred * H R; K P_pred * H / S; % 更新状态和协方差 x_est x_pred K * (z_meas - z_pred); P_est (eye(3) - K * H) * P_pred;这里z_pred h(x_pred)是预测量测值z_meas是含噪声的量测真值由真值加上高斯噪声生成。整个EKF的核心步骤写下来不到20行但它背后是对非线性系统、雅可比矩阵和贝叶斯滤波框架的完整理解。3.3 UKF核心代码实现与参数调优UKF的实现比EKF直观一些不需要求导即可完全依赖函数求值。核心步骤如下% 生成Sigma点 n 3; alpha 1e-3; kappa 0; beta 2; lambda alpha^2 * (n kappa) - n; [r, c] size(P); sqrtP chol((n lambda) * P, lower); Sigma zeros(n, 2*n1); Sigma(:,1) x; for i 1:n Sigma(:, i1) x sqrtP(:, i); Sigma(:, ni1) x - sqrtP(:, i); end % 计算权重 Wm [lambda/(nlambda), repmat(1/(2*(nlambda)), 1, 2*n)]; Wc Wm; Wc(1) Wm(1) (1 - alpha^2 beta);这里用的是Cholesky分解来求矩阵平方根所以要求协方差矩阵必须正定。在仿真过程中如果发现P矩阵变得病态一个常用技巧是加一个小的对角阵做正则化比如P 1e-12*eye(n)。实际运行中我发现当过程噪声Q设置过小时P矩阵会出现快速收缩到接近奇异的情况这时候Sigma点会非常集中滤波效果反而变差所以Q的选取也是保证UKF正常工作的关键。Sigma点通过状态方程传播后Sigma_pred zeros(n, 2*n1); for i 1:2*n1 Sigma_pred(:,i) f_func(Sigma(:,i)); end x_pred sum(Wm .* Sigma_pred, 2); P_pred Q; for i 1:2*n1 diff Sigma_pred(:,i) - x_pred; P_pred P_pred Wc(i) * (diff * diff); end量测更新类似把Sigma_pred再经过h函数传播得到Sigma_z然后计算预测量测均值z_pred、新息协方差S、互协方差Pxz最后用互协方差构造卡尔曼增益K Pxz / S更新状态和协方差。UKF调参这块我踩过不少坑。alpha选得太小Sigma点离均值太近导致量测更新对噪声不敏感alpha选得太大比如0.1Sigma点过于分散非线性传播的精度下降滤波轨迹可能出现尖刺。在我的仿真工况下alpha取1e-3表现最优beta取2高斯分布最优kappa取0即可。该组参数的敏感性分析输出在文末图中可以直观看到参数变化对RMSE的影响。4. 仿真结果分析与性能对比4.1 状态估计轨迹与收敛性对比为了公平起见EKF和UKF使用同样的系统模型、同样的初始值、同样的噪声序列Matlab中用rng固定随机种子。下图读者运行代码后得到展示了功角δ、角速度偏差ω和暂态电动势E′q三个状态的估计结果。整体看EKF和UKF都能跟踪系统状态的真实变化但在关键转折处存在差异。在功角的估计中EKF在校正阶段会出现较明显的滞后尤其在扰动发生后前几个采样点估计值与真值之间存在一个“延迟型”偏差。这是EKF线性化误差的典型表现——当状态变化率较大时一阶近似不足修正力度不够。UKF的响应明显更快因为Sigma点能更好地抓住状态的非线性传播特征。从数值上看运行代码得到的稳态RMSE对比大致如下以具体运行结果为准状态量EKF稳态RMSEUKF稳态RMSEδ (rad)0.00870.0042ω (p.u.)0.00150.0008E′q (p.u.)0.00630.0031UKF的精度整体优于EKF约50%这和理论预期一致UKF至少二阶精度EKF只有一阶。值得注意的是角速度偏差ω的估计噪声相对最低因为ω的动力学方程线性度最高两种滤波的表现差距最小。4.2 噪声统计特性对滤波性能的影响在状态估计中噪声统计特性的设定是“看不见但无处不在”的影响因素。如果实际噪声特性和滤波器假设的不匹配再好的算法也会退化成“看起来在滤波、实际上在发散”。我做了三组实验第一组Q和R完全匹配仿真噪声第二组R放大10倍模拟量测设备精度实际较差的情况第三组Q放大50倍模拟模型不确定性被严重低估的情况。结果发现在R失配时EKF和UKF都会变钝——增益变小估计轨迹变得平滑但滞后明显在Q失配时滤波器会出现过度信任量测、状态轨迹抖动加剧的现象。UKF在Q失配情况下仍然能保持稳定而EKF在强非线性工况下出现了轻微震荡。这说明UKF对参数失配的鲁棒性更好工程应用中这是一个值得考虑的加分项。对于实际工程中如何确定Q和R我的经验是先根据物理模型和量测仪表的标称精度确定初始值然后利用新息序列的统计特性做在线自适应。比如通过窗口化计算新息的协方差再反向调节Q和R这在移动窗口自适应卡尔曼滤波中是很成熟的做法。4.3 计算效率对比与实时性评估仿真运行在普通笔记本上Intel i5处理器、16GB内存、Matlab R2022b单步仿真耗时统计如下算法单步平均耗时msEKF0.18UKF0.42这里UKF因为要传播7个Sigma点2n17每个点都要计算一次非线性的f和h计算量大概是EKF的2.3倍。即便这样在0.01秒的采样周期下UKF的单步耗时仍然远小于实时约束完全可以跑在实时仿真系统里。如果你的系统规模扩大比如多机系统状态量达到几十维UKF的计算量会迅速上升。Sigma点数量是2n1意味着每增加一个状态量要额外传播两个点。这时可以考虑稀疏化Sigma点策略或者用无迹变换的降阶形式。对于常规电力系统动态状态估计场景10机39节点系统约60维状态UKF仍然在可接受范围内。5. 常见问题与排查技巧5.1 滤波器发散的五种典型原因滤波发散是状态估计中最让人头疼的问题。我总结过五种典型原因基本能覆盖90%以上的发散场景初始协方差P0设置不合理。P0过小导致滤波器对初始状态过度自信如果初始值又偏差较大滤波器会一直“坚持”错误状态无法纠正。解决办法是把P0放大1-2个量级让滤波器有足够的自由度去修正。过程噪声Q设置过小。Q描述的是系统模型的不确定性如果设置过小滤波器认为自己的状态方程很准导致增益很小、不信任新量测。在电力系统模型中负荷波动和模型简化带来的不确定性往往比我们预期的更大。建议先让Q比理论估计值大一个量级试试。雅可比矩阵求导出错。EKF中对导数符号、公式的细节要求极高一个符号错误导致的发散是“悄无声息”的——前面几步正常后面越来越偏。我在代码里用了符号工具箱和手工推导双重校验配合使用有限差分法交叉验证关键点确保F和H矩阵正确。量测方程和状态方程不一致。比如代码中f函数里计算Pe用的公式和量测方程h里计算Pe用的公式不是同一套会导致滤波器内部自洽性破坏滤波结果自然不对。请确保两个函数中的公共物理量如Pe、Id调用同一个子函数避免重复定义。Sigma点生成的Cholesky分解失败。协方差矩阵非正定时chol会报错或返回错误结果。常见原因是数值下溢或者P在迭代中失去对称性。解决方法是每次更新后强制对称化P (P P)/2并且在chol前加一个小的对角正则项。5.2 量测配置不同时EKF和UKF的表现差异我额外做了一组实验只保留一组量测仅Pe可测比较EKF和UKF的可观测性表现。结果很有意思单量测情况下EKF对ω的估计精度下降明显RMSE比三量测工况高了近一个量级而UKF仍然能保持相对稳定的估计只是收敛速度变慢了。这说明UKF在系统可观测性弱化时的鲁棒性更强因为Sigma点传播实际上利用了更多的状态分布信息即使在欠定测量下也能维持基本的估计能力。工程中的启示是如果你的系统只有SCADA量测、PMU覆盖不完整UKF可能是更稳妥的选择。当然严格的可观测性分析仍然是必要的——如果系统本身不可观测再强的滤波器也无法“无中生有”此时需要调整量测布置方案。5.3 调参心得快速定位最优参数的粗调-微调两步法参数调节不需要盲目穷举。我的经验是先粗调后微调。粗调阶段alpha取1e-3kappa取0beta取2P0根据初值可信度设一个量级Q先设成状态量典型幅值的1%的平方。跑一次仿真观察估计轨迹是否跟随真值变化。如果轨迹基本跟随但噪声大说明R偏小或者Q偏大把Q缩小一个量级再试。如果轨迹滞后明显、过度平滑说明R偏大或Q偏小把R缩小、Q放大再试。微调阶段在粗调结果附近细调Q的对角元关注RMSE指标的变化梯度。通常3-5次微调就能找到性能平台区再往下调收益不大、且对噪声失配更敏感。调参要避免一个误区不要追求RMSE绝对最小。工程应用中滤波器对噪声统计失配的鲁棒性往往比极端工况下的标称精度更重要。一个RMSE略大但参数变化时性能稳定的滤波器远好过一个RMSE极小但对参数极其敏感的滤波器。6. 扩展方向从单机系统到多机系统的迁移代码目前跑的是单机无穷大系统但框架是完全可扩展的。如果你想迁移到多机系统比如IEEE 39节点系统有几个关键点需要注意多机系统的状态维度大幅增加。以10机39节点系统为例每台发电机3个状态量总状态数达到30个。UKF的Sigma点会变成61个计算量指数级上升。EKF的优势在这个场景下会更明显——它的雅可比矩阵虽然是稀疏的但总计算量远低于UKF。实际工程中多机系统的在线动态状态估计常常选择EKF或者它的变体如平方根EKF因为实时性压力更大。另一个重要变化是量测方程更加复杂。多机系统中每台发电机的Pe不仅取决于自身的功角和E′q还跟其他发电机的状态量有关需要通过全网潮流方程耦合。这时候h函数的计算需要调用潮流计算或牛顿-拉夫逊法每步滤波都要做一次全网潮流计算负担很重。一些研究中采用灵敏度矩阵近似解耦把h简化成局部量对状态量的映射能大幅降低计算量但会损失一定精度。对于非线性程度更强的模型如发电机五阶模型、含励磁系统动态的完整模型UKF无需重新推导雅可比矩阵的优势会进一步放大。EKF需要推导更高维的雅可比矩阵推导出错的风险和调试成本都会上升。这也是我在新工科课题中更推荐用UKF作为基线算法的原因——你把模型改得再复杂UKF只需要调整f和h两个函数表达式不需要动滤波器主体逻辑。就我个人在实际操作中的体会来说EKF和UKF不是简单的谁替代谁的关系。如果你的系统模型简单、维度低、对实时性要求高而且你有耐心推导雅可比矩阵EKF完全可以胜任如果你面对的是高维强非线性系统或者你需要频繁调整模型结构UKF能让你从求导的泥潭里解脱出来把精力集中在系统建模和结果分析上。判断一个状态估计算法好不好不能只看RMSE还要看你对模型的掌控程度、调试成本和长期维护的便利性。这个项目把两者放在同一框架下对比最终目的就是让你在具体工程问题里能做出有依据的选择。