在写电力系统的状态估计算法时很多人一上来就扎进EKF和UKF的公式推导里结果往往卡在滤波器发散协方差矩阵非正定估计结果震荡这类问题上。这篇文章从我的实际仿真经验出发把扩展卡尔曼滤波和无迹卡尔曼滤波在电力系统动态状态估计里的建模、代码思路、参数整定以及踩坑过程完整过一遍所有代码段都基于Matlab实现可以直接参考改造。1. 动态状态估计解决的是什么问题从静态断面到动态追踪的关键跨越1.1 为什么传统静态估计在电力系统中不够用了电力系统的状态估计传统上指的都是静态状态估计也就是基于加权最小二乘法WLS对某一时刻的SCADA量测断面做一次快照式的求解得到当前电网各节点的电压幅值和相角。这个方法成熟、可靠工程上用了三四十年但它的假设前提是系统处于稳态——量测断面之间的时间尺度足够长系统状态几乎没有变化。可实际运行中电网的动态过程是连续发生的负荷波动、机组调速器动作、励磁系统调节、故障后的暂态过渡这些都是秒级甚至毫秒级的动态过程。SCADA的刷新周期在2到10秒左右而且各量测点的上传时刻并不严格同步拿这些不同步的快照去做静态估计在系统波动较快的时段会引入明显的模型误差。PMU同步相量测量单元普及之后量测数据的时间分辨率提升到了毫秒级、微秒级但量测点覆盖率又不足以支撑全网静态可观测。这个时候动态状态估计——用系统的动态模型把状态量的演变规律和实时量测结合起来——就成了必然选择。动态状态估计和静态估计的本质区别在于静态估计只用了量测方程也就是状态量到量测值的映射关系动态估计多用了状态转移方程也就是系统状态随时间演变的动态规律。通俗点说静态估计是拿当前时刻的多个量测做一次最小二乘拟合动态估计是拿上一时刻的状态预测出当前时刻的状态再用当前时刻的量测把这个预测修正一下。卡尔曼滤波家族就是做这件事的标准框架。1.2 动态状态估计的对象从发电机到节点电压需要先厘清一个概念电力系统动态状态估计估计的对象不是整个电网的全部节点电压那是静态状态估计干的活而是发电机组的动态状态。这一点我从一开始就理解偏了导致前几版仿真完全跑偏。标准的发电机三阶实用模型状态量通常取δ发电机功角radω发电机转子角速度标幺值同步速时为1.0Eqq轴暂态电动势标幺值微分方程组是dδ/dt ω0 × (ω - 1)2H/ω0 × dω/dt Tm - Te - D(ω - 1)Td0 × dEq/dt Efd - Eq - (Xd - Xd)×Id其中ω0是同步转速H是惯性时间常数Tm是机械功率Te是电磁功率D是阻尼系数Td0是励磁绕组暂态时间常数Efd是励磁电压Xd和Xd分别是d轴同步电抗和暂态电抗Id是d轴电流分量。量测方程则是从状态量映射到PMU/SCADA可测的物理量比如机端电压幅值Vt、输出有功功率Pe、无功功率Qe、转子角速度ω等。这里有一个非常关键的难点量测方程几乎全都是非线性的。Vt的求解需要用到Eq、δ以及网络方程Pe、Qe更是δ和Eq的三角函数组合。这就是EKF和UKF发挥作用的前提条件——线性卡尔曼滤波解决不了这类问题。2. EKF在电力系统状态估计里的实现逻辑线性化的代价与收益2.1 EKF的经典三步流程扩展卡尔曼滤波的思路非常直接既然状态转移和量测方程都是非线性的那就把它们在当前估计点附近做泰勒展开只保留一阶项把非线性系统近似成线性系统然后套用标准卡尔曼滤波的框架。这个思路在工程上应用极广因为它实现简单、计算量小几乎所有控制系统教材都会讲到。EKF在电力系统动态状态估计中的每一步是这样的预测步——先用状态转移方程算出先验状态估计x_pre f(x_est, u) wP_pre F × P_est × F Q这里的f就是上一节那个发电机微分方程组的离散化形式u是外部输入机械功率、励磁电压等w是过程噪声F是f对状态向量x的雅可比矩阵P是协方差矩阵Q是过程噪声协方差矩阵。更新步——用当前时刻的量测修正预测值K P_pre × H × (H × P_pre × H R)⁻¹x_upd x_pre K × (z - h(x_pre))P_upd (I - K × H) × P_pre其中H是量测函数h对状态向量x的雅可比矩阵z是当前量测向量R是量测噪声协方差矩阵。2.2 雅可比矩阵EKF最脆弱也最关键的部分EKF的整个精度几乎都押在了F和H这两个雅可比矩阵上。一阶线性化意味着系统的高阶信息被直接丢弃当系统的非线性程度较强、或者采样步长较大时这种丢弃会带来不可忽视的截断误差。我在Matlab里实现时雅可比矩阵的标准求法是解析求导。比如量测函数h(x)包含机端电压Vt sqrt((R×Id Xq×Iq Eq)² ...)这种表达式手推导数非常痛苦一不小心就推错。后来我学到一个更工程化的做法用复步长微分或者中心差分法做数值雅可比虽然计算量增加一些但至少避免了手推导数的低级错误。不过数值雅可比的精度非常依赖步长的选取。中心差分法里扰动步长取经验值h ≈ sqrt(eps) 量级eps是机器精度约2.2e-16在标幺值系统里我实测发现取1e-6到1e-8之间效果最稳定。步长太大则导数的有限差分误差明显步长太小则数值舍入误差占主导。EKF在电力系统里的实际表现很大程度上也取决于系统是否正处在强非线性区间。系统稳态时功角δ变化平缓量测函数接近线性EKF的估计精度和收敛性都很好但系统遭遇扰动、功角大范围摆动时三角函数的高阶项开始起作用EKF的一阶线性化近似就有些吃力了。我在后面第5章会给出具体的对比结果。3. UKF如何绕过线性化Sigma点策略与核心代码3.1 无迹变换的直觉理解UKF的全称是Unscented Kalman Filter中文一般译作无迹卡尔曼滤波。它的核心思想不是去逼近非线性函数而是去逼近状态量的概率分布。既然状态量的统计特性均值和协方差可以通过一组精心选取的采样点来传播那就不需要对非线性函数本身做任何线性化假设。这组采样点叫Sigma点。假设状态向量x是n维的均值为x_mean协方差为P那么取2n1个Sigma点χ₀ x_meanχᵢ x_mean (sqrt((n λ)P))ᵢi 1, ..., nχ_{in} x_mean - (sqrt((n λ)P))ᵢi 1, ..., n这里sqrt((n λ)P)表示协方差矩阵的Cholesky分解因子λ是一个尺度参数。直觉上理解这2n1个点就像是在均值附近沿着主成分方向对称撒开的侦察兵它们穿过非线性函数之后再由加权统计还原新分布的一阶矩和二阶矩。3.2 Sigma点的权重设置与参数选择UKF的权重分为均值权重Wm和协方差权重Wc具体为Wm₀ λ / (n λ)Wc₀ λ / (n λ) (1 - α² β)Wmᵢ Wcᵢ 1 / (2(n λ))i 1, ..., 2nλ α²(n κ) - n这里的参数选择直接关系到滤波的稳定性α控制Sigma点在均值周围的散布范围通常取1e-3到1之间的正数。α越大Sigma点离均值越远对非线性传播的捕捉能力增强但协方差更容易出现非正定α过小则会丢失高阶信息。我做电力系统仿真时通常取0.1到0.5之间既不保守也不激进。κ是一个次要尺度参数高斯分布背景下一般取0有些建议取3 - n。β与先验分布有关对高斯分布β 2是最优的这个值是文献里推导出来的直接写上即可。Sigma点生成时那个sqrt((n λ)P)在Matlab里用chol((n λ)P)来计算。需要注意这里的P必须是严格正定的否则chol会报错。我在调试初期遇到最多的就是这个错误——协方差矩阵迭代几轮之后出现负特征值导致chol崩溃。后面的排查链路里会详细讲。3.3 UKF的量测更新逻辑UKF的量测更新比EKF更简洁因为它不需要计算量测函数h的雅可比H。算法流程如下用状态转移函数f把每个Sigma点传播一步得到一组预测Sigma点对这组预测Sigma点加权求均值和协方差得到先验估计x_pre和P_pre把每个预测Sigma点代入量测函数h得到量测预测点集再加权算出量测预测均值z_pre、量测预测协方差S、以及互协方差C卡尔曼增益K C × S⁻¹最后更新状态和协方差。整个过程完全避开了雅可比矩阵非线性函数在Sigma点的传播过程中被原封不动地使用。这也是UKF相对EKF最大的优势不需要推导任何导数而且精度理论上可以达到三阶泰勒展开的水平而EKF只有一阶。4. Matlab代码实现一个可复现的仿真框架4.1 仿真系统与量测设置为了兼顾代码可读性和实现效果我采用一个单机无穷大系统SMIB作为仿真对象。这台发电机用三阶模型描述状态量是[δ; ω; Eq]量测量选四个发电机输出有功功率Pe、机端电压幅值Vt、转子角速度ω和功角δ。PMU可以直接测量功角和角速度有功功率和机端电压也都能以较高精度获取这个配置跟实际PMU量测场景是一致的。仿真步长取Δt 0.01s模拟时长为10s。过程噪声协方差Q取diag([1e-6; 1e-6; 1e-6])量测噪声协方差R取diag([1e-4; 1e-4; 1e-6; 1e-6])。这个取值后面会详细讲依据。核心代码如下%% 发电机三阶模型的状态转移函数 function x_next gen_smib_model(x, u, dt) % 状态量 x [delta; omega; Eq] % 输入 u [Pm; Efd] % 发电机参数 omega0 2*pi*50; % 同步转速(rad/s) H 4.0; % 惯性时间常数(s) D 2.0; % 阻尼系数 Td0 6.0; % 励磁绕组暂态时间常数(s) Xd 1.8; % d轴同步电抗(pu) Xd_prime 0.3; % d轴暂态电抗(pu) Xq 1.7; % q轴同步电抗(pu) Xq_prime 0.55; % q轴暂态电抗(pu) Rs 0.0; % 定子电阻(pu) delta x(1); omega x(2); Eq_prime x(3); Pm u(1); Efd u(2); % 网络方程单机无穷大系统 Vinf 1.0; % 无穷大母线电压 XdT Xd_prime Xd_prime; % 简化后的等值电抗 XqT Xq_prime Xq_prime; Id (Eq_prime - Vinf*cos(delta)) / XdT; Iq Vinf * sin(delta) / XqT; % 电磁功率 Pe Eq_prime * Iq (Xq_prime - Xd_prime) * Iq * Id; Vt sqrt((Vinf XqT*Iq)^2 (XdT*Id)^2); % 简化机端电压计算 % 微分方程组 ddelta omega0 * (omega - 1); domega (omega0/(2*H)) * (Pm - Pe - D*(omega - 1)); dE (Efd - Eq_prime - (Xd - Xd_prime)*Id) / Td0; x_dot [ddelta; domega; dE]; x_next x dt * x_dot; end %% 量测方程 function z gen_smib_h(x, u) delta x(1); omega x(2); Eq_prime x(3); Vinf 1.0; XdT 0.6; XqT 1.1; Id (Eq_prime - Vinf*cos(delta)) / XdT; Iq Vinf * sin(delta) / XqT; Pe Eq_prime * Iq (XqT - XdT) * Iq * Id; Vt sqrt((Vinf XqT*Iq)^2 (XdT*Id)^2); z [Pe; Vt; omega; delta]; end4.2 EKF的主循环核心代码EKF主循环里最需要谨慎的就是雅可比矩阵的处理。我在下面的代码里用的是解析式在工程里可以先用数值差分做交叉验证。%% EKF 主循环片断 for k 1:N % --- 预测步 --- x_pre gen_smib_model(x_est(:,k-1), u(:,k-1), dt); % 状态转移雅可比 F解析式此处示意 F eval_F_jacobian(x_est(:,k-1), u(:,k-1), dt); P_pre F * P_est * F Q; % --- 更新步 --- z_pre gen_smib_h(x_pre, u(:,k)); H eval_H_jacobian(x_pre, u(:,k)); K P_pre * H / (H * P_pre * H R); x_est(:,k) x_pre K * (z_meas(:,k) - z_pre); P_est (eye(3) - K * H) * P_pre; end这里有个小细节值得注意Matlab里求K我用的是/运算符也就是右除等价于K P_pre * H * inv(H * P_pre * H R)但数值上更稳定也更快。建议不要在代码里写inv尤其是矩阵维度增大之后inv的性能和稳定性都不理想。4.3 UKF的主循环核心代码UKF的核心是Sigma点的生成与传播我用一个单独的脚本把UT变换写清楚。%% UKF 主循环片断 n 3; % 状态维度 alpha 0.3; kappa 0; beta 2; lambda alpha^2 * (n kappa) - n; % 生成Sigma点 P_sqrt chol((n lambda) * P_est, lower); X_sig zeros(n, 2*n1); X_sig(:,1) x_est(:,k-1); for i 1:n X_sig(:, i1) x_est(:,k-1) P_sqrt(:,i); X_sig(:, ni1) x_est(:,k-1) - P_sqrt(:,i); end % 权重 Wm zeros(1, 2*n1); Wc zeros(1, 2*n1); Wm(1) lambda / (n lambda); Wc(1) lambda / (n lambda) (1 - alpha^2 beta); Wm(2:end) 1 / (2*(n lambda)); Wc(2:end) 1 / (2*(n lambda)); % 预测步Sigma点通过状态转移函数 X_pre zeros(n, 2*n1); for i 1:2*n1 X_pre(:,i) gen_smib_model(X_sig(:,i), u(:,k-1), dt); end x_pre X_pre * Wm; P_pre (X_pre - x_pre) * diag(Wc) * (X_pre - x_pre) Q; % 更新步Sigma点通过量测函数 Z_pre zeros(4, 2*n1); for i 1:2*n1 Z_pre(:,i) gen_smib_h(X_pre(:,i), u(:,k)); end z_pre Z_pre * Wm; S (Z_pre - z_pre) * diag(Wc) * (Z_pre - z_pre) R; C (X_pre - x_pre) * diag(Wc) * (Z_pre - z_pre); K C / S; x_est(:,k) x_pre K * (z_meas(:,k) - z_pre); P_est P_pre - K * S * K;4.4 初始协方差矩阵与噪声矩阵的整定在仿真启动阶段我踩过最深的坑就是P0和Q、R三个矩阵的取值完全凭感觉导致滤波器在第2到3个采样周期后剧烈发散。后来梳理出一套可操作的经验**P0初始状态协方差**反映对初始状态估计的不确定度。如果初始状态是从静态状态估计或潮流计算得到的可以取较小的对角值比如1e-4到1e-2。如果完全不知道初值就取大一些让滤波器有足够的自由度去收敛。我P0取diag([1e-3; 1e-3; 1e-3])。**Q过程噪声协方差**反映的是模型误差也就是发电机微分方程组没有完全刻画的实际扰动。Q太小会让滤波器过度相信模型量测的修正作用变得微弱Q太大则会让状态轨迹跟着量测噪声乱跳。Q的值本质上是模型不准确程度的高斯近似标幺值系统下取1e-6到1e-4之间是我实测比较稳的范围。**R量测噪声协方差**反映的是PMU和互感器测量链路的误差水平。PMU的相量测量误差在稳态时很小总相量误差TVE通常在1%以内对应到标幺值约1e-4量级。但要注意R完全取实际噪声方差有时候会让卡尔曼增益过大在模型误差存在时引发震荡。工程上往往把R稍微放大一点相当于给量测也留出一定冗余。这张表总结了我调试过程中整理的经验区间矩阵物理含义典型取值范围标幺值调参方向P0初始状态不确定度1e-4 ~ 1e-2初始误差小就取小完全无先验就取大Q模型不确定度/过程扰动1e-6 ~ 1e-4轨迹太滑就增大轨迹太碎就减小R量测误差水平1e-6 ~ 1e-4量测硬件好就取小量测噪声大就取大比例关系比绝对数值更重要。我反复试过的经验是Q与R的比值决定滤波器在信任模型和信任量测之间的权重。调参不能只看单个矩阵要把Q和R放在一起感受系统的动态响应。5. 实测对比收敛性、数值稳定性与常见问题排查5.1 两种滤波器的估计精度对比在同样一套仿真系统、同样的初值和噪声设置下EKF和UKF的估计结果差异在稳态区段并不大两者的误差带都在可接受范围内。但在系统发生阶跃扰动比如机械功率Pm在t4s时突增5%之后的过渡过程里差异变得明显EKF估计的功角曲线出现了大约0.03到0.08 rad的瞬时偏差且需要约0.3到0.5s才回到真值附近UKF的偏差明显更小基本在0.02 rad以内回稳速度也更快。这个结果符合理论预期。扰动发生后功角δ和Eq的变化速率较快状态轨迹在短时间内偏离线性化工作点EKF用一阶雅可比外推的误差被放大。UKF用Sigma点传递非线性信息相当于保留了更多的非线性特征因此在暂态过程中的优势更加突出。但UKF不是没有缺点。同样维度下UKF单步计算量约为EKF的2到3倍因为每个采样周期要传播2n1个Sigma点每个点都要完整调用一次状态转移函数和量测函数。在实际工程里如果系统中需要估计的发电机组数量较多、量测通道庞大UKF的计算负担是必须考虑的。好在现代电力系统动态状态估计通常是分站部署、单机或机组群独立滤波计算量问题在很多场景下可以接受。5.2 数值不稳定的完整排查链路仿真过程中最让人崩溃的就是滤波器莫名其妙发散。我把调试中遇到的现象、排查步骤和最终解法整理成一条完整的链路供后来者对照排查。第一步看协方差矩阵是否非正定。在Matlab里直接用chol分解Sigma点因子如果报Matrix must be positive definiteP矩阵已经有非正特征值了。此时先检查每一步的(P_pre - KSK)是否可能出现浮点舍入导致的轻微不对称。解决方案是强制对称化P_est (P_est P_est) / 2。第二步检查雅可比是否出错。EKF中解析雅可比与实际系统不匹配是常见问题。我用的验证方法是拿数值雅可比做对照用中心差分算出的H和解析H比较如果最大误差超过1e-4说明解析推导有问题。不要嫌麻烦这一步能省下大量debug时间。第三步检查量测函数是否使用了与真值生成相一致的模型版本。真值量的生成、量测函数、滤波器的量测预测这三处必须严格用同一个函数。我犯过低级错误仿真生成量测用的Pe公式里漏了一项交叉项滤波器里却用了完整的公式导致永远存在固定偏差。这类问题不通过代码审查很难发现。第四步检查单位是否统一。电力系统标幺制下功角δ的单位是rad转子速度ω用标幺值同步速为1.0但在有些论文里ω用的是偏差量Δω ω - 1。状态转移函数、量测函数、初值三者的定义必须完全一致我见过最隐蔽的问题是量测里对ω加了一个1.0的基准偏移而滤波器里没有对应的逆变换。第五步检查采样步长是否过大。EKF的一阶线性化精度和步长强相关。如果仿真步长超过0.02sEKF的雅可比近似误差就会显著上升。UKF对步长的容忍度稍好但也不能无限放大。如果无法减小步长考虑换用更高阶的离散化方法比如用四阶Runge-Kutta替代前向欧拉能明显改善精度。5.3 EKF与UKF的选型建议做了大量对比实验后我对这两种滤波器的选型形成了一个比较清晰的判断标准系统非线性程度弱、状态变化平缓、对计算速度敏感的场景比如稳态运行状态下的机组在线监测EKF完全够用。它的实现简单、参数少、运算速度快后期维护也容易。很多实际工程系统至今还在跑EKF不是因为落后而是因为够用。系统可能进入强非线性区间的场景比如故障暂态过程分析、重负荷下的电压稳定监测、机组启动或停机过程UKF的更优表现值得额外的计算开销。系统复杂度高、雅可比解析推导困难的场景UKF可以大幅减少建模工作量。不需要推导导数意味着模型换一套就得重写雅可比的问题被绕开了。我的实际建议是在Matlab里先同时实现两套滤波器共用量测函数和状态转移函数做一个如第4章那样的对比框架。前期调试成本不高换来的是对两种算法在具体系统中的差异有直观认知。后续若要工程部署再根据实际算力选择合适的版本。调试中还发现一个容易被忽视的细节UKF里chol分解前最好加一个特征值裁剪把P矩阵小于某个阈值的负特征值清零这能避免很多偶发性的崩溃。EKF的协方差更新如果改用Joseph形式也就是P_upd (I - KH) * P_pre * (I - KH) KRK比标准形式在数值上鲁棒得多。这些都是教科书不会写、但实际调仿真一定能用上的细节。最后补一句个人体会。这类滤波算法的实现难点从来不在卡尔曼增益公式本身而在状态方程和量测方程的建模准确度、噪声协方差整定和数值稳定性处理。这三块做扎实了EKF和UKF在大多数电力系统场景下都能稳定输出高质量的状态估计结果。我的建议是做仿真时不要只盯着估计误差曲线好不好看先把协方差矩阵的演变过程打印出来看看是否收敛、是否保持在合理范围内——这往往能提前暴露很多问题。