发电机最优励磁控制系统设计讲解

发电机最优励磁控制系统设计
1 线性最优控制理论
1.1 线性最优控制原理
线性最优控制理论所研究的核心问题是选取最优的控制规律，使控制系统在特定指标下的性能为最优。

控制系统框如图1所示。

图1 最优控制框图
线性定常系统状态空间方程的一般形式为：
()X t AX BU =+ （1）
式中，A 为状态系数矩阵；B 为控制系数矩阵；X 为n 维状态向量；U 为r 维控制向量。

如图1所示，如果要改善系统性能，可以引入状态反馈构成闭环系统。

反馈系统的状态向量为：
U KX =- （2）
式中，K 为状态反馈增益矩阵。

将式（1）和式（2）合并，可以得到：
()()X AX B KX A BK X =+-=- （3）
可见，最优控制的本质就是如何选取反馈矩阵K ，以使它在给定控制规律下达到特定条件下的最优。

1.2 二次型性能指标
假定()y t 为系统的实际响应，()t ξ为系统预期的响应。

最优控制性能指标是使两者的偏差最小，即
[]2
min 0
()()J t y t dt J ξ∞
=-=⎰ （4）
式中，J 是随函数()y t 而变化的一个泛函数，是在0∞时间内求取偏差平
方的积分，称为二次型指标。

如果以()X t 为实际状态向量，以ˆ()X t 为预期状态向量，则要求状态向量偏差最小的二次型性能指标为：
ˆˆ()()()()T
J X t X t X t X t dt ∞
⎡⎤⎡⎤
=--⎣⎦⎣⎦⎰
（5）
在二次型性能指标中，也需要引入对控制量的约束。

如果没有这个约束，所设计的控制器中，控制量的变化范围可能会很大，难以实现。

引入控制量约束的二次型指标为：
{
}
ˆˆ()()()()()()T
T J X t X t Q X t X t U t RU t dt ∞⎡⎤⎡⎤=--+⎣⎦⎣⎦
⎰
（6） 式中Q 和R 表示状态向量和控制向量的权矩阵。

为了便于分析，通常把系统平衡点置于状态空间的原点，即ˆ()0X t =，则式（6）可以变换为：
0()()()()T T
J X t QX t U t RU t dt ∞
⎡⎤=+⎣
⎦⎰ （7） 以上式作为性能指标设计最优控制系统，可以证明这个最优控制规律是存在且唯一的，表达式为：
1()()()T U t R B PX t KX t -=-=- （8）
式中，K 为最优反馈增益矩阵。

P 为n n ⨯维对称常矩阵，矩阵P 是黎卡梯方程：10T T A P PA PBR B P Q -+-+=的解。

2 单机无穷大系统数学模型
与无穷大系统并联运行的同步发电机的基本数学模型包括4个部分：与系统并联的同步发电机基本方程组，发电机电势、定子电压和电流之间的联系方程，功率方程以及转子绕组动态方程。

2.1与系统并联的同步发电机基本方程组
与系统并联的同步发电机基本方程组决定了运行中的发电机的稳态与过渡过程的全部动力学特性。

转子运动方程的形式如下：
22
2m e D H d P P P f dt δ
π=-- （9）
式中，H 用秒，δ用弧度，t 用秒，m P 、e P 、D P 均为标幺值。

在实用计算中，一般采用02D D d P f dt
δ
π=
，这里的D 是机组的阻尼系数，用标幺值，单位用弧度。

2.2发电机电势、定子电压和电流之间的联系方程
将观察发电机的坐标系从abc 坐标系转换到dq0坐标系，可以极大地简化发电机的数学模型。

这是因为abc 坐标系固定在三相对称的定子上，由于转子的旋转，对于每一坐标轴方向上的磁路磁导率都是时间t 的周期函数，不便于计算和分析。

而dq0坐标系是固定在转子上，坐标轴与磁路相对静止，这样磁导率就是一常数，因而惦记的各个参数也会成为与时间t 无关的常数。

