企业官网建设流程全解析

我最初接触电力系统动态状态估计（DSE）时，第一反应是拿EKF和UKF直接各跑一遍Matlab代码，比一比谁的误差更小。但真到动手写代码才发现，关键瓶颈根本不是滤波器本身，而是怎么把发电机动态模型和量测方程组织成可计算的形式。静态状态估计只需要解一个非线性最小二乘问题，动态状态估计却要同时处理状态预测和量测修正两条链路，模型稍微给错，EKF和UKF都会在几十个采样点内飘飘飘地发散掉。

这篇文章就围绕我用Matlab实现EKF和UKF做电力系统动态状态估计的完整过程展开。适合两类人：一类是刚接触动态状态估计、想要一份能跑通的代码框架的初学者；另一类是已经在用静态WLS（加权最小二乘状态估计）做断面分析，想往预测型估计器转的工程师。你不需要提前懂很多，但最好对卡尔曼滤波有一点点印象，否则建议先把线性卡尔曼的五条公式过一遍再来读。

1. 为什么静态状态估计满足不了动态分析的需求

1.1 静态估计的"冻结断面"假设错在哪里

传统电力系统状态估计是典型的静态问题：给定某一时刻的全部量测，求电压幅值和相角的最优解。SCADA系统大概每2到5秒推送一次遥测数据，调度端拿到这批数据后做一次WLS估计，输出一个断面。这个逻辑在运行方式变化缓慢的年代没有问题，但新能源出力波动、负荷快速变化、直流输电功率调整这些场景出现后，断面的间隔期里系统照样在演化，量测到了下一帧却还没传上来，调度员看到的状态就已经滞后了。

更麻烦的是SCADA的刷新率并不均匀，有的厂站遥测晚到几百毫秒甚至几秒，状态估计器只能等全部数据齐了再算。而动态状态估计的思路完全相反：用上一时刻的状态和系统模型去预测当前时刻的状态，再用最新量测修正预测值，量测不需要等齐，来一个修正一次，天然适合处理数据延迟和丢失。

1.2 动态模型从哪里来

做动态状态估计必须先有一个状态方程。电力系统里最经典的动态模型是发电机的摇摆方程，它描述了转子转速与电磁功率之间的平衡关系。如果手头有发电机的详细参数，状态量可以取功角、转速、暂态电动势；如果只做母线级的估计，文献里也有一种更务实的做法——把状态限定为节点电压幅值与相角，状态转移用随机游走或一阶惯性方程近似。

我在Matlab实现里采用的是后者，状态向量直接取除参考母线外的所有节点电压幅值和相角。这样做有两个原因：一是通用性强，换一个系统拓扑不需要重新推导发电机模型；二是量测方程直接复用潮流计算的映射关系，省去了从发电机内电势到网侧状态量的复杂转换。代价是状态转移模型的物理意义弱一些，Q矩阵需要认真调，这部分留到后面的调参章节细讲。

1.3 量测方程与状态方程的分工

动态状态估计能成立，靠的是状态方程和量测方程各管一段。状态方程负责描述系统的连续演化规律，量测方程负责描述量测设备和状态之间的映射关系。在母线级模型里，量测方程就是节点注入功率和线路潮流的非线性表达式，和静态估计完全一样。

这两条链路的配合很有意思：状态方程预测不准，量测修正还能往回拉；量测数不够，状态预测也能顶一阵子。冗余度低的可观测系统在静态估计里很难收敛，动态估计反而能利用模型信息维持可观性，这一点在实际工程里很有价值。

2. EKF方案：在非线性系统里做卡尔曼的代价

2.1 从线性卡尔曼到扩展卡尔曼

卡尔曼滤波本身是线性系统的产物，五条公式里最核心的一步是误差协方差的传播：P(k+1|k)等于A乘以P乘以A转置再加Q。一旦状态方程或量测方程变成非线性，这套传播公式就不能直接用。EKF的做法很简单也很暴力：在工作点附近对非线性函数做一阶泰勒展开，求出雅可比矩阵，然后把雅可比当成线性卡尔曼里的系统矩阵A和量测矩阵H。

电力系统的量测方程是潮流方程，非线性程度中等偏强，EKF的性能很大程度上取决于线性化工作点选得准不准。状态预测一步完成之后，量测更新之前的预测状态就是线性化点，预测误差越小，雅可比越准确，滤波效果越好。这正是EKF的软肋：遇到强突变工况，预测误差大，线性化点偏离真实状态远，雅可比矩阵失真，然后在下一拍继续劣化。

2.2 电力系统里的雅可比矩阵构造

