MATLAB实验电力系统暂态稳定分析报告.docx
《MATLAB实验电力系统暂态稳定分析报告.docx》由会员分享，可在线阅读，更多相关《MATLAB实验电力系统暂态稳定分析报告.docx（17页珍藏版）》请在冰豆网上搜索。
MATLAB实验电力系统暂态稳定分析报告
实验三电力系统暂态稳定分析
电力系统暂态稳定计算实际上就是求解发电机转子运动方程的初值问题，从而得出δ-t和ω-t的关系曲线。
每台发电机的转子运动方程是两个一阶非线性的常微分方程。
因此，首先介绍常微分方程的初值问题的数值解法。
、常微分方程的初值问题
一）问题及求解公式的构造方法
我们讨论形如式（3-1）的一阶微分方程的初值问题
3-1）
y（x）f（x,y）,axby（x0）y0
设初值问题（3-1）的解为y（x），为了求其数值解而采取离散化方法，在求解区间[a,b]上
取一组节点
ax0x1xixi1xnb
称hixi1xi（i0,1,,n1）为步长。
在等步长的情况下，步长为
ba
h
n
用yi表示在节点xi处解的准确值y（xi）的近似值。
设法构造序列yi所满足的一个方程（称为差分方程）
yi1yih（xi,yi,h）（3-2）
作为求解公式，这是一个递推公式，从（x0，y0）出发，采用步进方式，自左相右逐步算
出y（x）在所有节点xi上的近似值yi（i1,2,,n）。
在公式（3-2）中，为求yi1只用到前面一步的值yi，这种方法称为单步法。
在公式（3-2）中的yi1由yi明显表示出，称为显式公式。
而形如（3-3）
yi1yih（xi,yi,yi1,h）（3-3）
的公式称为隐式公式，因为其右端中还包括yi1。
如果由公式求yi1时，不止用到前一个节点的值，则称为多步法。
由式（3-1）可得
dy=f（x,y）dx（3-4）
两边在[xi，xi1]上积分，得
xi1
y（xi1）y（xi）f（x,y（x））dx（3-5）
xi
由此可以看出，如果想构造求解公式，就要对右端的积分项作某种数值处理。
这种求解公式
的构造方法叫做数值积分法。
（二）一般的初值问题的解法
1．欧拉法和改进欧拉法
对于初值问题（3-1），采用数值积分法，从而得到（3-5）。
对于（3-5）右端的积分用矩形公式（取左端点），则得到
xi1f（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）右端的积分用梯形公式
1
xix
h
（x,y（x））dxh2（f（xi,y（xi））f（xi1,y（xi1）））
则可以得到初值问题（3-1）的梯形求解公式如式（3-7）
yi1yih2f（xi,yi）f（xi1,yi1）（i=0，1，2，n-1）（3-7）
式（3-7）是个隐式公式。
可以采取先用欧拉格式求一个y（xi1）的初步近似值，记作yi1，
称之为预报值，然后用预报值yi1替代式（3-7）右端的yi1，再计算得到yi1，称之为校正值，这样建立起来的预报－校正方法称为改进欧拉格式
3-8）
yi1yihf（xi,yi）
yi1yihf（xi,yi）f（xi1,yi1）2
2．龙格—库塔方法
在单步法中，应用最广泛的是龙格－库塔（Runge-kutta）法，简称R－K法。
下面直接
给出一种四阶的龙格－库塔法的计算公式（3-9）
1
yi1yi（K12K22K3K4）
6
K1hf（xi,yi）
3-9）
h1
K2hf（xi,yiK1）
22
h1
K3hf（xi,yiK2）
22
K4hf（xih,yiK3）
它也称为标准（古典）龙格－库塔法。
例3-1研究下列微分方程的初值问题
12y22y
1x2y（0）0
解：
这是一个特殊的微分方程，其解的解析式可以给出，为
1x
应用龙格－库塔法，取h=0.25，根据式（3-9）编写一段程序，由零开始自左相右逐步算出
y（x）在所有节点xi上的近似值yi。
计算结果见表3-1。
计算结果表明，四阶龙格－库塔方法的精度是较高的。
表3-1
xn
yn
y（xn）yn
2.0
0.39995699
4.3e-5
4.0
0.23529159
2.5e-6
6.0
0.16216179
3.7e-7
8.0
0.12307683
9.2e-8
实际上，MATLAB为常微分方程提供了很好的解题指令，使得求解常微分方程变得很容
易，并且能将问题及解答表现在图形上。
因此，我们可以不用根据式（3-9）编写较复杂的
程序，而只需应用MATLAB提供的常微分方程解题器来解决问题。
下面给出用MATLAB编写的
解题程序。
首先编写描述常微分方程的ODE文件，文件名为ˊmyfunˊ，便于解题器调用它。
functiondy=myfun（x,y）
dy=zeros（1,1）;
dy=1/（1+x^2）-2*y^2;
再编写利用解题器指令求解y的程序。
clear
x0=0;
fori=1:
4
xm=2*i;
y0=0;
[x,y]=ode45（'myfun',[x0xm],[y0]）;formatlong
y（length（y））
end
plot（x,y,'-'）
运行上述程序，在得到几个点的函数值的同时，也得到函数y的曲线，如图3-1所示。
简单电力系统的暂态稳定性
一）物理过程分析
这时发电机发出的电磁功率可表示为
势E与无限大系统之间的联系电抗为
在故障情况下发电机输出的电磁功率为
在短路故障发生之后，线路继电保护装置将迅速断开故障线路两端的断路器，如图
3-2（c）所示。
此时发电机电势E与无限大系统间的联系电抗为
发电机输出的功率为
（c）
图3-2简单电力系统及其等值电路
P0。
因此，可
3-3所示。
于P0。
假定不计故障后几秒种之内调速器的作用，即认为机械功率始终保持以得到此简单电力系统正常运行、故障期间及故障切除后的功率特性曲线如图
P
P
图3-3简单系统正常运行、故障期间及故障切除后的功率特性曲线
P
PTP0
0kcm
但是，实际工作需
对于上述简单电力系统，我们可以根据等面积定则求得极限切除角。
要知道在多少时间之内切除故障线路，也就是要知道与极限切除角对应的极限切除时间。
解决这个问题，必须求解发电机的转子运动方程。
二）求解发电机的转子运动方程
求解发电机转子运动方程可以得出δ-t和ω-t的关系曲线。
其中δ-t曲线一般称为摇
摆曲线。
在上述简单电力系统中故障期间的转子运动方程为
d
（1）1
ddt1（PPsin）（3-16）
ddtT1J（PTPMsin）
式中，——功率角，其单位为弧度；——转子角速度，标幺值；1——转子的同步角
速度，即1=2f=314.16，其单位为弧度/秒；TJ——发电机的惯性时间常数，其单位为秒；
PT、PM――分别为机械和电磁功率，标幺值。
这是两个一阶的非线性常微分方程，它的起始条件是已知的，即
故障切除后，由于系统参数改变，以致发电机功率特性发生变化，必须开始求解另一组微分方程：
d
（1）1
dt（3-17）
d1
dtT（PTPMsin）
式中变量含义同前述，其中PM也为标幺值。
这组方程的起始条件为
t=tc；=c；=c
其中tc为给定的切除时间；c、c为与tc时刻对应的和，它们可由故障期间的δ-t和ω-t的关系曲线求得（和都是不突变的）。
一般来说，在计算故障发生后几秒种的过程中，如果δ始终不超过180o，而且振荡幅值越来越小，则系统是暂态稳定的。
当发电机与无限大系统之间发生振荡或失去同步时，在发电机的转子回路中，特别是阻尼绕组中将有感应电流而形成阻尼转矩（也称为异步转矩）。
当作微小振荡时，阻尼功率可表达为：
PD=D=D
（1）（3-18）
式中，D称为阻尼功率系数；为转子角速度的偏移量，标幺值；为转子角速度，标幺
值。
阻尼功率系数D除了与发电机的参数有关外，还和原始功角、的振荡频率有关。
在
一般情况下它是正数。
在原始功角较小，或者定子回路中有串联电容使定子回路总电阻相对于总电抗较大时，D可能为负数。
如果考虑阻尼功率的影响，则故障后的转子运动方程又可表达为
ddt
（1）1
d1[PTD
（1）PMsin]dtTJ
3-19）
电力系统暂态稳定计算包括两类问题，一类是应用数值计算法得出故障期间的曲线后，根据曲线找到与极限切除角对应的极限切除时间，此时只需要求解微分方程（3-16）；另一类是已知故障切除时间，需要求出摇摆曲线来判断系统的稳定性，此时需要分段分别求解微分方程（3-16）和（3-17）。
如果考虑阻尼转矩的影响，则此时需要分段分别求解微分方程（3-16）和（3-19）。
三、例题
例3-2某简单电力系统如图3-4所示，取基准值SB=220MVA，UB=209KV。
换算后的参数已经标在图中，其中一回线的电抗xL=0.486，TJ=8.18秒。
设电力线路某一回的始端发生两相接地短路。
假定E=常数。
（1）计算保持暂态稳定而要求的极限切除角。
（2）计算极限切除时间，并且作出在0.15秒切除故障时的δ-t曲线。
U=1.0
图3-4某简单电力系统的接线图
解：
计算系统正常运行方式，决定E和0。
由3-3（a）的正序网络可得，此时系统的总电抗为
x=0.295+0.138+0.243+0.122=0.798
发电机的暂态电势为：
E=（1.00.20.798）2（1.00.798）2=1.41
2）故障后的功率特性
又由3-3（b）的负序、零序网络可得故障点的负序、零序等值电抗为
（0.4320.138）（0.2430.122）
==0.222
（0.4320.138）（0.2430.122）
x=0.138（0.9720.122）=0.123x0=0.138（0.9720.122）
所以在正序网络故障点上的附加电抗为：
0.079
x0.2220.123
x0.2220.123
于是故障时等值电路如图3-3（c）所示，则
因此，故障期间发电机的最大功率为：
3）故障切除后的功率特性
故障切除后的等值电路如图3-3（d）所示
x0.2950.1380.4860.1221.041
此时最大功率为
j0.138j0.2434j0.122
（b）
f（0）
U=1.0
E1.41j0.295j0.138j0.243j0.122
U=1.0
E1.41j0.295j0.138j0.486j0.122
（d）
故障切除后等值电路
图3-5例題7-12的等值电路
（a）正常运行等值电路；（b）负序和零序等值电路；（c）故障时等值电路；
4）计算极限切除角
coscm
PT（0h）PMcoshPMcos0
00
1.0（132.234.53）1.35cos132.200.504cos34.530180
1.350.504
=0.458
cm62.740
（5）找出极限切除时间tcm
根据（3-16），首先计算初值
34.53
00.6027，01.0
01800
令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/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',[t0tm],[d0w0]）;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给出短路发生后0秒到0.25秒期间的δ-t计算曲线，根据最大切除角cm
（62.740）找到极限切除时间tcm为0.194秒。
由图3-6可见，如果故障切除时间大于0.194秒，则发电机的功角将不断地增大，最终失去暂态稳定。
在极限切除时间之前切除故障，发电机的摇摆曲线的状况将在下面作计算、分析。
（6）不考虑阻尼转矩的影响，当故障切除时间为0.15秒时通过计算得出δ-t曲线
首先编写描述故障期间转子运动方程的ODE文件，文件名为”myfun01”。
functiondy=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”。
functiondy=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',[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,'-',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所示。
实际上，由于阻尼转矩的影响，振荡的幅度是逐渐衰减的，功角最终运行在k=47.8o。
因此，发电机能够保持暂态稳定。
7）考虑阻尼转矩的影响，当故障切除时间为0.15秒时通过计算得出δ-t曲线
描述故障期间转子运动方程的ODE文件与（6）相同，文件名也为”myfun01”。
重新编写描述故障切除后转子运动方程的ODE文件，文件名为”myfun03”，阻尼功率系
数D取为15。
functiondy=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）））;
text（0.28,43,'\it{\delta}_{c}=51.71\circ'xlabel（'\it{t}'）ylabel（'\it{\delta}'）
t
图3-8不考虑阻尼转矩影响，当0.15秒切除故障时发电机的t曲线计算结果表明，功角沿着故障切除后的功角特性曲线作减幅振荡，如图3-8所示。
功角最终运行在k=47.8o。
因此，发电机能够保持暂态稳定。
升级会员 