MATLAB实验电力系统暂态稳定分析.docx
《MATLAB实验电力系统暂态稳定分析.docx》由会员分享，可在线阅读，更多相关《MATLAB实验电力系统暂态稳定分析.docx（16页珍藏版）》请在冰豆网上搜索。
MATLAB实验电力系统暂态稳定分析
实验三电力系统暂态稳定分析
电力系统暂态稳定计算实际上就是求解发电机转子运动方程的初值问题，从而得出3-t
和3-t的关系曲线。
每台发电机的转子运动方程是两个一阶非线性的常微分方程。
因此，首先介绍常微分方程的初值问题的数值解法。
-、常微分方程的初值问题
（一）问题及求解公式的构造方法
我们讨论形如式（3-1）的一阶微分方程的初值问题
（3-1）
y（x）f（x,y）,axby（xo）yo
设初值问题（3-1）的解为y（x），为了求其数值解而采取离散化方法，在求解区间[a,b]上
取一组节点
axox1xixi1xnb
称hixi1xi（i0,1,,n1）为步长。
在等步长的情况下，步长为
ba
h
n
用yi表示在节点Xi处解的准确值y（Xi）的近似值。
设法构造序列yi所满足的一个方程（称为差分方程）
yi1yih（Xi,yi,h）（3-2）
作为求解公式，这是一个递推公式，从（Xo,yo）出发，采用步进方式，自左相右逐步算
出y（x）在所有节点x上的近似值yi（i1,2,,n）。
在公式（3-2）中，为求yi1只用到前面一步的值yi,这种方法称为单步法。
在公式（3-2）
中的yi1由yi明显表示出，称为显式公式。
而形如（3-3）
yi1yih（Xi,yi,yi1,h）（3-3）
的公式称为隐式公式，因为其右端中还包括y1。
如果由公式求y1时，不止用到前一个节点的值，则称为多步法。
由式（3-1）可得
（3-4）
dy=f（x,y）dx
两边在[x,Xii]上积分，得
xi1
y（Xi1）y（Xi）f（x,y（x））dx（3-5）
xi
由此可以看出，如果想构造求解公式，就要对右端的积分项作某种数值处理。
这种求解公式
的构造方法叫做数值积分法。
（二）一般的初值问题的解法
1.欧拉法和改进欧拉法
对于初值问题（3-1），采用数值积分法，从而得到（3-5）。
对于（3-5）右端的积分用
矩形公式（取左端点），则得到
xi1
f（x,y（x））dxhf（Xi,y（Xi））
Xi
进而得到（3-1）的求解公式（3-2）
yi1yihf（Xi,yi）（i=0,1,2,n-1）（3-6）
此公式称为欧拉（Euler）格式。
如果对式（3-5）右端的积分用梯形公式
为1h
f（x,y（x））dx（f（Xi,y（Xi））f（x1,y（Xi1）））
Xi2
则可以得到初值问题（3-1）的梯形求解公式如式（3-7）
h
yi1yif（Xi,yi）f（x1,yi1）（i=0,1,2,n-1）（3-7）
2
式（3-7）是个隐式公式。
可以采取先用欧拉格式求一个y（xi1）的初步近似值，记作yi1,
称之为预报值，然后用预报值yi1替代式（3-7）右端的yi1，再计算得到yi1，称之为校正值，这样建立起来的预报-校正方法称为改进欧拉格式
y1yihf（Xi,yi）
h-（3-8）
yi1yif（Xi,yi）f（x1,yi1）
2
2.龙格一库塔方法
在单步法中，应用最广泛的是龙格—库塔（Runge-kutta）法，简称R—K法。
下面直接
给出一种四阶的龙格一库塔法的计算公式（3-9）
1
yi丄（Ki2K22K36
hf（Xi,yJ
它也称为标准（古典）龙格-库塔法。
例3-1研究下列微分方程的初值问题
2y
1
y2
1X
y（0）0
解:
这是一个特殊的微分方程，其解的解析式可以给出，为
应用龙格—库塔法，取h=,根据式（3-9）编写一段程序，由零开始自左相右逐步算出y（x）
在所有节点Xi上的近似值yi。
计算结果见表3-1。
计算结果表明，四阶龙格—库塔方法的精度是较高的。
表3-1
Xn
yn
y（Xn）yn
0.
0.
0.
0.
实际上，MATLAB为常微分方程提供了很好的解题指令，使得求解常微分方程变得很容
易，并且能将问题及解答表现在图形上。
因此，我们可以不用根据式（3-9）编写较复杂的
程序，而只需应用MATLAB供的常微分方程解题器来解决问题。
下面给出用MATLAB编写的
解题程序。
functiondy=myfun（x,y）
dy=zeros（1,1）;
dy=1/（1+xA2）-2*yA2;
再编写利用解题器指令求解y的程序。
clear
x0=0;
fori=1:
4
xm=2*i;
y0=0;
[x,y]=ode45（'myfun',[x0xm],[yO]）;
formatlong
y（length（y））
end
plot（x,y,'-'）
运行上述程序，在得到几个点的函数值的同时，也得到函数y的曲线，如图3-1所示。
x
图3-1根据运算结果画出y的曲线
二简单电力系统的暂态稳定性
（一）物理过程分析
某简单电力系统如图3-2（a）所示，正常运行时发电机经过变压器和双回线路向无限大
系统供电。
发电机用电势E作为其等值电势，则电势E与无限大系统间的电抗为
这时发电机发出的电磁功率可表示为
势E与无限大系统之间的联系电抗为
Xl
（XdXti）（XT2）
2
在故障情况下发电机输出的电磁功率为
在短路故障发生之后，线路继电保护装置将迅速断开故障线路两端的断路器，如图
T2U=c
EjXdjXT1
-jXLjXT2U
jXL
jx
（b）
L
U
GT1
发电机输出的功率为
EjXdjXT1jXLjXT2
（c）
图3-2简单电力系统及其等值电路
（a）正常运行方式及其等值电路；（b）故障情况及其等值电路；（c）故障切除后及其等值电路
Po。
因此，可
于Po。
假定不计故障后几秒种之内调速器的作用，即认为机械功率始终保持
以得到此简单电力系统正常运行、故障期间及故障切除后的功率特性曲线如图3-3所示。
图3-3简单系统正常运行、故障期间及故障切除后的功率特性曲线
对于上述简单电力系统，我们可以根据等面积定则求得极限切除角。
但是，实际工作需
要知道在多少时间之内切除故障线路，也就是要知道与极限切除角对应的极限切除时间。
要
解决这个问题，必须求解发电机的转子运动方程。
（二）求解发电机的转子运动方程
求解发电机转子运动方程可以得出3-t和3-t的关系曲线。
其中3-t曲线一般称为摇
摆曲线。
在上述简单电力系统中故障期间的转子运动方程为
1）
速度，即1=2f=，其单位为弧度/秒；Tj――发电机的惯性时间常数，其单位为秒；Pt、
Pm——分别为机械和电磁功率，标幺值。
这是两个一阶的非线性常微分方程，它的起始条件是已知的，即
故障切除后，由于系统参数改变，以致发电机功率特性发生变化，必须开始求解另一组微分方程:
-
（1）1
d1（3-17）
（Ptpmsin）
dtTj
式中变量含义冋前述，其中
PM也为标幺值。
这组方程的起始条件为
t=tc；=c；=c
其中tc为给定的切除时间；c、c为与tc时刻对应的和，它们可由故障期间的3-t和
3-t的关系曲线求得（和都是不突变的）。
一般来说，在计算故障发生后几秒种的过程中，如果3始终不超过1800,而且振荡幅值越来越小，则系统是暂态稳定的。
当发电机与无限大系统之间发生振荡或失去同步时，在发电机的转子回路中，特别是阻
尼绕组中将有感应电流而形成阻尼转矩（也称为异步转矩）。
当作微小振荡时，阻尼功率可
表达为：
Pd=D=D
（1）（3-18）
式中，D称为阻尼功率系数；为转子角速度的偏移量，标幺值；为转子角速度，标幺
值。
阻尼功率系数d除了与发电机的参数有关外，还和原始功角、的振荡频率有关。
在
般情况下它是正数。
在原始功角较小，或者定子回路中有串联电容使定子回路总电阻相对
dtd
1）
（3-19）
dt
TJ[Pt
D（
1）Pmsin]
于总电抗较大时，d可能为负数。
如果考虑阻尼功率的影响，则故障后的转子运动方程又可表达为
电力系统暂态稳定计算包括两类问题，一类是应用数值计算法得出故障期间的曲线后,
根据曲线找到与极限切除角对应的极限切除时间，此时只需要求解微分方程（3-16）;另
类是已知故障切除时间，需要求出摇摆曲线来判断系统的稳定性，此时需要分段分别求解微
分方程（3-16）和（3-17）。
如果考虑阻尼转矩的影响，则此时需要分段分别求解微分方程
（3-16）和（3-19）。
三、例题
例3-2某简单电力系统如图3-4所示，取基准值Sb=220MVAUb=209KV。
换算后的参数
已经标在图中，其中一回线的电抗Xl=,Tj=秒。
设电力线路某一回的始端发生两相接地短
路。
假定E=常数。
（1）计算保持暂态稳定而要求的极限切除角。
（2）计算极限切除时间,
并且作出在秒切除故障时的3-t曲线。
图3-4某简单电力系统的接线图
解：
计算系统正常运行方式，决定E和0。
由3-3（a）的正序网络可得，此时系统的总电抗为
x=+++=
（1.00.798）2
1.0
0.798
=0
1.00.20.798
0.20.798）2
1.0
1
发电机的暂态电势为:
（2）故障后的功率特性
又由3-3（b）的负序、零序网络可得故障点的负序、零序等值电抗为
=（0.4320.138）（0.2430.122）
X2=（0.4320.138）（0.2430.122）
0.138（0.9720.122）_
Xo==
0.138（0.9720.122）
所以在正序网络故障点上的附加电抗为:
x
于是故障时等值电路如图
°222°1230.079
0.2220.123
3-3（c）所示，则
0.4330.365
x
0.433
0.365
2.80
0.079
因此，故障期间发电机的最大功率为:
（3）故障切除后的功率特性
此时最大功率为
0.295
0.1380.4860.122
1.041
EU
竺卫1.35
1.041
1800sin1型
1.35
1322°
j0.138j0.2434j0.122
I'
（c）
计算极限切除角
f（0）
（b）U=1.0
1.41j0.295j0.138j0.243j0.122
n_/"vv%
j0.079
U=1.0
1.41j0.295j0.138j0.486j0.122
:
-r^r-^r\
（d）
图3-5例題7-12的等值电路
（a）正常运行等值电路;
（b）
负序和零序等值电路；（C）故障时等值电路;
（d）故障切除后等值电路
COScm
PT（0h）PMCOShPMCOS0
PmPm
0
0.504cos34.53
0
「°面（13223453）「35cos132.2
1.350.504
cm
62.740
（5）
找出极限切除时间tcm
根据
（3-16），首先计算初值
34.53
0.6027,
1.0
令y
（1）=,y
（2）=。
编写描述故障期间转子运动方程的ODE文件，文件名为/myequ'。
functiondy=myequ（t,y）
dy=zeros（2,1）;
f=50;w1=2*pi*f;
dy
（1）=（y
（2）-1）*w1;
dy
（2）=（1/*再编写利用解题器指令求解y的程序。
clear
t0=0;tm=;
d0=180）*pi;w0=1;
[T,Y]=ode45（'myequ',[t0tm],[dOwO]）;
plot（T,（Y（:
1）/pi）*180,'-',,,'*'）
text,60,'\delta_{cmax}=\circ','FontSize',10）
text,56,'t_{cmax}=','FontSize',10）
图3-6例題7-12的S-t曲线
图3-6给出短路发生后0秒到秒期间的S-t计算曲线，根据最大切除角cm（62.740）
找到极限切除时间tcm为秒。
由图3-6可见，如果故障切除时间大于秒，则发电机的功角将不断地增大，最终失去暂态稳定。
在极限切除时间之前切除故障，发电机的摇摆曲线的状况
将在下面作计算、分析。
（6）不考虑阻尼转矩的影响，当故障切除时间为秒时通过计算得出S-t曲线
首先编写描述故障期间转子运动方程的ODE文件，文件名为”myfun01
文件名为”myfun02”。
'-',tc,（dc/pi）*180,'*'）
f=50;w1=2*pi*f;
TJ=;Pt=;P2m=;
dy=zeros（2,1）;
dy
（1）=（y
（2）-1）*w1;
dy
（2）=（1/TJ）*（Pt-P2m*sin（y
（1）））;
再编写描述故障切除后转子运动方程的ODE文件,
functiondy=myfun02（t,y）
f=50;w1=2*pi*f;
TJ=;Pt=;P3m=;
dy=zeros（2,1）;
dy
（1）=（y
（2）-1）*w1;
dy
（2）=（1/TJ）*（Pt-P3m*sin（y
（1）））;
编写利用解题器指令求解y的小程序。
clear
t0=0;tc=;tm=;
d0=180）*pi;w0=;
[T1,Y1]=ode45（'myfun01',[t0tc],[d0w0]）;
dc=Y1（length（Y1）,1）;
wc=Y1（length（Y1）,2）;
[T2,Y2]=ode45（'myfun02',[tctm],[dcwc]）;plot（T1,（Y1（:
1）/pi）*180,'-',T2,（Y2（:
1）/pi）*180,text,50,'\it{t}_{c}=','FontSize',8）
text,43,'\it{\delta}_{c}=\circ','FontSize',8）xlabel（'\it{t}'）
ylabel（'\it{\delta}'）
计算结果表明，功角沿着故障切除后的功角特性曲线根据等面积定则作等幅振荡，
最终运行在
图3-7所示。
实际上，由于阻尼转矩的影响，振荡的幅度是逐渐衰减的，功角
图3-7不考虑阻尼转矩影响，当秒切除故障时发电机的t曲线
k=o。
因此，发电机能够保持暂态稳定。
（7）考虑阻尼转矩的影响，当故障切除时间为秒时通过计算得出3-t曲线
描述故障期间转子运动方程的ODE文件与（6）相同，文件名也为”myfunOI”。
重新编写描述故障切除后转子运动方程的ODE文件，文件名为”myfun03”，阻尼功率系
数D取为15。
functiondy=myfunO3（t,y）
f=50;w1=2*pi*f;
TJ=;Pt=;P3m=;
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=;tm=4;
d0=*180;w0=1;
[T1,Y1]=ode45（'myfun01',[t0tc],[d0w0]）;
dc=Y1（length（Y1）,1）;
wc=Y1（length（Y1）,2）;
[T2,Y2]=ode45（'myfun03',[tctm],[dcwc]）;
plot（T1,（Y1（:
1）/pi）*180,'-',T2,（Y2（:
1）/pi）*180,'-',tc,（dc/pi）*180,'*'）
text,50,'\it{t}_c='，‘FontSize',8）
text,43,'\it{\delta}_{c}=\circ','FontSize',8）
xlabel（'\it{t}'）
ylabel（'\it{\delta}'）
t
图3-8不考虑阻尼转矩影响，当秒切除故障时发电机的t曲线
计算结果表明，功角沿着故障切除后的功角特性曲线作减幅振荡，如图3-8所示。
功角最终运行在k=o。
因此，发电机能够保持暂态稳定。
升级会员 