由发电机电势与电流矢量图可以直接得到如下关系式组：
()
cos sin q sq d d
q
q d d d q
tq d d sq q d d sq q d d sd q q sq s sd s s E V I x E E I x x E V I x V E I x V E I x V I x V V V V V δδ∑∑∑=+⎧⎪''=+-⎪''⎪=+⎪
=-⎪⎪''
=-⎨⎪=⎪
=⎪⎪=⎪
⎪=⎩
（10） 由上述的发电机定子电势、电压和电流的关系式组，可以解出定子电流为：
cos q s d d E V I x δ
∑
-=
（11）
sin s q q V I x δ
∑
=
（12） cos q
s d d
E V I x δ∑''-=
' （13） 由矢量图又可以得到：
cos tq s d s V V I x δ=+ （14） sin tq s q s V V I x δ=+ （15）
t V = （16）
代入式（11）和（12），可以解得：
t V =
（17）
2.3 功率方程
与无穷大系统并联的发电机有功功率的表达式为：
2sin sin 22q s d q s e q sq d sd d d q E V x x V P I V I V x x x δδ∑∑
∑∑∑
-=+=
+ （18）
以q
E '和δ为自变量的电磁功率e P 表达式如下： 2sin sin 22q s d
q s e d d
q E V x x V P x x x δδ∑∑∑∑∑''-=+'' （19）
2.4 转子绕组动态方程
在列写转子绕组的实用动态方程时，忽略阻尼绕组，只列写励磁绕组动态方程。

另外假设定子绕组中磁链式可以突变的，励磁绕组磁链守恒，这是因为定子绕组中暂态电流的非周期分量比励磁绕组非周期分量衰减快得多。

带负载的发电机在转子d 轴上的总磁链为：
0fd d fl ad f ad f fl d ad i L i L I L ψψψψ=+-=+- （20）
式中，0d ψ是有用磁通；fl ψ是转子漏磁通；ad ψ是电枢反应磁通；ad L 是励磁绕组与定子d 轴绕组的互感；fl L 是励磁绕组的漏感。

假定过渡过程中转速标幺值为*1ω=，式（20）可以写为：
fd f ad f fl d ad i x i x I x ψ=+- （21）
可以求得：
fd d ad
f f
I x i x ψ+=
（22）
此处，f ad fl x x x =+为励磁绕组的自感抗。

假设定子d 轴总磁链与端电压q 轴分量标幺值相等。

类似于转子绕组，可以得到：
2
()ad ad
tq fd
d d f f
x x V I x x x ψ=-- （23） 这里d ad l x x x =+为发电机d 轴同步电抗，其中l x 是定子漏电抗；2
ad
d d f
x x x x '=-为
d 轴暂态电抗；ad
q
fd f
x E x ψ'=称为d 轴暂态电抗后电势。

由励磁绕组等值电路可得回路方程式：
f f f f f
di V i R x dt
=+ （24）
以
ad
f
x R 乘以上式，可得在空载情况下励磁绕组的动态方程： 0
q f q d dE E E T dt
=+ （25）
上式中，0f d f
x T R =
为定子绕组开路时励磁绕组的时间常数。

若以
ad
f
x R 乘以式（24），可得： 0
q f q d dE E E T dt
'=+ （26）
以上就是得到了主发电机励磁绕组的动态方程。

