实验三 电力系统暂态稳定分析
电力系统暂态稳定计算实际上就是求解发电机转子运动方程的初值问题，从而得出δ-t 和ω-t 的关系曲线。每台发电机的转子运动方程是两个一阶非线性的常微分方程。因此，首先介绍常微分方程的初值问题的数值解法。
一、常微分方程的初值问题 （一）问题及求解公式的构造方法
我们讨论形如式（3-1）的一阶微分方程的初值问题
??
?=≤≤='00
)(),,()(y x y b
x a y x f x y （3-1） 设初值问题（3-1）的解为)(x y ，为了求其数值解而采取离散化方法，在求解区间[b a ,]上取一组节点
b x x x x x a n i i =<<<<<<=+ 110
称i i i x x h -=+1（1,,1,0-=n i ）为步长。在等步长的情况下，步长为
n
a
b h -=
用i y 表示在节点i x 处解的准确值)(i x y 的近似值。
设法构造序列{}i y 所满足的一个方程（称为差分方程）
),,(1h y x h y y i i i i ??+=+ （3-2）
作为求解公式，这是一个递推公式，从（0x ，0y ）出发，采用步进方式，自左相右逐步算出)(x y 在所有节点i x 上的近似值i y （n i ,,2,1 =）。
在公式（3-2）中，为求1+i y 只用到前面一步的值i y ，这种方法称为单步法。在公式（3-2）中的1+i y 由i y 明显表示出，称为显式公式。而形如（3-3）
),,,(11h y y x h y y i i i i i ++?+=ψ （3-3）
的公式称为隐式公式，因为其右端ψ中还包括1+i y 。
如果由公式求1+i y 时，不止用到前一个节点的值，则称为多步法。 由式（3-1）可得
dy =dx y x f ),( （3-4）
两边在[i x ，1+i x ]上积分，得
?
++
=+1
))(,()()(1i i
x x i i dx x y x f x y x y （3-5）
由此可以看出，如果想构造求解公式，就要对右端的积分项作某种数值处理。这种求解公式的构造方法叫做数值积分法。
（二）一般的初值问题的解法
1． 欧拉法和改进欧拉法
对于初值问题（3-1），采用数值积分法，从而得到（3-5）。对于（3-5）右端的积分用矩形公式(取左端点)，则得到
?
+?≈1
))(,())(,(i i
x x i i x y x f h dx x y x f
进而得到（3-1）的求解公式（3-2）
),(1i i i i y x f h y y ?+=+ （i =0，1，2，n-1） （3-6）
此公式称为欧拉（Euler ）格式。
如果对式（3-5）右端的积分用梯形公式
)))(,())(,((2
))(,(111
++?
++?≈
i i x x i i x y x f x y x f h
dx x y x f i i
则可以得到初值问题（3-1）的梯形求解公式如式（3-7）
[]),(),(2
111=+++?+
=i i i i i i y x f y x f h
y y （i =0，1，2，n-1） （3-7） 式（3-7）是个隐式公式。可以采取先用欧拉格式求一个)(1+i x y 的初步近似值，记作1+i y ，称之为预报值，然后用预报值1+i y 替代式（3-7）右端的1+i y ，再计算得到1+i y ，称之为校正值，这样建立起来的预报－校正方法称为改进欧拉格式
[]??
?
??+?+=?+=++++),(),(2),(1111i i i i i i i i i i y x f y x f h y y y x f h y y （3-8）
2． 龙格—库塔方法
在单步法中，应用最广泛的是龙格－库塔（Runge-kutta ）法，简称R －K 法。下面直接给出一种四阶的龙格－库塔法的计算公式（3-9）
????
???????
++?=++?=++?=?=++++==)
,()
21,2()
21
,2()
,()22(61342
3121
43211K y h x f h K K y h x f h K K y h x f h K y x f h K K K K K y y i i i i i i i i i i （3-9） 它也称为标准（古典）龙格－库塔法。 例3-1 研究下列微分方程的初值问题
?????=-+='0
)0(2112
2y y x
y 解：
这是一个特殊的微分方程，其解的解析式可以给出，为
2
1x x
y +=
应用龙格－库塔法，取h =0.25，根据式（3-9）编写一段程序，由零开始自左相右逐步算出)(x y 在所有节点i x 上的近似值i y 。计算结果见表3-1。计算结果表明，四阶龙格－库塔方法的精度是较高的。
表3-1
实际上，MATLAB 为常微分方程提供了很好的解题指令，使得求解常微分方程变得很容易，并且能将问题及解答表现在图形上。因此，我们可以不用根据式（3-9）编写较复杂的程序，而只需应用MATLAB 提供的常微分方程解题器来解决问题。下面给出用MATLAB 编写的解题程序。
首先编写描述常微分方程的ODE 文件，文件名为ˊmyfun ˊ，便于解题器调用它。
function dy = myfun(x,y) dy = zeros(1,1); dy=1/(1+x^2)-2*y^2;
再编写利用解题器指令求解y 的程序。
clear x0=0; for i=1:4 xm=2*i;
y0=0;
[x,y] = ode45('myfun',[x0 xm],[y0]); format long y(length(y)) end
plot(x,y,'-')
运行上述程序，在得到几个点的函数值的同时，也得到函数y 的曲线，如图3-1所示。
图3-1 根据运算结果画出y 的曲线
二、简单电力系统的暂态稳定性 （一）物理过程分析
某简单电力系统如图3-2(a)所示，正常运行时发电机经过变压器和双回线路向无限大系统供电。发电机用电势E
' 作为其等值电势，则电势E ' 与无限大系统间的电抗为 212
T L
T d
x x x x x +++'=I （3-10） 这时发电机发出的电磁功率可表示为
δδsin sin M P x U
E P I I
I ='=
（3-11） 如果突然在一回输电线路始端发生不对称短路，如图3-2(b)所示。故障期间发电机电势E
' 与无限大系统之间的联系电抗为 ?
II ++'++++'=x x x x x x x x x x T L
T d
T L
T d
)2)(()2
()(2121 （3-12）
在故障情况下发电机输出的电磁功率为
δδsin sin M P x U
E P II II
II ='=
（3-13） 在短路故障发生之后，线路继电保护装置将迅速断开故障线路两端的断路器，如图3-2(c)所示。此时发电机电势E
' 与无限大系统间的联系电抗为
21T L T d
x x x x x +++'=III （3-14） 发电机输出的功率为
δδsin sin M P x U
E P III III
III ='=
（3-15） U =c
E
~
G T 1
T 2
L
~
jx jx jx x j 'U
E
' (c )
L
jx 1
T jx 2T jx d
x j 'U
图3-2 简单电力系统及其等值电路
（a ）正常运行方式及其等值电路；（b ）故障情况及其等值电路；（c ）故障切除后及其等值电路
如果正常时发电机向无限大系统输送的有功功率为0P ，则原动机输出的机械功率T P 等于0P 。假定不计故障后几秒种之调速器的作用，即认为机械功率始终保持0P 。因此，可以得到此简单电力系统正常运行、故障期间及故障切除后的功率特性曲线如图3-3所示。
0δk δc
δm δh δ0
P P T =P
δ
图3-3 简单系统正常运行、故障期间及故障切除后的功率特性曲线
对于上述简单电力系统，我们可以根据等面积定则求得极限切除角。但是，实际工作需要知道在多少时间之切除故障线路，也就是要知道与极限切除角对应的极限切除时间。要解决这个问题，必须求解发电机的转子运动方程。 （二）求解发电机的转子运动方程
求解发电机转子运动方程可以得出δ-t 和ω-t 的关系曲线。其中δ-t 曲线一般称为摇摆曲线。在上述简单电力系统中故障期间的转子运动方程为
???
????-=-=II )sin (1
)1(1
δωωωδ
M T J P P T dt d dt
d （3-16） 式中，δ——功率角，其单位为弧度；ω——转子角速度，标幺值；1ω——转子的同步角速度，即1ω=f π2=314.16，其单位为弧度/秒；J T ——发电机的惯性时间常数，其单位为秒；T P 、M P II ――分别为机械和电磁功率，标幺值。
这是两个一阶的非线性常微分方程，它的起始条件是已知的，即
t =0t =0； ω=0ω=1.0；δ=0δ=M
T
P P I -1
sin 故障切除后，由于系统参数改变，以致发电机功率特性发生变化，必须开始求解另一组微分方程：
???
????-=-=III )sin (1
)1(1
δωωωδ
M T J P P T dt d dt
d （3-17） 式中变量含义同前述，其中M P III 也为标幺值。这组方程的起始条件为
t =c t ；δ=c δ；ω=c ω
其中c t 为给定的切除时间；c δ、c ω为与c t 时刻对应的δ和ω，它们可由故障期间的δ-t 和ω-t 的关系曲线求得（δ和ω都是不突变的）。一般来说，在计算故障发生后几秒种的过程中，如果δ始终不超过180o，而且振荡幅值越来越小，则系统是暂态稳定的。
当发电机与无限大系统之间发生振荡或失去同步时，在发电机的转子回路中，特别是阻尼绕组中将有感应电流而形成阻尼转矩（也称为异步转矩）。当作微小振荡时，阻尼功率可表达为：
D P =ω?D =)1(-ωD （3-18）
式中，D 称为阻尼功率系数；ω?为转子角速度的偏移量，标幺值；ω为转子角速度，标幺值。阻尼功率系数D 除了与发电机的参数有关外，还和原始功角、δ?的振荡频率有关。在一般情况下它是正数。在原始功角较小，或者定子回路中有串联电容使定子回路总电阻相对于总电抗较大时，D 可能为负数。如果考虑阻尼功率的影响，则故障后的转子运动方程又可表达为
???
????---=-=III ]sin )1([1
)1(1
δωωωωδ
M T J P D P T dt d dt
d （3-19） 电力系统暂态稳定计算包括两类问题，一类是应用数值计算法得出故障期间的曲线后，根据曲线找到与极限切除角对应的极限切除时间，此时只需要求解微分方程（3-16）；另一类是已知故障切除时间，需要求出摇摆曲线来判断系统的稳定性，此时需要分段分别求解微分方程（3-16）和（3-17）。如果考虑阻尼转矩的影响，则此时需要分段分别求解微分方程（3-16）和（3-19）。 三、例题
例3-2 某简单电力系统如图3-4所示，取基准值B S =220MV A ，B U =209KV 。换算后的参
数已经标在图中，其中一回线的电抗L x =0.486，J T =8.18秒。设电力线路某一回的始端发生两相接地短路。假定E '=常数。（1）计算保持暂态稳定而要求的极限切除角。（2）计算极限切除时间，并且作出在0.15秒切除故障时的δ-t 曲线。
~
G
T 1
T 2
L
k
d
x '1T x L x 2T x 0P 0Q U =0.122=1.0=0.486
=0.138
=0.432=0.295L
L x x 40==1.0
=0.2
2x
图3-4 某简单电力系统的接线图
解：计算系统正常运行方式，决定E '和0δ。 由3-3(a)的正序网络可得，此时系统的总电抗为
I x =0.295+0.138+0.243+0.122=0.798
发电机的暂态电势为：
E '=2
2)0
.1798.00.1()0.1798.02.00.1(?+?+
=1.41 0δ=1
-tg 798
.02.00.1798
.0?+=34.53o
（2）故障后的功率特性
又由3-3(b)的负序、零序网络可得故障点的负序、零序等值电抗为
∑2x =
)
122.0243.0()138.0432.0()
122.0243.0()138.0432.0(++++?+=0.222
∑0x =
)
122.0972.0(138.0)
122.0972.0(138.0+++=0.123
所以在正序网络故障点上的附加电抗为：
079.0123
.0222.0123
.0222.0=+?=
?x
于是故障时等值电路如图3-3(c)所示，则
80.2079
.0365
.0433.0365.0433.0=++
+=II x
因此，故障期间发电机的最大功率为：
504.08
.20
.141.1=?='=
II II x U E P M （3）故障切除后的功率特性
故障切除后的等值电路如图3-3(d)所示
041.1122.0486.0138.0295.0=+++=III x
此时最大功率为
35.1041
.10
.141.1=?='=
III III x U E P M 01
02.13235
.10
.1sin 180=-=-h δ
2
.00=Q ='E 0=δ
j 0.295j 0.138j 0.486j 0.122
(a )(b )
(c )
)
(d ='E 41.1='E U =1.0
图3-5 例題7-12的等值电路
（a ）正常运行等值电路；（b ）负序和零序等值电路；（c ）故障时等值电路；（d ）故障切除后等值电路
（4）计算极限切除角
M
M M h M h T cm P P P P P II III II III --+-=
0cos cos )(cos δδδδδ
=
504
.035.153.34cos 504.02.132cos 35.1)53.342.132(180
0.10
0--+-?
π
=0.458
cm δ074.62=
（5）找出极限切除时间cm t 根据（3-16），首先计算初值
6027.0180
53
.340=?=
πδ，0.10=ω 令y(1)=δ，y(2)=ω。编写描述故障期间转子运动方程的ODE 文件，文件名为ˊmyequ ˊ。
function dy = myequ(t,y) dy = zeros(2,1); f=50;w1=2*pi*f; dy(1) = (y(2)-1)*w1;
dy(2) = (1/8.18)*(1.0-0.504*sin(y(1)));
再编写利用解题器指令求解y 的程序。
clear
t0=0;tm=0.25;
d0=(34.53/180)*pi;w0=1;
[T,Y] = ode45('myequ',[t0 tm],[d0 w0]); plot(T,(Y(:,1)/pi)*180,'-',0.194,62.76,'*')
text(0.194,60,'\delta_{cmax}=62.76\circ','FontSize',10) text(0.194,56,'t_{cmax}=0.194s','FontSize',10)
图3-6 例題7-12的δ-t 曲线
图3-6给出短路发生后0秒到0.25秒期间的δ-t 计算曲线，根据最大切除角cm
δ
（0
t为0.194秒。由图3-6可见，如果故障切除时间大于0.194 62
=）找到极限切除时间
.
74
cm
秒，则发电机的功角将不断地增大，最终失去暂态稳定。在极限切除时间之前切除故障，发电机的摇摆曲线的状况将在下面作计算、分析。
（6）不考虑阻尼转矩的影响，当故障切除时间为0.15秒时通过计算得出δ-t曲线
首先编写描述故障期间转子运动方程的ODE文件，文件名为”myfun01”。
function dy = myfun01(t,y)
f=50; w1=2*pi*f;
TJ=8.18; Pt=1.0; P2m=0.504;
dy = zeros(2,1);
dy(1) = (y(2)-1)*w1;
dy(2) = (1/TJ)*(Pt-P2m*sin(y(1)));
再编写描述故障切除后转子运动方程的ODE文件，文件名为”myfun02”。
function dy = myfun02(t,y)
f=50; w1=2*pi*f;
TJ=8.18; Pt=1.0; P3m=1.35;
dy = zeros(2,1);
dy(1) = (y(2)-1)*w1;
dy(2) = (1/TJ)*(Pt-P3m*sin(y(1)));
编写利用解题器指令求解y的小程序。
clear
t0=0; tc=0.15; tm=2.0;
d0=(34.53/180)*pi; w0=1.0;
[T1,Y1] = ode45('myfun01',[t0 tc],[d0 w0]);
dc=Y1(length(Y1),1);
wc=Y1(length(Y1),2);
[T2,Y2] = ode45('myfun02',[tc tm],[dc wc]);
plot(T1,(Y1(:,1)/pi)*180,'-',T2,(Y2(:,1)/pi)*180,'-',tc,(dc/pi)*180,' *')
text(0.28,50,'\it{t}_{c}=0.15s','FontSize',8)
text(0.28,43,'\it{\delta}_{c}=51.71\circ','FontSize',8)
xlabel('\it{t}')
ylabel('\it{\delta}')
计算结果表明，功角δ沿着故障切除后的功角特性曲线根据等面积定则作等幅振荡，如图3-7所示。实际上，由于阻尼转矩的影响，振荡的幅度是逐渐衰减的，功角δ最终运行在
δ=47.8o。因此，发电机能够保持暂态稳定。
k
δ曲线
图3-7 不考虑阻尼转矩影响，当0.15秒切除故障时发电机的t-
（7）考虑阻尼转矩的影响，当故障切除时间为0.15秒时通过计算得出δ-t曲线描述故障期间转子运动方程的ODE文件与（6）相同，文件名也为”myfun01”。
重新编写描述故障切除后转子运动方程的ODE文件，文件名为”myfun03”，阻尼功率系数D取为15。
function dy = myfun03(t,y)
f=50; w1=2*pi*f;
TJ=8.18; Pt=1.0; P3m=1.35;
D=15;
dy = zeros(2,1);
dy(1) = (y(2)-1)*w1;
dy(2) = (1/TJ)*(Pt-D*(y(2)-1)-P3m*sin(y(1)));
再编写利用解题器指令求解y的小程序。
clear
t0=0; tc=0.15; tm=4;
d0=34.53*3.14/180; w0=1;
[T1,Y1] = ode45('myfun01',[t0 tc],[d0 w0]);
dc=Y1(length(Y1),1);
wc=Y1(length(Y1),2);
[T2,Y2] = ode45('myfun03',[tc tm],[dc wc]);
plot(T1,(Y1(:,1)/pi)*180,'-',T2,(Y2(:,1)/pi)*180,'-',tc,(dc/pi)*180,' *')
text(0.3,50,'\it{t}_c=0.15s','FontSize',8)
text(0.28,43,'\it{\delta}_{c}=51.71\circ','FontSize',8)
xlabel('\it{t}')
ylabel('\it{\delta}')
图3-8 不考虑阻尼转矩影响，当0.15秒切除故障时发电机的t -δ曲线
计算结果表明，功角δ沿着故障切除后的功角特性曲线作减幅振荡，如图3-8所示。功角δ最终运行在k δ=47.8o。因此，发电机能够保持暂态稳定。