以母线注入有功为例，极坐标下第i条母线的有功注入可以写成：

P_i = V_i * Σ_j V_j (G_ij * cos(θ_i - θ_j) + B_ij * sin(θ_i - θ_j))

对状态向量里的每个电压幅值和相角求偏导，才会得到量测矩阵H。我在代码里用符号工具箱先推导解析表达式，再转成函数句柄。一开始偷懒试过数值差分求雅可比，跑IEEE 39节点系统时误差不大，但换到重负荷算例后协方差矩阵经常出现非正定，后来全部改成解析雅可比才稳定下来。

EKF的量测更新公式没有变化，仍然是卡尔曼增益K、状态修正和协方差更新三条主线，只不过H矩阵变成了当前工作点下的雅可比。需要注意的坑是：雅可比矩阵每一列都要对应状态向量里的一个元素，如果状态里混入了常数项或者参考母线的相角，维度对不上，Matlab里矩阵乘法的报错会排山倒海一样来。

2.3 EKF的已知短板

EKF最大的问题是一阶线性化在大扰动场景下不够用。振荡、甩负荷这些工况会让状态量跑出线性化区间，雅可比矩阵的局部近似失效。还有一个工程上很痛的点：雅可比矩阵的解析表达式和系统拓扑强相关，改一次接线方式，所有偏导公式都要连带修改，维护成本极高。

但我依然建议初学者先跑通EKF再上UKF。原因很简单，EKF的每一行代码都能和卡尔曼五条公式对应上，出错时容易排查。UKF虽然数学上更优雅，调试起来反而更黑盒，协方差出现问题不好定位。

3. UKF方案：绕开雅可比矩阵的确定性采样

3.1 sigma点传播的核心思想

UKF的思路和EKF完全不同，它不去求导，而是用一组精心挑选的sigma点来逼近非线性函数的统计特性。假设状态维度是n，UKF会生成2n+1个采样点，这些点经过非线性函数传播之后，用加权均值和加权协方差来近似真实的后验分布。在线性卡尔曼里，协方差传播严格精确；在EKF里是一阶近似；在UKF里则是"无迹变换"的近似，精度至少达到泰勒展开的二阶项，对强非线性系统的适应性明显更好。

sigma点的生成逻辑是：在状态均值周围按照协方差矩阵的Cholesky分解结果做偏移，偏移量由尺度参数决定。Matlab里直接调chol函数就能得到下三角矩阵，但要注意如果P矩阵因为数值误差变得非正定，chol会报错，所以必须在每次协方差更新后加一个对称化和正定性修正。

3.2 权重设计与参数选择

sigma点的均值权重和协方差权重由几个参数控制，通常取alpha=1e-3，beta=2，kappa=0或3-n。alpha决定sigma点离均值的距离，取太小可能让数值稳定性变差；beta在高斯分布假设下取2是最优的；kappa用于保证协方差矩阵的半正定性。这些参数不是金银细软随便调的，alpha取多少直接关系到高阶项误差的放大程度，我个人的经验是alpha落在1e-4到1e-2之间比较稳妥，再小就要注意数值精度问题。

UKF的不变性是一个容易被忽略的优点：由于不需要求导，换量测函数、换网络拓扑都只需要改函数本身，滤波器的骨架完全不用动。对于EKF必须重新推导雅可比矩阵的场景，UKF省掉的工程量非常大。

3.3 为什么UKF在电力系统里更稳

电力系统的量测函数是三角函数和乘积项的叠加，在重负荷或者电压偏低的时候非线性程度接近临界，EKF的一阶线性化会产生明显截断误差。UKF的sigma点传播等效于用多个采样点去拟合真实映射，即便工作点移动，精度也不会瞬间恶化。我在算例里观察到一个典型现象：电压幅值从0.95pu往下掉的过程中，EKF的误差开始周期性振荡，UKF还能平稳跟踪。

需要泼一点冷水：UKF不是银弹。sigma点数量是2n+1，状态维度一旦上百，每一拍的计算量会明显上涨。另外UKF对协方差的数值品质更敏感，尤其在角度类状态上，很小的舍入误差可能让权重出现负值，最后导致协方差非正定。所以做UKF一定要在更新后强制检查P矩阵。

4. Matlab代码实现：状态推进、量测更新与整体架构

4.1 数据组织方式

我把整套代码按模块拆分，最顶层是一个主脚本run_dse.m，负责加载系统参数、生成量测、循环调用滤波器的predict和update方法、最后画图存结果。系统参数统一放在一个结构体里，包括母线导纳矩阵Ybus、节点编号、量测位置索引，这样换算例只改数据文件，不需要动滤波逻辑。