2.5 基本方程组的偏差化与线性化
由式（9）可以得到偏差方程：
22m e D m e H D
P P P P P f f ωδππ∆=∆-∆-∆=∆-∆-∆ （27） 根据式（18）和式（19），将电磁功率方程进行偏差化和线性化，可得：
e E E q P S R E δ∆=∆+∆ （28）
e E E q P S R E δ'''∆=∆+∆ （29）
式中，2
cos cos2q s d q E s
d d q E V x x S V
x x x δδ∑∑∑
∑∑-=
-；2cos cos2q s d q E s d d q E V x x S V x x x δδ∑
∑
'∑∑∑
''-=-''；
sin s E d V R x δ∑=
；sin s E d
V
R x δ∑='。

将式（28）和（29）代入式（27），得到线性化的转子运动方程：
0022m E E q H D P S R E f f ωδωππ∆=∆-∆-∆-∆ （30） 00
22m E E q H D
P S R E f f ωδωππ'''∆=∆-∆-∆-∆ （31）
由t t t q q
V V
V E E δδ∂∂∆=
∆+∆∂∂，可得到运行角偏差δ∆和机端电压偏差t V ∆表示的电磁功率偏差：
e V V t P S R V δ∆=∆+∆ （32）
其中，t V E V
V S S R δ∂=-∂，E
V t q
R R V E =∂∂。

