拓冰网站定制 
1. 从静态断面到动态递推电力系统动态状态估计到底解决了什么问题在电力系统里状态估计这个词绝大多数场景下指的是调度中心SCADA/EMS里的加权最小二乘静态估计把一个电网断面里的母线电压幅值、相角、支路潮流和注入功率作为量测量解一个线性或弱非线性方程组得到一个稳态快照。这套办法在正常运行方式下非常可靠调度员看的就是那张大屏上的潮流分布。但问题来了当系统真的发生扰动——线路跳闸、机组跳机、大负荷突增——功角和转速开始围绕稳态点来回摆动的那几秒到几十秒静态状态估计的假设直接被击穿。断面量测还没算完实际状态已经跑了十万八千里。真要等到稳态恢复再算对广域保护、紧急控制来说黄花菜都凉了。所以我第一次接触电力系统动态状态估计(DSE)时心里很清楚这不是静态估计的升级版而是完全换了一套数学工具把发电机本身的微分方程拿进来用滤波器的思路逐拍递推估计功角、转速、暂态电动势这类真正随动态过程变化的状态变量。支撑DSE落地的硬件基础是这些年在电网里铺开的PMU同步相量测量单元。PMU以10到60帧每秒的速度上报带时标的相量和频率时间分辨率远高于传统RTU遥测。有了高频量测才谈得上动态估计的闭环校正。进一步说DSE的核心算法就是我们的标题主角扩展卡尔曼滤波EKF以及它之后发展出的无迹卡尔曼滤波UKF。EKF的思路最直观把非线性方程在工作点附近做一阶泰勒展开用雅可比矩阵硬套线性卡尔曼滤波公式UKF则绕开求导用一组sigma点去传播概率分布对强非线性模型更宽容。这篇文章我会从一个可复现的Matlab仿真入手把动态状态估计这整套链路讲清楚怎么搭发电机三阶模型、怎么构造量测方程、EKF和UKF的代码到底怎么写、哪些参数会让滤波发散、怎么从单机算例扩展到多机系统。适合的读者包括正在做电力系统课程设计或毕业设计的学生也包括刚接手广域测量系统数据分析、需要在现有PMU数据上试算动态估计的研究人员和工程师。我会把我在调试中踩过的坑一并写出来你们可以照着抄作业也可以拿来当避坑清单。需要先说明的是这是基于我过往项目经验的通用实践如果你要结合真实电网的详细模型需要在这个框架上修改参数矩阵和量测配置但核心的滤波更新流程基本不变。2. EKF与UKF的原理拆解线性化与sigma点的取舍2.1 EKF把非线性问题局部线性化的工程妥协卡尔曼滤波原本是为线性系统设计的在线性高斯假设下它是最小方差意义下的最优估计器。但电力系统动态模型是明显的非线性常微分方程组比如功角运动方程里有sin(δ)励磁方程里还叠着电压和电流的乘积项根本不可能用线性Kalman那套严格推导。EKF做的妥协很直白在每个时刻的工作点附近对状态方程和量测方程求雅可比矩阵用一阶偏导把非线性问题强行压成线性问题再套标准卡尔曼递推。EKF的标准递推分五步我直接写出来。假设离散状态方程是 x_k f(x_{k-1}, u_{k-1}) w量测方程是 z_k h(x_k) v过程噪声协方差Q量测噪声协方差R预测x_pred f(x_prev, u)P_pred F * P_prev * F Q其中F是状态方程雅可比在x_prev处求值。更新卡尔曼增益K P_pred * H * (H * P_pred * H R)^{-1}其中H是量测方程雅可比在x_pred处求值滤波值x_est x_pred K * (z - h(x_pred))协方差更新P_est (I - K*H) * P_pred。EKF的优点非常明显实现简单计算量在众多非线性滤波算法里几乎最小一个状态维度为n的系统主要代价是生成两个n×n雅可比矩阵并做两次矩阵求逆。我在Matlab里跑单机系统的时候EKF单步耗时基本在毫秒级完全能满足PMU的实时性要求。但EKF的短板也很致命当系统非线性较强、或者工作点变化剧烈时一阶泰勒展开的截断误差会迅速累积。最典型的情况是功角在暂态过程中摆动幅度很大sin(δ)的高阶项不再能忽略EKF的协方差估计就开始失真严重时直接发散。还有一个工程上的痛点雅可比矩阵推导容易出错尤其在多机系统里状态方程涉及所有发电机的互作用项手推F矩阵推错一两个符号滤波结果在视觉上看着还行但协方差早就病态了。2.2 UKF不用求导用sigma点传播统计量UKF的设计动机非常聪明既然直接求非线性函数的导数容易出错且只保到一阶那干脆用一组精心挑选的采样点sigma点去感受非线性函数的形状这些点在原分布的均值附近按协方差信息展开经过非线性函数传播后从它们的新位置重新计算均值和协方差。这样一来完全不需要推导雅可比矩阵而且统计精度理论上能捕获到非线性函数的二阶项比EKF的一阶局部分解更贴合真实分布。生成sigma点的标准做法是对于n维状态一共取2n1个点。设状态均值为x_m协方差为P引入比例因子λ α²(n κ) - n其中α通常取0.01到1之间的小正数κ一般取0那么sigma点构成为第0个点x_0 x_m第i个点i1,...,nx_i x_m sqrt((nλ)P)的第i列第in个点x_i x_m - sqrt((nλ)P)的第i列。权重系数方面均值权重为W_m0 λ/(nλ)其余点权重为W_mi 1/(2(nλ))协方差权重的第0项额外叠加一个(1-α²β)修正项β一般取2其余项与均值权重相同。把这组sigma点代入非线性状态函数和量测函数之后用加权平均重构预测均值和协方差再用与标准卡尔曼相同的增益公式做更新。UKF的实际计算量大概是EKF的2到3倍因为每个时刻要传播2n1组sigma点每组都是一次完整的非线性函数求值。但代价换来的是强非线性场景下明显更稳的结果而且代码里不需要手推任何雅可比矩阵。换句话说EKF的难点在建模时的推导UKF的难点在实现时的矩阵运算细节两者对工程人员的要求确实不同。2.3 两种滤波器的性能对照我用下面这张表总结两者差异这是我在实际仿真中最直观的感受也符合多数文献里的共识对比维度EKFUKF是否需要雅可比矩阵需要推导易错不需要零求导非线性处理精度一阶泰勒近似约二阶精度单步计算量较小适合高维状态约2~3倍EKF随维数线性增长强非线性暂态场景可能出现协方差失稳更稳残差更小实现复杂程度公式简洁但调试花时间代码略长但调参逻辑清晰适合场景弱非线性、模型维数高强非线性、模型故障切换频繁从我自己的调试经验看单机或小规模系统上UKF比EKF的优越性不一定体现在快或准上更多体现在你不需要把每一个偏导数都算对。省下的推导时间远比多跑的那几毫秒有价值。这也是我在文末会强烈建议初学者先跑通UKF的原因。3. 仿真模型与Matlab参数设计先把发电机三阶模型搭对3.1 单机无穷大母线系统建模基准做动态状态估计的仿真我不建议一上来就搭完整多机系统工程习惯是先跑通一个单机无穷大母线模型。它结构简单却保留了发电机动态的核心非线性特征转子摆动方程和励磁绕组的暂态过程。所谓无穷大母线就等价于一个幅值、频率恒定不变的大电网发电机挂在这上面功角摆动会直接反映电磁功率的变化正是我们想用滤波器跟踪的动态。我采用的经典三阶单轴发电机模型状态向量选为x [δ; Δω; Eq]三个状态分别对应δ发电机功角单位弧度。Δω转子转速相对同步转速的标幺偏差量纲为p.u.稳态时为0。Eqq轴暂态电动势单位p.u.。连续时间状态方程写成下面的形式dδ/dt ωb × ΔωdΔω/dt (Pm - Pe - D × Δω) / MdEq/dt (Efd - Eq - (xd - xd) × Id) / Td0其中ωb 2π×50 ≈ 314.159 rad/sM为惯性时间常数我取M8s相当于惯性常数H4sD为阻尼系数Pm为机械功率Efd为励磁电压Td0为励磁绕组暂态时间常数。电磁功率Pe和d轴电流Id的表达式为Pe Eq × V∞ × sin(δ) / XdΣId (Eq - V∞ × cos(δ)) / XdΣ这里V∞是无穷大母线电压XdΣ是发电机暂态电抗加上外部线路、变压器电抗的总和单位均为标幺值p.u.。这一步是整个动态估计的状态方程来源滤波器的预测步靠的就是这组微分方程。我把仿真参数列在下面如果你手头有别的电力系统教科书例子可以按对应位置替换但要注意保证稳态平衡点自洽。参数符号物理含义数值f0系统额定频率50 Hzωb同步转速314.159 rad/sM惯性时间常数8 sD阻尼系数2 p.u.xdd轴同步电抗1.8 p.u.xdd轴暂态电抗0.3 p.u.Td0励磁绕组时间常数6 sV∞无穷大母线电压1.0 p.u.XdΣ暂态总电抗0.5 p.u.Pm稳态机械功率0.8 p.u.Efd稳态励磁电压1.1 p.u.这组参数下稳态功角大约在δ≈0.42radEq≈0.96p.u.左右发电机的初始平衡点是自洽的。你可以自己算一遍稳态时PePm0.8由Pe公式反推Eqsinδ0.4再由励磁方程反推EfdEq(xd-xd)×Id两边对照一下就能确认参数没矛盾。3.2 量测方程与噪声模型PMU到底告诉我们什么动态状态估计里最理想的量测来源是PMU它能直接或间接得到我们需要的观测信号。我在这套算例中假设量测向量为y [δ; ω; Pe]对应的物理含义是功角可以由PMU相角与参考母线相角做差得到、转速标幺值频率可得、电磁功率机组有功功率。量测方程写出来h(x) [δ; 1 Δω; Eq × V∞ × sin(δ) / XdΣ]注意转速量测这里用的是1Δω即实际转速的标幺形式数值在1附近小幅波动功角量和暂态电动势则不需要额外换算。噪声建模方面我采用高斯白噪声假设。量测噪声协方差R在Matlab里设置为对角阵我实际用的值是R diag([1e-4; 1e-6; 1e-4])对角线三个值对应的物理标准差大约是功角误差0.01弧度转速误差0.001p.u.有功误差0.01p.u.。这在PMU的典型性能范围内但我故意把它调得比某些理想指标大一点目的是让滤波器的降噪效果在仿真曲线上看得更明显。如果你希望更贴近工程实际可以把第一项降到1e-6量级但滤波曲线会更贴近量测噪声本身视觉对比效果差一些。过程噪声协方差Q代表的是模型本身的不可靠程度我取Q diag([1e-5; 1e-6; 1e-5])这个量级表示我对发电机模型比较信任但功角和暂态电动势有一定的模型扰动。后面调参会专门说Q和R的比例关系直接决定滤波是更信模型还是更信量测这是整个动态估计调试的命门。3.3 离散化与采样周期滤波器的时钟基准上面给出的都是连续时间方程而卡尔曼滤波家族处理的是离散时间序列所以仿真时要做离散化处理。最简单的做法是前向欧拉法即x_{k1} x_k Ts × f(x_k)采样周期Ts取0.01s对应PMU常用上报率。欧拉法实现简单在小步长下精度足够演示但如果后续要接高精度的暂态稳定程序对比我更推荐在生成仿真数据时用四阶龙格-库塔(RK4)做积分这一步对滤波器本身的收敛性影响不大主要影响的是真值轨迹的精度。在滤波器内部EKF和UKF都需要在每个采样时刻完成一遍预测和更新。这里的采样周期与你实际使用PMU数据的刷新率直接挂钩——工程上如果PMU是30帧每秒Ts就取0.033s。从我的经验看Ts过大时EKF的线性化误差会快速放大同样一组Q/R参数在0.01s下UKF和EKF差距不大换成0.05s后EKF就开始明显滞后于真值而UKF还能维持可接受的跟随精度。这就是强非线性场景下HQ? UKF的优势在实际采样率受限时更有意义的直接原因。4. EKF/UKF核心代码实现初始化、预测、更新一手调通4.1 仿真数据生成先造一套带噪声的PMU数据写滤波器之前得先有真值和量测数据。下面的Matlab代码片段用于生成一条带扰动的轨迹前1秒系统稳定在初始平衡点1秒时刻机械功率Pm从0.8阶跃到1.0模拟一次扰动让功角和暂态电动势真正摆动起来然后再在真值上叠加热噪声作为量测。% 系统参数 sys.wb 2*pi*50; sys.M 8; sys.D 2; sys.xd 1.8; sys.xdp 0.3; sys.Td0p 6; sys.Vinf 1.0; sys.XdpSum 0.5; sys.EqVinfX sys.Vinf / sys.XdpSum; sys.invXdp 1 / sys.XdpSum; % 稳态初值可先由潮流/时域仿真得到 x0_true [0.42; 0; 0.96]; % 状态方程函数 function dx f_gen(x, u, sys) delta x(1); dw x(2); Eq x(3); Pm u(1); Efd u(2); Pe sys.EqVinfX * Eq * sin(delta); Id sys.invXdp * (Eq - sys.Vinf * cos(delta)); dx zeros(3,1); dx(1) sys.wb * dw; dx(2) (Pm - Pe - sys.D * dw) / sys.M; dx(3) (Efd - Eq - (sys.xd - sys.xdp) * Id) / sys.Td0p; end % 量测方程函数 function y h_gen(x, sys) delta x(1); dw x(2); Eq x(3); Pe sys.EqVinfX * Eq * sin(delta); y [delta; 1 dw; Pe]; end % 离散仿真主循环 Ts 0.01; N 500; % 5秒数据 t (0:N-1) * Ts; x_true zeros(3, N); z zeros(3, N); x_true(:,1) x0_true; for k 1:N-1 if k 100 k 300 % 1s开始机械功率阶跃 Pm_now 1.0; else Pm_now 0.8; end x_true(:,k1) x_true(:,k) Ts * f_gen(x_true(:,k), [Pm_now; 1.1], sys); end % 叠加量测噪声 R diag([1e-4; 1e-6; 1e-4]); for k 1:N z(:,k) h_gen(x_true(:,k), sys) sqrt(R) * randn(3,1); end这里的函数写法在Matlab脚本里可以直接放在文件末尾作为局部函数或者保存成独立函数文件。我实际更推荐保存成独立文件这样后续扩展多机系统时每个模型文件都可以独立维护。4.2 EKF核心更新代码雅可比矩阵是调试重点EKF的实现重心就是两个雅可比矩阵状态转移雅可比F和量测雅可比H。对于三阶单机模型F可以由f_gen逐项求偏导得到H也可以直接推导。我给出Matlab里的实现% 状态方程雅可比 F df/dx function F F_gen(x, u, sys) delta x(1); dw x(2); Eq x(3); Pm u(1); Efd u(2); Pe sys.EqVinfX * Eq * sin(delta); Id sys.invXdp * (Eq - sys.Vinf * cos(delta)); F zeros(3,3); F(1,2) sys.wb; F(2,1) -sys.EqVinfX * Eq * cos(delta) / sys.M; F(2,2) -sys.D / sys.M; F(2,3) -sys.EqVinfX * sin(delta) / sys.M; F(3,1) (sys.xd - sys.xdp) * sys.invXdp * sys.Vinf * sin(delta) / sys.Td0p; F(3,3) (-1 - (sys.xd - sys.xdp) * sys.invXdp) / sys.Td0p; end % 量测方程雅可比 H dh/dx function H H_gen(x, sys) delta x(1); Eq x(3); H zeros(3,3); H(1,1) 1; H(2,2) 1; H(3,1) sys.EqVinfX * Eq * cos(delta); H(3,3) sys.EqVinfX * sin(delta); endEKF主循环如下% 滤波器初始化 x_ekf [0.45; 0.005; 0.98]; % 故意给偏差看收敛能力 P_ekf diag([0.05^2; 0.01^2; 0.05^2]); Q diag([1e-5; 1e-6; 1e-5]); for k 1:N-1 if k 100 k 300 Pm_now 1.0; else Pm_now 0.8; end u_now [Pm_now; 1.1]; % 预测 x_pred x_ekf Ts * f_gen(x_ekf, u_now, sys); F F_gen(x_ekf, u_now, sys); P_pred F * P_ekf * F Q; % 更新 H H_gen(x_pred, sys); K P_pred * H / (H * P_pred * H R); x_ekf x_pred K * (z(:,k1) - h_gen(x_pred, sys)); P_ekf (eye(3) - K * H) * P_pred; end这里的更新公式用的是Matlab的右除 / 代替了显式的矩阵求逆数值上更稳。从我在实际调试中踩过的坑看EKF出问题最常见的原因不是公式写错而是F矩阵里的偏导符号写错尤其是F(3,1)和F(3,3)这类交叉项符号一错P_pred会越推越歪最终K算出来完全背离真实量测。4.3 UKF核心更新代码sigma点传播与权重重构UKF的核心是生成sigma点并传播。对于状态维数n3一共生成7个点权重按照无迹变换公式计算。不需要求雅可比但需要小心sqrtm(P)的数值问题。代码如下% UKF参数 alpha 0.1; beta 2; kappa 0; n 3; lambda alpha^2 * (n kappa) - n; w_m [lambda/(nlambda); 1/(2*(nlambda))*ones(2*n,1)]; w_c [lambda/(nlambda) (1-alpha^2beta); 1/(2*(nlambda))*ones(2*n,1)]; % 初始化 x_ukf [0.45; 0.005; 0.98]; P_ukf diag([0.05^2; 0.01^2; 0.05^2]); for k 1:N-1 if k 100 k 300 Pm_now 1.0; else Pm_now 0.8; end u_now [Pm_now; 1.1]; % 生成sigma点 S sqrtm((nlambda) * P_ukf); X_sig zeros(n, 2*n1); X_sig(:,1) x_ukf; for i 1:n X_sig(:,i1) x_ukf S(:,i); X_sig(:,in1) x_ukf - S(:,i); end % 状态传播预测 X_pred zeros(n, 2*n1); for i 1:2*n1 X_pred(:,i) X_sig(:,i) Ts * f_gen(X_sig(:,i), u_now, sys); end x_pred zeros(n,1); for i 1:2*n1 x_pred x_pred w_m(i) * X_pred(:,i); end P_pred Q; for i 1:2*n1 diff X_pred(:,i) - x_pred; P_pred P_pred w_c(i) * (diff * diff); end P_pred (P_pred P_pred) / 2; % 量测传播 Y_pred zeros(3, 2*n1); for i 1:2*n1 Y_pred(:,i) h_gen(X_pred(:,i), sys); end y_pred zeros(3,1); for i 1:2*n1 y_pred y_pred w_m(i) * Y_pred(:,i); end Pyy R; Pxy zeros(n, 3); for i 1:2*n1 dy Y_pred(:,i) - y_pred; dx X_pred(:,i) - x_pred; Pyy Pyy w_c(i) * (dy * dy); Pxy Pxy w_c(i) * (dx * dy); end Pyy (Pyy Pyy) / 2; K Pxy / Pyy; x_ukf x_pred K * (z(:,k1) - y_pred); P_ukf P_pred - K * Pyy * K; P_ukf (P_ukf P_ukf) / 2; end这段代码里我刻意没有把alpha取到极端小的0.001而是取了0.1。原因很实际lambda接近-n时权重项的负值会很大sqrtm里(nlambda)P的缩放矩阵会变得非常小数值上对舍入误差更敏感。对于三阶单机模型alpha0.1已经能覆盖非线性传播如果你的系统状态维数更高可以考虑改成alpha1配合beta2效果也足够好。4.4 仿真结果对比与调参经验把EKF和UKF跑同一组数据画在同一张图上最直观的结果是在稳态阶段两条曲线都贴在真值附近肉眼看差别不大但在1秒扰动发生的暂态段EKF的跟踪会出现一个明显的超调或滞后而UKF能更快地回到真值附近。这不是偶然本质原因是EKF的线性化截断误差在功角大范围摆动时被放大而UKF的sigma点传播保留了更多的非线性信息。我建议你们在跑代码时不只画状态量还要画每个时刻的滤波残差z减去量测预测。残差序列如果大致随机、没有明显趋势说明Q和R的比例是合适的如果残差在扰动后持续偏离零很多拍说明Q取得太小滤波器过度信任模型追不上真实变化。我自己的调参顺序是先固定一个R调Q决定跟随速度和光滑度的平衡再反过来调R迭代两次就能收敛到合理范围。后面的调试章节会展开讲。5. 调试中最容易踩的5个坑协方差、初值与Q/R平衡5.1 sqrtm报错协方差矩阵非正定UKF代码里调用sqrtm(P)时Matlab偶尔会报Matrix must be square positive definite之类的错误。原因通常是数值误差让P_ukf失去对称正定性。解决的办法我在代码里已经埋了伏笔在每次更新后强制对称化P_ukf (P_ukf P_ukf) / 2同时也可以加上一个极小的正则项P_ukf P_ukf 1e-10 * eye(n);这里加正则项的幅度不能太大否则会过度抬高量测的信任度让滤波结果出现额外噪声。这个处理本质上是给协方差矩阵垫一层数值地板避免它被舍入误差击穿。5.2 初值给得不对滤波器直接发散动态状态估计不是从零开始的盲估计它非常依赖合理的初始状态和初始协方差。我在这套算例里故意把初值从[0.42;0;0.96]改成[0.45;0.005;0.98]目的就是展示滤波器的收敛能力。对于EKF如果初值偏差超过一定幅度比如功角偏差超过0.1rad在暂态扰动时F矩阵的线性化误差会把P_pred带偏滤波器很可能收敛不到真值附近而是围绕一个错误状态小幅振荡。工程上的正确做法是先用静态状态估计或潮流计算得到各台发电机的功角和电动势作为初始值并把P0设置得略大一些比如对角元取对应状态稳态方差的5到10倍让滤波器在前几个采样周期有足够自由度去修正初始偏差。不要真的把P0设成单位阵甚至0P0太小会让滤波器死死咬住错误的初值。5.3 滤波结果抖成心电图多半是R设小了在调参阶段我见过太多人把R直接设成几千万分之一结果滤波残差很小但估计曲线疯狂抖动根本没法用。原理很简单R越小卡尔曼增益K越大滤波器越相信每一拍量测甚至把量测噪声原样放进估计结果R越大滤波曲线越平滑但对真实动态的跟随会变慢产生滞后。直观的判断办法把仿真量测z的波形和滤波结果的波形叠画如果滤波曲线几乎和带噪声的量测重合抖动幅度跟噪声幅度一个量级就是R太小如果扰动发生后滤波曲线要好几拍才能拐弯追上真值就是Q太小或者R太大。对三阶单机模型我给的R、Q对角值就是相对合理的起点你可以在此基础上做±一个数量级的微调。5.4 状态量数值差异过大导致协方差病态这是很多仿真里隐蔽的一个坑功角δ的单位是弧度数值在0到1左右暂态电动势Eq标幺值在1附近而转速偏差Δω的数值可能是0.001到0.01这一档。三个状态放在同一个协方差矩阵里量级相差上百倍P矩阵天然接近病态再经过多次平方根运算或雅可比乘法数值问题就更容易暴露。我通常在建模时就用标幺形式表示转速偏差而不是用rad/s这样可以缩小状态量级差。如果你确实需要以rad/s为单位那就必须在Q里按量级差异设置不同的对角值比如转速项Q降到1e-7甚至1e-8同时P0里对应的对角元也要跟着调节。否则滤波器的权重分配就完全乱套了。5.5 从单机算例扩展到多机系统时的注意点单机模型跑通之后很多人想直接推多机。这里要提醒一下几个工程细节多机系统的状态方程里每台发电机的功角是相对的需要选定参考机否则状态方程不可观量测方程里PMU量测到的有功从发电机出口功率换到内部电动势端要考虑变压器和线路电抗的分压状态维数n变大后UKF的sigma点数量变成2n1虽然增长可控但每次传播都是完整的非线性系统模型求值计算量会明显上升。这时可以先用EKF跑在线闭环再用UKF离线对比精度工程上更省算力。我自己的一个实用建议是扩展多机时不要用完整的发电机六阶模型起步先用三阶模型把滤波框架跑通确认每台机的状态可观测、量测配置充足再逐步加详细励磁系统和速度调节器模型。这样即使出现问题也能很容易地定位到具体某台机甚至某个状态方程。至少我在实际项目里这种先三阶、再六阶的路径比一开始就上全模型要顺滑得多。最后再分享一个小技巧。调试EKF和UKF时我习惯在同一个循环里同时保存EKF、UKF和真值并实时计算两个滤波器的均方根误差RMSE。这样不仅能看到谁更准还能看到误差是出在稳态还是暂态。如果误差主要在暂态段问题大概率在非线性和Q设置如果稳态也有系统性偏差那就要回头检查量测方程是不是有模型结构错误。这个检查习惯救过我很多次你们可以直接用在现有代码的循环里。 电力系统动态状态估计：EKF与UKF仿真实现与调试实战
1. 从静态断面到动态递推电力系统动态状态估计到底解决了什么问题在电力系统里状态估计这个词绝大多数场景下指的是调度中心SCADA/EMS里的加权最小二乘静态估计把一个电网断面里的母线电压幅值、相角、支路潮流和注入功率作为量测量解一个线性或弱非线性方程组得到一个稳态快照。这套办法在正常运行方式下非常可靠调度员看的就是那张大屏上的潮流分布。但问题来了当系统真的发生扰动——线路跳闸、机组跳机、大负荷突增——功角和转速开始围绕稳态点来回摆动的那几秒到几十秒静态状态估计的假设直接被击穿。断面量测还没算完实际状态已经跑了十万八千里。真要等到稳态恢复再算对广域保护、紧急控制来说黄花菜都凉了。所以我第一次接触电力系统动态状态估计(DSE)时心里很清楚这不是静态估计的升级版而是完全换了一套数学工具把发电机本身的微分方程拿进来用滤波器的思路逐拍递推估计功角、转速、暂态电动势这类真正随动态过程变化的状态变量。支撑DSE落地的硬件基础是这些年在电网里铺开的PMU同步相量测量单元。PMU以10到60帧每秒的速度上报带时标的相量和频率时间分辨率远高于传统RTU遥测。有了高频量测才谈得上动态估计的闭环校正。进一步说DSE的核心算法就是我们的标题主角扩展卡尔曼滤波EKF以及它之后发展出的无迹卡尔曼滤波UKF。EKF的思路最直观把非线性方程在工作点附近做一阶泰勒展开用雅可比矩阵硬套线性卡尔曼滤波公式UKF则绕开求导用一组sigma点去传播概率分布对强非线性模型更宽容。这篇文章我会从一个可复现的Matlab仿真入手把动态状态估计这整套链路讲清楚怎么搭发电机三阶模型、怎么构造量测方程、EKF和UKF的代码到底怎么写、哪些参数会让滤波发散、怎么从单机算例扩展到多机系统。适合的读者包括正在做电力系统课程设计或毕业设计的学生也包括刚接手广域测量系统数据分析、需要在现有PMU数据上试算动态估计的研究人员和工程师。我会把我在调试中踩过的坑一并写出来你们可以照着抄作业也可以拿来当避坑清单。需要先说明的是这是基于我过往项目经验的通用实践如果你要结合真实电网的详细模型需要在这个框架上修改参数矩阵和量测配置但核心的滤波更新流程基本不变。2. EKF与UKF的原理拆解线性化与sigma点的取舍2.1 EKF把非线性问题局部线性化的工程妥协卡尔曼滤波原本是为线性系统设计的在线性高斯假设下它是最小方差意义下的最优估计器。但电力系统动态模型是明显的非线性常微分方程组比如功角运动方程里有sin(δ)励磁方程里还叠着电压和电流的乘积项根本不可能用线性Kalman那套严格推导。EKF做的妥协很直白在每个时刻的工作点附近对状态方程和量测方程求雅可比矩阵用一阶偏导把非线性问题强行压成线性问题再套标准卡尔曼递推。EKF的标准递推分五步我直接写出来。假设离散状态方程是 x_k f(x_{k-1}, u_{k-1}) w量测方程是 z_k h(x_k) v过程噪声协方差Q量测噪声协方差R预测x_pred f(x_prev, u)P_pred F * P_prev * F Q其中F是状态方程雅可比在x_prev处求值。更新卡尔曼增益K P_pred * H * (H * P_pred * H R)^{-1}其中H是量测方程雅可比在x_pred处求值滤波值x_est x_pred K * (z - h(x_pred))协方差更新P_est (I - K*H) * P_pred。EKF的优点非常明显实现简单计算量在众多非线性滤波算法里几乎最小一个状态维度为n的系统主要代价是生成两个n×n雅可比矩阵并做两次矩阵求逆。我在Matlab里跑单机系统的时候EKF单步耗时基本在毫秒级完全能满足PMU的实时性要求。但EKF的短板也很致命当系统非线性较强、或者工作点变化剧烈时一阶泰勒展开的截断误差会迅速累积。最典型的情况是功角在暂态过程中摆动幅度很大sin(δ)的高阶项不再能忽略EKF的协方差估计就开始失真严重时直接发散。还有一个工程上的痛点雅可比矩阵推导容易出错尤其在多机系统里状态方程涉及所有发电机的互作用项手推F矩阵推错一两个符号滤波结果在视觉上看着还行但协方差早就病态了。2.2 UKF不用求导用sigma点传播统计量UKF的设计动机非常聪明既然直接求非线性函数的导数容易出错且只保到一阶那干脆用一组精心挑选的采样点sigma点去感受非线性函数的形状这些点在原分布的均值附近按协方差信息展开经过非线性函数传播后从它们的新位置重新计算均值和协方差。这样一来完全不需要推导雅可比矩阵而且统计精度理论上能捕获到非线性函数的二阶项比EKF的一阶局部分解更贴合真实分布。生成sigma点的标准做法是对于n维状态一共取2n1个点。设状态均值为x_m协方差为P引入比例因子λ α²(n κ) - n其中α通常取0.01到1之间的小正数κ一般取0那么sigma点构成为第0个点x_0 x_m第i个点i1,...,nx_i x_m sqrt((nλ)P)的第i列第in个点x_i x_m - sqrt((nλ)P)的第i列。权重系数方面均值权重为W_m0 λ/(nλ)其余点权重为W_mi 1/(2(nλ))协方差权重的第0项额外叠加一个(1-α²β)修正项β一般取2其余项与均值权重相同。把这组sigma点代入非线性状态函数和量测函数之后用加权平均重构预测均值和协方差再用与标准卡尔曼相同的增益公式做更新。UKF的实际计算量大概是EKF的2到3倍因为每个时刻要传播2n1组sigma点每组都是一次完整的非线性函数求值。但代价换来的是强非线性场景下明显更稳的结果而且代码里不需要手推任何雅可比矩阵。换句话说EKF的难点在建模时的推导UKF的难点在实现时的矩阵运算细节两者对工程人员的要求确实不同。2.3 两种滤波器的性能对照我用下面这张表总结两者差异这是我在实际仿真中最直观的感受也符合多数文献里的共识对比维度EKFUKF是否需要雅可比矩阵需要推导易错不需要零求导非线性处理精度一阶泰勒近似约二阶精度单步计算量较小适合高维状态约2~3倍EKF随维数线性增长强非线性暂态场景可能出现协方差失稳更稳残差更小实现复杂程度公式简洁但调试花时间代码略长但调参逻辑清晰适合场景弱非线性、模型维数高强非线性、模型故障切换频繁从我自己的调试经验看单机或小规模系统上UKF比EKF的优越性不一定体现在快或准上更多体现在你不需要把每一个偏导数都算对。省下的推导时间远比多跑的那几毫秒有价值。这也是我在文末会强烈建议初学者先跑通UKF的原因。3. 仿真模型与Matlab参数设计先把发电机三阶模型搭对3.1 单机无穷大母线系统建模基准做动态状态估计的仿真我不建议一上来就搭完整多机系统工程习惯是先跑通一个单机无穷大母线模型。它结构简单却保留了发电机动态的核心非线性特征转子摆动方程和励磁绕组的暂态过程。所谓无穷大母线就等价于一个幅值、频率恒定不变的大电网发电机挂在这上面功角摆动会直接反映电磁功率的变化正是我们想用滤波器跟踪的动态。我采用的经典三阶单轴发电机模型状态向量选为x [δ; Δω; Eq]三个状态分别对应δ发电机功角单位弧度。Δω转子转速相对同步转速的标幺偏差量纲为p.u.稳态时为0。Eqq轴暂态电动势单位p.u.。连续时间状态方程写成下面的形式dδ/dt ωb × ΔωdΔω/dt (Pm - Pe - D × Δω) / MdEq/dt (Efd - Eq - (xd - xd) × Id) / Td0其中ωb 2π×50 ≈ 314.159 rad/sM为惯性时间常数我取M8s相当于惯性常数H4sD为阻尼系数Pm为机械功率Efd为励磁电压Td0为励磁绕组暂态时间常数。电磁功率Pe和d轴电流Id的表达式为Pe Eq × V∞ × sin(δ) / XdΣId (Eq - V∞ × cos(δ)) / XdΣ这里V∞是无穷大母线电压XdΣ是发电机暂态电抗加上外部线路、变压器电抗的总和单位均为标幺值p.u.。这一步是整个动态估计的状态方程来源滤波器的预测步靠的就是这组微分方程。我把仿真参数列在下面如果你手头有别的电力系统教科书例子可以按对应位置替换但要注意保证稳态平衡点自洽。参数符号物理含义数值f0系统额定频率50 Hzωb同步转速314.159 rad/sM惯性时间常数8 sD阻尼系数2 p.u.xdd轴同步电抗1.8 p.u.xdd轴暂态电抗0.3 p.u.Td0励磁绕组时间常数6 sV∞无穷大母线电压1.0 p.u.XdΣ暂态总电抗0.5 p.u.Pm稳态机械功率0.8 p.u.Efd稳态励磁电压1.1 p.u.这组参数下稳态功角大约在δ≈0.42radEq≈0.96p.u.左右发电机的初始平衡点是自洽的。你可以自己算一遍稳态时PePm0.8由Pe公式反推Eqsinδ0.4再由励磁方程反推EfdEq(xd-xd)×Id两边对照一下就能确认参数没矛盾。3.2 量测方程与噪声模型PMU到底告诉我们什么动态状态估计里最理想的量测来源是PMU它能直接或间接得到我们需要的观测信号。我在这套算例中假设量测向量为y [δ; ω; Pe]对应的物理含义是功角可以由PMU相角与参考母线相角做差得到、转速标幺值频率可得、电磁功率机组有功功率。量测方程写出来h(x) [δ; 1 Δω; Eq × V∞ × sin(δ) / XdΣ]注意转速量测这里用的是1Δω即实际转速的标幺形式数值在1附近小幅波动功角量和暂态电动势则不需要额外换算。噪声建模方面我采用高斯白噪声假设。量测噪声协方差R在Matlab里设置为对角阵我实际用的值是R diag([1e-4; 1e-6; 1e-4])对角线三个值对应的物理标准差大约是功角误差0.01弧度转速误差0.001p.u.有功误差0.01p.u.。这在PMU的典型性能范围内但我故意把它调得比某些理想指标大一点目的是让滤波器的降噪效果在仿真曲线上看得更明显。如果你希望更贴近工程实际可以把第一项降到1e-6量级但滤波曲线会更贴近量测噪声本身视觉对比效果差一些。过程噪声协方差Q代表的是模型本身的不可靠程度我取Q diag([1e-5; 1e-6; 1e-5])这个量级表示我对发电机模型比较信任但功角和暂态电动势有一定的模型扰动。后面调参会专门说Q和R的比例关系直接决定滤波是更信模型还是更信量测这是整个动态估计调试的命门。3.3 离散化与采样周期滤波器的时钟基准上面给出的都是连续时间方程而卡尔曼滤波家族处理的是离散时间序列所以仿真时要做离散化处理。最简单的做法是前向欧拉法即x_{k1} x_k Ts × f(x_k)采样周期Ts取0.01s对应PMU常用上报率。欧拉法实现简单在小步长下精度足够演示但如果后续要接高精度的暂态稳定程序对比我更推荐在生成仿真数据时用四阶龙格-库塔(RK4)做积分这一步对滤波器本身的收敛性影响不大主要影响的是真值轨迹的精度。在滤波器内部EKF和UKF都需要在每个采样时刻完成一遍预测和更新。这里的采样周期与你实际使用PMU数据的刷新率直接挂钩——工程上如果PMU是30帧每秒Ts就取0.033s。从我的经验看Ts过大时EKF的线性化误差会快速放大同样一组Q/R参数在0.01s下UKF和EKF差距不大换成0.05s后EKF就开始明显滞后于真值而UKF还能维持可接受的跟随精度。这就是强非线性场景下HQ? UKF的优势在实际采样率受限时更有意义的直接原因。4. EKF/UKF核心代码实现初始化、预测、更新一手调通4.1 仿真数据生成先造一套带噪声的PMU数据写滤波器之前得先有真值和量测数据。下面的Matlab代码片段用于生成一条带扰动的轨迹前1秒系统稳定在初始平衡点1秒时刻机械功率Pm从0.8阶跃到1.0模拟一次扰动让功角和暂态电动势真正摆动起来然后再在真值上叠加热噪声作为量测。% 系统参数 sys.wb 2*pi*50; sys.M 8; sys.D 2; sys.xd 1.8; sys.xdp 0.3; sys.Td0p 6; sys.Vinf 1.0; sys.XdpSum 0.5; sys.EqVinfX sys.Vinf / sys.XdpSum; sys.invXdp 1 / sys.XdpSum; % 稳态初值可先由潮流/时域仿真得到 x0_true [0.42; 0; 0.96]; % 状态方程函数 function dx f_gen(x, u, sys) delta x(1); dw x(2); Eq x(3); Pm u(1); Efd u(2); Pe sys.EqVinfX * Eq * sin(delta); Id sys.invXdp * (Eq - sys.Vinf * cos(delta)); dx zeros(3,1); dx(1) sys.wb * dw; dx(2) (Pm - Pe - sys.D * dw) / sys.M; dx(3) (Efd - Eq - (sys.xd - sys.xdp) * Id) / sys.Td0p; end % 量测方程函数 function y h_gen(x, sys) delta x(1); dw x(2); Eq x(3); Pe sys.EqVinfX * Eq * sin(delta); y [delta; 1 dw; Pe]; end % 离散仿真主循环 Ts 0.01; N 500; % 5秒数据 t (0:N-1) * Ts; x_true zeros(3, N); z zeros(3, N); x_true(:,1) x0_true; for k 1:N-1 if k 100 k 300 % 1s开始机械功率阶跃 Pm_now 1.0; else Pm_now 0.8; end x_true(:,k1) x_true(:,k) Ts * f_gen(x_true(:,k), [Pm_now; 1.1], sys); end % 叠加量测噪声 R diag([1e-4; 1e-6; 1e-4]); for k 1:N z(:,k) h_gen(x_true(:,k), sys) sqrt(R) * randn(3,1); end这里的函数写法在Matlab脚本里可以直接放在文件末尾作为局部函数或者保存成独立函数文件。我实际更推荐保存成独立文件这样后续扩展多机系统时每个模型文件都可以独立维护。4.2 EKF核心更新代码雅可比矩阵是调试重点EKF的实现重心就是两个雅可比矩阵状态转移雅可比F和量测雅可比H。对于三阶单机模型F可以由f_gen逐项求偏导得到H也可以直接推导。我给出Matlab里的实现% 状态方程雅可比 F df/dx function F F_gen(x, u, sys) delta x(1); dw x(2); Eq x(3); Pm u(1); Efd u(2); Pe sys.EqVinfX * Eq * sin(delta); Id sys.invXdp * (Eq - sys.Vinf * cos(delta)); F zeros(3,3); F(1,2) sys.wb; F(2,1) -sys.EqVinfX * Eq * cos(delta) / sys.M; F(2,2) -sys.D / sys.M; F(2,3) -sys.EqVinfX * sin(delta) / sys.M; F(3,1) (sys.xd - sys.xdp) * sys.invXdp * sys.Vinf * sin(delta) / sys.Td0p; F(3,3) (-1 - (sys.xd - sys.xdp) * sys.invXdp) / sys.Td0p; end % 量测方程雅可比 H dh/dx function H H_gen(x, sys) delta x(1); Eq x(3); H zeros(3,3); H(1,1) 1; H(2,2) 1; H(3,1) sys.EqVinfX * Eq * cos(delta); H(3,3) sys.EqVinfX * sin(delta); endEKF主循环如下% 滤波器初始化 x_ekf [0.45; 0.005; 0.98]; % 故意给偏差看收敛能力 P_ekf diag([0.05^2; 0.01^2; 0.05^2]); Q diag([1e-5; 1e-6; 1e-5]); for k 1:N-1 if k 100 k 300 Pm_now 1.0; else Pm_now 0.8; end u_now [Pm_now; 1.1]; % 预测 x_pred x_ekf Ts * f_gen(x_ekf, u_now, sys); F F_gen(x_ekf, u_now, sys); P_pred F * P_ekf * F Q; % 更新 H H_gen(x_pred, sys); K P_pred * H / (H * P_pred * H R); x_ekf x_pred K * (z(:,k1) - h_gen(x_pred, sys)); P_ekf (eye(3) - K * H) * P_pred; end这里的更新公式用的是Matlab的右除 / 代替了显式的矩阵求逆数值上更稳。从我在实际调试中踩过的坑看EKF出问题最常见的原因不是公式写错而是F矩阵里的偏导符号写错尤其是F(3,1)和F(3,3)这类交叉项符号一错P_pred会越推越歪最终K算出来完全背离真实量测。4.3 UKF核心更新代码sigma点传播与权重重构UKF的核心是生成sigma点并传播。对于状态维数n3一共生成7个点权重按照无迹变换公式计算。不需要求雅可比但需要小心sqrtm(P)的数值问题。代码如下% UKF参数 alpha 0.1; beta 2; kappa 0; n 3; lambda alpha^2 * (n kappa) - n; w_m [lambda/(nlambda); 1/(2*(nlambda))*ones(2*n,1)]; w_c [lambda/(nlambda) (1-alpha^2beta); 1/(2*(nlambda))*ones(2*n,1)]; % 初始化 x_ukf [0.45; 0.005; 0.98]; P_ukf diag([0.05^2; 0.01^2; 0.05^2]); for k 1:N-1 if k 100 k 300 Pm_now 1.0; else Pm_now 0.8; end u_now [Pm_now; 1.1]; % 生成sigma点 S sqrtm((nlambda) * P_ukf); X_sig zeros(n, 2*n1); X_sig(:,1) x_ukf; for i 1:n X_sig(:,i1) x_ukf S(:,i); X_sig(:,in1) x_ukf - S(:,i); end % 状态传播预测 X_pred zeros(n, 2*n1); for i 1:2*n1 X_pred(:,i) X_sig(:,i) Ts * f_gen(X_sig(:,i), u_now, sys); end x_pred zeros(n,1); for i 1:2*n1 x_pred x_pred w_m(i) * X_pred(:,i); end P_pred Q; for i 1:2*n1 diff X_pred(:,i) - x_pred; P_pred P_pred w_c(i) * (diff * diff); end P_pred (P_pred P_pred) / 2; % 量测传播 Y_pred zeros(3, 2*n1); for i 1:2*n1 Y_pred(:,i) h_gen(X_pred(:,i), sys); end y_pred zeros(3,1); for i 1:2*n1 y_pred y_pred w_m(i) * Y_pred(:,i); end Pyy R; Pxy zeros(n, 3); for i 1:2*n1 dy Y_pred(:,i) - y_pred; dx X_pred(:,i) - x_pred; Pyy Pyy w_c(i) * (dy * dy); Pxy Pxy w_c(i) * (dx * dy); end Pyy (Pyy Pyy) / 2; K Pxy / Pyy; x_ukf x_pred K * (z(:,k1) - y_pred); P_ukf P_pred - K * Pyy * K; P_ukf (P_ukf P_ukf) / 2; end这段代码里我刻意没有把alpha取到极端小的0.001而是取了0.1。原因很实际lambda接近-n时权重项的负值会很大sqrtm里(nlambda)P的缩放矩阵会变得非常小数值上对舍入误差更敏感。对于三阶单机模型alpha0.1已经能覆盖非线性传播如果你的系统状态维数更高可以考虑改成alpha1配合beta2效果也足够好。4.4 仿真结果对比与调参经验把EKF和UKF跑同一组数据画在同一张图上最直观的结果是在稳态阶段两条曲线都贴在真值附近肉眼看差别不大但在1秒扰动发生的暂态段EKF的跟踪会出现一个明显的超调或滞后而UKF能更快地回到真值附近。这不是偶然本质原因是EKF的线性化截断误差在功角大范围摆动时被放大而UKF的sigma点传播保留了更多的非线性信息。我建议你们在跑代码时不只画状态量还要画每个时刻的滤波残差z减去量测预测。残差序列如果大致随机、没有明显趋势说明Q和R的比例是合适的如果残差在扰动后持续偏离零很多拍说明Q取得太小滤波器过度信任模型追不上真实变化。我自己的调参顺序是先固定一个R调Q决定跟随速度和光滑度的平衡再反过来调R迭代两次就能收敛到合理范围。后面的调试章节会展开讲。5. 调试中最容易踩的5个坑协方差、初值与Q/R平衡5.1 sqrtm报错协方差矩阵非正定UKF代码里调用sqrtm(P)时Matlab偶尔会报Matrix must be square positive definite之类的错误。原因通常是数值误差让P_ukf失去对称正定性。解决的办法我在代码里已经埋了伏笔在每次更新后强制对称化P_ukf (P_ukf P_ukf) / 2同时也可以加上一个极小的正则项P_ukf P_ukf 1e-10 * eye(n);这里加正则项的幅度不能太大否则会过度抬高量测的信任度让滤波结果出现额外噪声。这个处理本质上是给协方差矩阵垫一层数值地板避免它被舍入误差击穿。5.2 初值给得不对滤波器直接发散动态状态估计不是从零开始的盲估计它非常依赖合理的初始状态和初始协方差。我在这套算例里故意把初值从[0.42;0;0.96]改成[0.45;0.005;0.98]目的就是展示滤波器的收敛能力。对于EKF如果初值偏差超过一定幅度比如功角偏差超过0.1rad在暂态扰动时F矩阵的线性化误差会把P_pred带偏滤波器很可能收敛不到真值附近而是围绕一个错误状态小幅振荡。工程上的正确做法是先用静态状态估计或潮流计算得到各台发电机的功角和电动势作为初始值并把P0设置得略大一些比如对角元取对应状态稳态方差的5到10倍让滤波器在前几个采样周期有足够自由度去修正初始偏差。不要真的把P0设成单位阵甚至0P0太小会让滤波器死死咬住错误的初值。5.3 滤波结果抖成心电图多半是R设小了在调参阶段我见过太多人把R直接设成几千万分之一结果滤波残差很小但估计曲线疯狂抖动根本没法用。原理很简单R越小卡尔曼增益K越大滤波器越相信每一拍量测甚至把量测噪声原样放进估计结果R越大滤波曲线越平滑但对真实动态的跟随会变慢产生滞后。直观的判断办法把仿真量测z的波形和滤波结果的波形叠画如果滤波曲线几乎和带噪声的量测重合抖动幅度跟噪声幅度一个量级就是R太小如果扰动发生后滤波曲线要好几拍才能拐弯追上真值就是Q太小或者R太大。对三阶单机模型我给的R、Q对角值就是相对合理的起点你可以在此基础上做±一个数量级的微调。5.4 状态量数值差异过大导致协方差病态这是很多仿真里隐蔽的一个坑功角δ的单位是弧度数值在0到1左右暂态电动势Eq标幺值在1附近而转速偏差Δω的数值可能是0.001到0.01这一档。三个状态放在同一个协方差矩阵里量级相差上百倍P矩阵天然接近病态再经过多次平方根运算或雅可比乘法数值问题就更容易暴露。我通常在建模时就用标幺形式表示转速偏差而不是用rad/s这样可以缩小状态量级差。如果你确实需要以rad/s为单位那就必须在Q里按量级差异设置不同的对角值比如转速项Q降到1e-7甚至1e-8同时P0里对应的对角元也要跟着调节。否则滤波器的权重分配就完全乱套了。5.5 从单机算例扩展到多机系统时的注意点单机模型跑通之后很多人想直接推多机。这里要提醒一下几个工程细节多机系统的状态方程里每台发电机的功角是相对的需要选定参考机否则状态方程不可观量测方程里PMU量测到的有功从发电机出口功率换到内部电动势端要考虑变压器和线路电抗的分压状态维数n变大后UKF的sigma点数量变成2n1虽然增长可控但每次传播都是完整的非线性系统模型求值计算量会明显上升。这时可以先用EKF跑在线闭环再用UKF离线对比精度工程上更省算力。我自己的一个实用建议是扩展多机时不要用完整的发电机六阶模型起步先用三阶模型把滤波框架跑通确认每台机的状态可观测、量测配置充足再逐步加详细励磁系统和速度调节器模型。这样即使出现问题也能很容易地定位到具体某台机甚至某个状态方程。至少我在实际项目里这种先三阶、再六阶的路径比一开始就上全模型要顺滑得多。最后再分享一个小技巧。调试EKF和UKF时我习惯在同一个循环里同时保存EKF、UKF和真值并实时计算两个滤波器的均方根误差RMSE。这样不仅能看到谁更准还能看到误差是出在稳态还是暂态。如果误差主要在暂态段问题大概率在非线性和Q设置如果稳态也有系统性偏差那就要回头检查量测方程是不是有模型结构错误。这个检查习惯救过我很多次你们可以直接用在现有代码的循环里。