滤波器用Matlab的类来封装，基类定义predict和update两个抽象接口，EKF类和UKF类分别继承实现。这个设计让我后来加平方根UKF和CKF变得非常轻松，只需要多写一个类文件，主循环一行不需要改。对初学者来说面向对象不是必需品，但如果你想反复换滤波器做对比研究，这个设计能省很多重复劳动。

4.2 量测数据生成

动态状态估计仿真需要先有"真值"，才能评价滤波器好坏。我先用潮流计算在若干时间断面上生成一组基准状态，再把基准状态代入量测方程得到无噪声量测，最后叠加上指定方差的高斯噪声。这里有一个关键技术细节：状态真值的时间演化必须由状态方程生成，而不是随便离散几条曲线。做法是先用状态方程从初始状态递推出真实状态序列，再把每拍状态映射到量测空间。否则滤波器里的状态方程和仿真里的真实动态不一致，评估结果毫无意义。

量测类型我采用PMU风格的电压相量加上少量注入功率，PMU采样率设定为每周期50帧，即采样间隔0.02秒。这个时间尺度能体现动态估计的价值，SCADA风格的秒级数据在这种算法演示里看不出太大区别。

4.3 滤波主循环与EKF/UKF的分叉点

每拍滤波循环中的逻辑可以浓缩成下面这段伪代码结构：