于是转子的运动偏差方程有：
00
22m V V t H D
P S R V f f ωδωππ∆=∆-∆-∆-∆ （33） 励磁绕组的动态偏差方程为：
0f q d q
E E T E '∆=∆+∆ （34） 3 最优励磁控制系统设计
由式（33）和式（34）可以得到系统状态方程为：
00
0001
0000011001100
E E E E R q q d E d d f f e e D S R H H H S S
V E E T R T T E E T T ω
ωδδωω'
'
⎡⎤⎢⎥⎡⎤∆⎡⎤⎢⎥∆⎡⎤---
⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥∆∆⎢⎥⎢⎥⎢⎥-⎢⎥=+∆-
⎢⎥⎢⎥'⎢⎥∆'∆⎢⎥'⎢⎥⎢⎥⎢⎥⎢⎥
∆∆⎢⎥⎢⎥⎢⎥⎣⎦⎣
⎦⎢⎥⎣⎦⎢⎥-⎢⎥⎣
⎦
（35） 由于q
E '不便于测量，所以采用机端电压t V ∆来代替q E '，即 E E V q
t E E S S R
E V R R δ'''
-'∆=∆+∆ （36）
又由于在实际系统中采用的是e P ∆、ω∆、t V ∆和f E ∆的一组状态量，通过状态变量置换，再结合式（35），可以得到e P ∆、ω∆、t V ∆和f E ∆所表示的状态方程：
00
00000
001100
E V
V E E E d V
d V d
e e R
E V E V
E E t t d V V V d V d V f f e e S S R S R S T S T S T P P D H H V S S S S S R V V T R S R T S T R E E T T ωωω''''-⎡⎤-
⎢⎥
''⎢⎥⎡⎤∆⎡⎤∆⎡⎤⎢⎥⎢⎥--⎢⎥⎢⎥⎢⎥⎢⎥∆∆⎢⎥⎢⎥⎢
⎥⎢⎥=+∆⎢⎥--⎢⎥∆∆⎢⎥⎢⎥-
⎢⎥⎢⎥⎢
⎥⎢⎥'∆∆⎢⎥⎢⎥⎣⎦⎣⎦⎢⎥⎢⎥⎣⎦⎢⎥
-⎢⎥⎣
⎦
（37） 如果发电机采用可控硅励磁系统，则可以忽略励磁时间常数0e T ≈，那么控制量U 就是励磁绕组电压f E ∆，系统状态方程可以表示为X AX BU =+形式的三阶形式：
000000E V V E E E d V d V d e e f t t E E V
E V
E d V d V V
V d V S S R S S R T S
T S T P P D E H H V V R S S S S S
T R T R S R T S ωωω''''-⎡⎤
-
⎡⎤
⎢⎥
''⎢⎥
⎢⎥⎡⎤∆∆⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥∆=--
∆+∆⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥∆∆⎣⎦⎣⎦⎢⎥--⎢⎥-⎢⎥⎣⎦⎢⎥'⎣⎦
（38） 最优控制量U kX =-
4 最优控制设计算例
发电机电气参数：隐极机纵轴电抗 2.5d X =，纵轴暂态电抗0.5d X '=，纵轴
开环时间常数010d T s =，8H s =，综合阻尼系数 5.0D =；变压器参数：0.01T X =；
线路参数：1
1.52
L X =⨯；系统参数：1s V =。

所研究的发输电系统的示意图见图
2。

s
图2 系统示意图
设计最优运行点为065δ=︒，00.55P =时的最优控制器。

4.1 计算最优控制参数
通过式（18）和式（19），根据最优运行点065δ=︒和00.55P =，计算得到最优运行点的相关参数如下： 由式（18）可得：00.55(2.50.010.75)
1.9784sin 1sin 65e d q s P x E V δ∑⨯++=
==⨯︒
由式（19）可得：20(sin 2) 1.02392sin d q d s
q e d q s x x x V E P x x V δδ
∑∑∑∑∑'-''=-='
由式（17）
可得：0 1.0487t V =
=
数据代入到相应公式，可以求得
-0.16950.6564-0.1422-39.2699-0.6250 0-0.1064-0.0545-0.0893A ⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦； 0.071900.0452B ⎡⎤
⎢⎥=⎢⎥
⎢⎥⎣⎦
选取(1,100,5000)Q diag =和1R =作为权矩阵，进而可以得到状态反馈矩阵：
[]46.5370-6.156068.76(,,)1,7K lqr A B Q R ==
可知最优控制量58.7480 5.176730.0717f e t u E KX P V ω=∆=-=-∆+∆-∆。

上述计算过程使用MATLAB 完成，见报告附录中的程序1。

4.2 验证最优励磁控制系统性能
由于励磁控制系统参数是采用线性化状态方程而得到的，而实际系统的状态方程是非线性的，所以使用MATLAB 的ode45算法，对单机系统中发电机的运动微分方程进行计算，通过观察发电机的转子角和输出电磁功率的变化情况，
得到最优励磁控制系统的控制性能。

单机系统的发电机运动方程如下：
0δωω=-
20
00()(sin sin 2)2q s d
q s m d d
q E V x x D V P H H H x x x ωωωωωδδ∑∑∑∑∑''-=---+''
00111cos d
d q
f q s d d d d
x x E E E V T T T x δ∑'-''=-+''
此运动方程以转子角δ、角速度ω和q 轴暂态电势q
E '作为三个状态量。

微分方程中，02100f ωππ==；0.55m P =；f E 作为输入量，引入最优状态反馈矩阵，即用三个状态量进行表示，为
0000000()(58.7480 5.176746.5370() 6.30.0711560()68.7617()
)7)(f q f q q e t e t q t E E E E KX E P P P V V V E ωωω=+∆=+⨯-+⨯--⨯--=+-∆+∆
-∆=+-
其中，e P 和t V 作为与q E '和δ有关的变量，分别用式（19）和式（17）表示。

4.2.1 静态稳定性
系统在运行至5s 时刻，加入一个时间长达1秒的电磁功率的扰动，模拟负荷的突然变化，即5s 到6s 之间e P 从原来的0.55增加到0.6，6s 之后e P 回到最初设计的稳态运行点。

进行小扰动的仿真。

仿真程序见报告附录中的程序2~7和13。

仿真的功角和电磁功率波形图如下：
图3 静态稳定功角仿真图
从波形图可以看出，在发生小扰动后，采用最优控制理论所设计的励磁系
统，能够快速地平息系统的振荡，转子角和电磁功率均能够迅速地回到最初设计的最优运行点，提高了电力系统的静态稳定性。

可见，发电机加入最优励磁控制器可以有效地改善系统稳定性。

4.2.2 暂态稳定性
系统运行至5s的时刻，在其中一条线路首端即图2中的d处，发生三相短路故障，故障0.1s后保护动作切除故障线路。

程序见报告附录中的程序2、3和8~12。

该过程中的功角变化波形图如下：
图4 暂态稳定功角仿真图
由图4可以看出，系统在5s时刻发生三相短路后功角和电磁功率急剧上升，5.1s切除故障线路后系统迅速地恢复稳定，整个过程不超过5秒钟。

因此，最优励磁控制器也可以很好地维持电力系统的暂态稳定性。

5 结论
本报告中先介绍了最优控制理论的相关原理，然后通过对与电网并联的发电机基本方程组的研究从而建立了发电机与励磁系统的线性化状态方程，采用MATLAB计算最优状态反馈参数并代入到系统的非线性微分方程组，然后对系统的非线性微分方程组进行了时域仿真。

仿真结果表明，最优励磁控制器可以很好地提高系统的静态稳定性和暂态稳定性。

可见，最优励磁控制器有效地增加了发电机励磁系统的阻尼，抑制了系统受到扰动时的振荡，加快系统恢复稳定的速度。

附录
程序1：
%输入发电机与系统的相关参数
Xd=2.5;
Xdz=0.5;
XT=0.01;
Xq=Xd;
XL=0.5*1.5;
Td0=10;
Vs=1;
D=5;
H=8;
Pe=0.55;
ang=65*pi/180;
f0=50;
%中间计算过程
Xe=XT+XL;
Xds=Xd+XT+XL;
Xdzs=Xdz+XT+XL;
Xqs=Xq+XT+XL;
Eq=Pe*Xds/(Vs*sin(ang))
Eqz=(Pe-0.5*Vs^2*(Xdzs-Xqs)/(Xdzs*Xqs)*sin(2*ang))*Xdzs/(Vs*sin(ang))
Td=Td0*Xdzs/Xds;
SE=Eq*Vs*cos(ang)/Xds+Vs*Vs*(Xds-Xqs)/(Xds*Xqs)*cos(2*ang);
SEz=Eqz*Vs*cos(ang)/Xdzs+Vs*Vs*(Xdzs-Xqs)/(Xdzs*Xqs)*cos(2*ang);
RE=Vs*sin(ang)/Xds;
REz=Vs*sin(ang)/Xdzs;
Vt=sqrt((Eq^2*Xe^2+(Vs*cos(ang)*Xd)^2+2*Xe*Xd*Eq*Vs*cos(ang))/Xds^2+(Vs*sin(ang)* Xq/Xqs)^2)
dVtEq=((Eq*Xe^2+Xe*Xd*Vs*cos(ang))/(Xds^2))*((Eq^2*Xe^2+Vs^2*(cos(ang))^2*Xd^2+2 *Xe*Xd*Eq*Vs*cos(ang))/(Xds^2)+(Vs*sin(ang)*Xq/Xqs)^2)^(-0.5);
dVtang=0.5*(Vs^2*Xq^2*sin(2*ang)/Xqs^2-(Vs^2*Xd^2*sin(2*ang)+2*Xe*Xd*Eq*Vs*sin(an g))/Xds^2)*((Eq^2*Xe^2+(Vs*cos(ang)*Xd)^2+2*Xd*Xe*Eq*Vs*cos(ang))/Xds^2+(Vs*sin(an g)*Xq/Xqs)^2)^(-0.5);
RV=RE/dVtEq;
SV=SE-RV*dVtang;
w0=2*pi*f0;
%计算得到A,B,K矩阵
Q=[1 0 0;0 100 0;0 0 5000];
R=1;
A=[(SE-SV)/(Td*SV) SEz -RV*SE/(Td*SV);-w0/H -D/H 0;(SE-SV)/(Td*RV*SV) (SEz-SV)/RV -SE/(Td*SV)]
B=[REz/Td0;0;REz/(Td0*RV)]
K=lqr(A,B,Q,R)
程序2：
%故障前或扰动前的系统状态
function dy=gz1(t,y);
Xd=2.5;Xdz=0.5
XT=0.01;
XL=0.75;
Xe=XT+XL;
Xds=Xd+Xe;
Xdzs=Xdz+Xe;
Vs=1;
Eq0=1.9784;
Pe=y(3)*Vs*sin(y(1))/Xdzs+0.5*Vs^2*(Xdzs-Xds)/(Xdzs*Xds)*sin(2*y(1));
Eq=Xds*(y(3)/Xdzs+Vs*(Xdzs-Xds)/(Xdzs*Xds)*cos(y(1)));
Vt=1/Xds*sqrt(Eq^2*Xe^2+Vs^2*Xd^2+2*Xe*Xd*Eq*Vs*cos(y(1)));
Td0=10;
Td=Td0*Xdzs/Xds;
dy(1)=y(2)-100*pi;
dy(2)=2*pi*50*0.55/8-5/8*(y(2)-2*pi*50)-2*pi*50/8*Pe;
dy(3)=(-46.537*(Pe-0.55)+6.156*(y(2)-2*pi*50)-68.7617*(Vt-1.0487)+Eq0)/Td0-y(3)/Td+1/Td 0*(Xd-Xdz)/Xdzs*Vs*cos(y(1));
dy=[dy(1);dy(2);dy(3)];
程序3：
%计算故障前或扰动前的电磁功率
function Pe1=dcgl1(y)
Pe1=zeros(length(y(:,1)),1);
Xd=2.5;Xdz=0.5;
XT=0.01;
XL=0.75;
Xe=XT+XL;
Xds=Xd+Xe;
Xdzs=Xdz+Xe;
Vs=1;
for i=1:length(y(:,1));
Pe1(i)=y(i,3)*Vs*sin(y(i,1))/Xdzs+0.5*Vs^2*(Xdzs-Xds)/(Xdzs*Xds)*sin(2*y(i,1));
End
程序4：
%某线路首端三相短路故障中的系统状态
function dy=gz2(t,y);
Xd=2.5;Xdz=0.5
XT=0.01;
Xe=XT;
Xds=Xd+Xe;
Xdzs=Xdz+Xe;
Vs=1;
Pe=0;
Eq=Xds*(y(3)/Xdzs+Vs*(Xdzs-Xds)/(Xdzs*Xds)*cos(y(1)));
Vt=1/Xds*sqrt(Eq^2*Xe^2+Vs^2*Xd^2+2*Xe*Xd*Eq*Vs*cos(y(1)));
Eq0=1.9784;
Td0=10;
Td=Td0*Xdzs/Xds;
dy(1)=y(2)-100*pi;
dy(2)=2*pi*50*0.55/8-5/8*(y(2)-2*pi*50)-2*pi*50/8*Pe;
dy(3)=(-46.537*(Pe-0.55)+6.156*(y(2)-2*pi*50)-68.7617*(Vt-1.0487)+Eq0)/Td0-y(3)/Td+1/Td 0*(Xd-Xdz)/Xdzs*Vs*cos(y(1));
dy=[dy(1);dy(2);dy(3)];
程序5：
%计算故障中的电磁功率
function Pe4=dcgl4(y)
Pe4=zeros(length(y(:,1)),1);
for i=1:length(y(:,1));
Pe4(i)=0;
End
程序6：
%切除故障线路后的系统状态
function dy=gz3(t,y);
Xd=2.5;Xdz=0.5
XT=0.01;
XL=1.5;
Xe=XT+XL;
Xds=Xd+Xe;
Xdzs=Xdz+Xe;
Vs=1;
Pe=y(3)*Vs*sin(y(1))/Xdzs+0.5*Vs^2*(Xdzs-Xds)/(Xdzs*Xds)*sin(2*y(1));
Eq=Pe*Xds/(Vs*sin(y(1)));
Vt=1/Xds*sqrt(Eq^2*Xe^2+Vs^2*Xd^2+2*Xe*Xd*Eq*Vs*cos(y(1)));
Eq0=1.9784;
Td0=10;
Td=Td0*Xdzs/Xds;
dy(1)=y(2)-100*pi;
dy(2)=2*pi*50*0.55/8-5/8*(y(2)-2*pi*50)-2*pi*50/8*Pe;
dy(3)=(-46.537*(Pe-0.55)+6.156*(y(2)-2*pi*50)-68.7617*(Vt-1.0487)+Eq0)/Td0-y(3)/Td+1/Td 0*(Xd-Xdz)/Xdzs*Vs*cos(y(1));
dy=[dy(1);dy(2);dy(3)];
程序7：
%计算故障后的电磁功率
function Pe5=dcgl5(y)
Pe5=zeros(length(y(:,1)),1);
Xd=2.5;Xdz=0.5;
XT=0.01;
XL=1.5;
Xe=XT+XL;
Xds=Xd+Xe;
Xdzs=Xdz+Xe;
Vs=1;
for i=1:length(y(:,1));
Pe5(i)=y(i,3)*Vs*sin(y(i,1))/Xdzs+0.5*Vs^2*(Xdzs-Xds)/(Xdzs*Xds)*sin(2*y(i,1));
end
程序8：
%发生扰动的系统状态
function dy=rd(t,y);
Xd=2.5;Xdz=0.5;
XT=0.01;
XL=0.75;
Xe=XT+XL;
Xds=Xd+Xe;
Xdzs=Xdz+Xe;
Vs=1;
Eq0=1.9784;
Pe=y(3)*Vs*sin(y(1))/Xdzs+0.5*Vs^2*(Xdzs-Xds)/(Xdzs*Xds)*sin(2*y(1));
Eq=Xds*(y(3)/Xdzs+Vs*(Xdzs-Xds)/(Xdzs*Xds)*cos(y(1)));
Vt=1/Xds*sqrt(Eq^2*Xe^2+Vs^2*Xd^2+2*Xe*Xd*Eq*Vs*cos(y(1)));
Td0=10;
Td=Td0*Xdzs/Xds;
dy(1)=y(2)-100*pi;
dy(2)=2*pi*50*0.6/8-5/8*(y(2)-2*pi*50)-2*pi*50/8*Pe;
dy(3)=(-46.537*(Pe-0.6)+6.156*(y(2)-2*pi*50)-68.7617*(Vt-1.0487)+Eq0)/Td0-y(3)/Td+1/Td0 *(Xd-Xdz)/Xdzs*Vs*cos(y(1));
dy=[dy(1);dy(2);dy(3)];
程序9：
%计算扰动中的电磁功率
function Pe2=dcgl2(y)
Pe2=zeros(length(y(:,1)),1);
Xd=2.5;Xdz=0.5;
XT=0.01;
XL=0.75;
Xe=XT+XL;
Xds=Xd+Xe;
Xdzs=Xdz+Xe;
Vs=1;
for i=1:length(y(:,1));
Pe2(i)=y(i,3)*Vs*sin(y(i,1))/Xdzs+0.5*Vs^2*(Xdzs-Xds)/(Xdzs*Xds)*sin(2*y(i,1));
end
程序10：
%扰动后的系统状态
function dy=rdh(t,y);
Xd=2.5;Xdz=0.5;
XT=0.01;
XL=0.75;
Xe=XT+XL;
Xds=Xd+Xe;
Xdzs=Xdz+Xe;
Vs=1;
Eq0=1.9784;
Pe=y(3)*Vs*sin(y(1))/Xdzs+0.5*Vs^2*(Xdzs-Xds)/(Xdzs*Xds)*sin(2*y(1));
Eq=Xds*(y(3)/Xdzs+Vs*(Xdzs-Xds)/(Xdzs*Xds)*cos(y(1)));
Vt=1/Xds*sqrt(Eq^2*Xe^2+Vs^2*Xd^2+2*Xe*Xd*Eq*Vs*cos(y(1)));
Td0=10;
Td=Td0*Xdzs/Xds;
dy(1)=y(2)-100*pi;
dy(2)=2*pi*50*0.55/8-5/8*(y(2)-2*pi*50)-2*pi*50/8*Pe;
dy(3)=(-46.537*(Pe-0.55)+6.156*(y(2)-2*pi*50)-68.7617*(Vt-1.0487)+Eq0)/Td0-y(3)/Td+1/Td 0*(Xd-Xdz)/Xdzs*Vs*cos(y(1));
dy=[dy(1);dy(2);dy(3)];
程序11：
%计算扰动后的电磁功率
function Pe3=dcgl3(y)
Pe3=zeros(length(y(:,1)),1);
Xd=2.5;Xdz=0.5;
XT=0.01;
XL=0.75;
Xe=XT+XL;
Xds=Xd+Xe;
Xdzs=Xdz+Xe;
Vs=1;
for i=1:length(y(:,1));
Pe3(i)=y(i,3)*Vs*sin(y(i,1))/Xdzs+0.5*Vs^2*(Xdzs-Xds)/(Xdzs*Xds)*sin(2*y(i,1)); end
程序12：
%设定仿真时间为40秒,5秒时电磁功率受到扰动,6秒时刻扰动结束
t0=0;
t1=5;
t2=6;
t3=40;
%扰动前，以选定的最优运行点为初值
[T1,Y1]=ode45(@gz1,[t0 t1],[65*pi/180 100*pi 1.0487]);
Pe1=dcgl1(Y1);
ac=Y1(length(Y1),1);
wc=Y1(length(Y1),2);
ec=Y1(length(Y1),3);
%扰动中，以扰动前状态量的最后一组解为扰动中的初值
[T2,Y2]=ode45(@rd,[t1 t2],[ac wc ec]);
Pe2=dcgl2(Y2);
ac2=Y2(length(Y2),1);
wc2=Y2(length(Y2),2);
ec2=Y2(length(Y2),3);
%扰动后，以扰动中状态量的最后一组解为扰动后的初值
[T3,Y3]=ode45(@rdh,[t2 t3],[ac2 wc2 ec2]);
Pe3=dcgl3(Y3);
subplot(2,1,1);
plot(T1,Y1(:,1)/pi*180,T2,Y2(:,1)/pi*180,T3,Y3(:,1)/pi*180);xlabel('t/s');ylabel('angle'); subplot(2,1,2);
plot(T1,Pe1,T2,Pe2,T3,Pe3);xlabel('t/s');ylabel('Pe');
程序13：
%设定仿真时间为40秒,5秒时发生三相接地故障，故障0.1秒后故障切除
t0=0;
t1=5;
t2=5.1;
t3=40;
%故障前，以选定的最优运行点为初值
[T1,Y1]=ode45(@gz1,[t0 t1],[65*pi/180 100*pi 1.0487]);
Pe1=dcgl1(Y1);
ac=Y1(length(Y1),1);
wc=Y1(length(Y1),2);
ec=Y1(length(Y1),3);
%故障中，以故障前的状态量的最后一组解为故障中的初值
[T2,Y2]=ode45(@gz2,[t1 t2],[ac wc ec]);
Pe4=dcgl4(Y2);
ac2=Y2(length(Y2),1);
wc2=Y2(length(Y2),2);
ec2=Y2(length(Y2),3);
%故障后，以故障中的状态量的最后一组解为故障中的初值
[T3,Y3]=ode45(@gz3,[t2 t3],[ac2 wc2 ec2]);
Pe5=dcgl5(Y3);
subplot(2,1,1);
plot(T1,Y1(:,1)/pi*180,T2,Y2(:,1)/pi*180,T3,Y3(:,1)/pi*180);xlabel('t/s');ylabel('angle'); subplot(2,1,2);
plot(T1,Pe1,T2,Pe4,T3,Pe5);xlabel('t/s');ylabel('Pe');。