% 状态预测 x_pred = f(x_est, u); if strcmp(filterType, 'ekf') [A, ~] = stateJacobian(x_est, u); P_pred = A * P_est * A' + Q; else P_pred = sigmaPropagate(x_est, P_est, @f); end % 量测更新 if strcmp(filterType, 'ekf') H = measJacobian(x_pred); K = P_pred * H' / (H * P_pred * H' + R); else K = ukfUpdate(x_pred, P_pred, z, @h); end x_est = x_pred + K * (z - h(x_pred)); P_est = P_pred - K * (H * P_pred)'; % EKF P_est = (P_est + P_est') / 2; % 强制对称

EKF的分叉点主要在雅可比的计算和常规卡尔曼更新；UKF的分叉点在sigma点的生成和传播。核心量测更新公式的形态一致，这让比较两者时的变量控制很干净，除了滤波器机理不同，其他环节完全相同。

4.4 误差评估与绘图

我习惯同时保存每拍的估计误差、协方差迹以及卡尔曼增益的范数。单看误差曲线会误导人，因为某一次随机噪声的波动可能很大，一定要配合协方差迹观察滤波器是否"自信"得过火。绘图时用三张图：电压幅值跟踪曲线、相角跟踪曲线、RMSE随时间变化曲线。相角曲线里要注意角度归一化，不然从+179度变到-179度的那一拍会被误判成巨大误差。

5. 算例测试：IEEE 39节点系统上的EKF与UKF对比

5.1 测试场景设置

我在IEEE 39节点系统上做了测试，状态量取除参考母线外的38个相角加上39个电压幅值，共77维。量测配置为一部分母线配置PMU电压相量，再加若干注入有功和无功功率，保证全网可观测。状态方程的动态模型采用一阶惯性加随机游走的混合形式，时间常数按系统惯量大致整定。

噪声协方差R按PMU精度设定，电压幅值标准差0.001pu，相角标准差0.001弧度。过程噪声协方差Q我花了不少时间调，最后电压幅值对应元素取1e-6量级，相角对应元素取1e-5量级。Q太大会让状态预测完全不信任模型，滤波退化成逐拍静态估计；Q太小又会让滤波器过于相信模型，量测失灵时协方差收缩过快，后面专门讲。

5.2 实测曲线与误差行为

以小扰动场景为例：系统稳定运行到第200拍时，某条母线负荷阶跃增加5%。两个滤波器都跟住了状态变化，但细节差异明显。EKF在扰动发生后的第1至第3拍出现一个误差尖峰，随后收敛回来；UKF从第一拍起误差就控制在较小范围，几乎看不到尖峰。这符合理论预期：负荷阶跃瞬间状态轨迹快速弯曲，一阶线性化的局部假设短暂失效。

我还跑了一个更容易发散的场景，把某台机组出力按正弦曲线大幅波动，EKF产生了几次幅度较大的协方差收缩，P矩阵一度接近奇异，UKF则一直保持平滑。需要注意的是，这两个都是单次仿真结果，不同随机种子下的表现会有波动，应该用Monte Carlo多跑几十遍再下结论。

5.3 数值指标与时间开销

单次仿真结束以后，我统计了RMSE和单拍平均耗时。典型结果EKF的电压幅值RMSE在0.0031pu左右，UKF在0.0022pu左右；相角RMSE两个滤波器差距更明显，EKF约0.0028rad，UKF约0.0017rad。时间上EKF的单拍耗时要低一些，因为一次量测更新只需要计算一个雅可比矩阵，UKF要传播155个sigma点。

但这不代表EKF永远更快，解析雅可比推导需要人工成本，如果算上推导时间，UKF反而更划算。我在项目里已经把解析雅可比提前算好存成函数文件，所以单拍耗时才压得比较低。

5.4 什么时候EKF够用

如果系统运行平稳、量测冗余度高、没有频繁的大扰动，EKF的精度和UKF差距并不大，有时甚至因为数值更稳定反而表现更好。UKF的优势在强非线性或扰动频繁的场景下才能体现出来。所以我不建议盲目追求UKF，先把系统的动态特性想清楚再选。

6. 调参心得与滤波器发散避坑

6.1 Q矩阵和R矩阵怎么定

R矩阵相对好办，按量测设备的精度指标来就可以。Q矩阵是动态估计里最让人头疼的参数。它没有明确的物理标定方法，只能从模型误差的估计出发：你的状态方程离真实动态有多远，Q就设多大。模型越粗糙，Q应该越大，给量测修正留出余地。

我在代码里写了一个辅助函数，可以在一定范围内扫描Q的乘子系数，自动输出不同系数下的估计RMSE。这样调参不是盲目的，而是先看RMSE谷值在哪一段，再在附近做细化。对于77维状态，全矩阵扫描不可能，我把它拆成幅值子块和相角子块，分别给出统一的乘子。经验上电压幅值子块和相角子块的Q系数可以差一个数量级，混在一起调很难收敛。

提示：如果滤波曲线出现"量测修正过头"的表现——估计值高频抖动、误差方差忽大忽小，优先减小Q；如果出现"量测拉不动"的表现——真值突变后估计值缓慢爬行，优先增大Q。这两个症状是调Q的可靠风向标。

6.2 协方差矩阵的非正定问题

P矩阵非正定在EKF和UKF里都出现过，是动态估计调试期最常遇到的坑。症状是Matlab报错说chol输入矩阵必须正定，或者RMSE曲线突然爆表。我做了三件加固：第一，每次P更新后强制对称化，也就是(P+P')/2；第二，加一个小的对角修正项，类似于Levenberg-Marquardt的做法，让最小特征值保持在某个阈值之上；第三，定期检查P的特征值，如果出现显著负值，说明模型发散了，靠数值修补救不回来，必须回头调Q。

UKF对P矩阵的正定性要求比EKF更高，因为sigma点生成依赖Cholesky分解。我最开始在39节点系统上跑UKF发散，十次里有七八次都是这里出的问题。修正之后稳定了很多，但也要记住UKF发散的路径往往比EKF更剧烈，调参时要边跑边看协方差迹曲线。

6.3 采样周期与离散化的陷阱

动态状态估计的采样周期直接决定状态转移模型的形式。PMU的50帧每秒采样对应20毫秒间隔，在这个尺度上很多连续动态可以用简单的一阶离散公式近似；但如果数据源是SCADA的秒级刷新，状态方程里的时间常数就必须重新整定，否则一拍之内状态变化太大，预测根本不准。

还有角度处理的细节：相角是周期量，预测值和量测值之间做差一定要归一到[-pi, pi]区间，否则遇到180度附近的跳变会产生假的大误差，进而拉歪卡尔曼增益。这个坑在Matlab里尤其隐蔽，因为Matlab的sin和cos函数不会提醒你角度差异异常。

6.4 我做这套代码后的几条总结

滤波发散时先别动滤波器结构，先把Q、R两个矩阵和初始协方差P0过一遍，大部分问题都出在这三个参数上。P0设得太小会让滤波器一开始就过度自信，量测稍微不稳就发散；P0设太大又会让前几十拍的估计噪声偏大。我一般取对角元素0.1到1之间的对角阵，让滤波器先通过前几十拍把协方差降下来。

如果你也是从静态估计切到动态估计的，我的建议是先别急着上UKF。老老实实把EKF在简单系统上跑通，体会状态预测和量测修正各占多少分量，再把系统换到39节点，最后才换成UKF做对比。这样每一步出了偏差你都清楚该去查哪一段代码。这套Matlab框架本身可以直接拿去做二次开发，换数据文件、换量测配置、换滤波器类型都不需要重写主循环。我后来加扩展卡尔曼的强跟踪版本，也只用了一天时间，架构上省下来的功夫相当可观。