MATLAB实验电力系统暂态稳定分析.docx

1、MATLAB实验电力系统暂态稳定分析实验三电力系统暂态稳定分析电力系统暂态稳定计算实际上就是求解发电机转子运动方程的初值问题， 从而得出3 -t和3 -t的关系曲线。每台发电机的转子运动方程是两个一阶非线性的常微分方程。因此， 首先介绍常微分方程的初值问题的数值解法。-、常微分方程的初值问题(一)问题及求解公式的构造方法我们讨论形如式(3-1 )的一阶微分方程的初值问题(3-1 )y (x) f(x, y), a x b y(xo) yo设初值问题(3-1 )的解为y(x)，为了求其数值解而采取离散化方法，在求解区间 a,b上取一组节点a xo x1 xi xi 1 xn b称hi xi 1

2、xi ( i 0,1, ,n 1 )为步长。在等步长的情况下，步长为b ahn用yi表示在节点Xi处解的准确值y(Xi)的近似值。设法构造序列yi所满足的一个方程(称为差分方程)yi 1 yi h (Xi,yi,h) (3-2 )作为求解公式，这是一个递推公式，从( Xo, yo)出发，采用步进方式，自左相右逐步算出y(x)在所有节点x上的近似值yi ( i 1,2, , n )。在公式(3-2 )中，为求yi 1只用到前面一步的值 yi ,这种方法称为单步法。在公式(3-2 )中的yi 1由yi明显表示出，称为显式公式。而形如( 3-3)yi 1 yi h (Xi, yi, yi 1, h)

3、 (3-3 )的公式称为隐式公式，因为其右端 中还包括y 1。如果由公式求y 1时，不止用到前一个节点的值，则称为多步法。由式(3-1 )可得(3-4 )dy = f (x, y)dx两边在x , Xi i 上积分，得xi 1y(Xi 1) y(Xi) f (x, y(x)dx (3-5 )xi由此可以看出，如果想构造求解公式，就要对右端的积分项作某种数值处理。 这种求解公式的构造方法叫做数值积分法。(二)一般的初值问题的解法1.欧拉法和改进欧拉法对于初值问题(3-1 )，采用数值积分法，从而得到( 3-5 )。对于(3-5 )右端的积分用矩形公式(取左端点)，则得到xi 1f (x, y(x

4、)dx h f(Xi,y(Xi)Xi进而得到(3-1 )的求解公式(3-2 )yi 1 yi h f (Xi, yi) ( i =0, 1, 2, n-1 ) ( 3-6)此公式称为欧拉(Euler )格式。如果对式(3-5 )右端的积分用梯形公式为1 hf (x, y(x)dx (f(Xi ,y(Xi) f(x 1,y(Xi 1)Xi 2则可以得到初值问题(3-1 )的梯形求解公式如式(3-7 )hyi 1 yi f (Xi, yi) f (x 1, yi 1) (i =0, 1, 2, n-1 ) (3-7 )2式(3-7 )是个隐式公式。可以采取先用欧拉格式求一个 y(xi 1)的初步近

5、似值，记作 yi 1 ,称之为预报值，然后用预报值 yi 1替代式(3-7 )右端的yi 1，再计算得到yi 1，称之为校正 值，这样建立起来的预报-校正方法称为改进欧拉格式y 1 yi h f (Xi, yi)h - (3-8 )yi 1 yi f (Xi, yi) f (x 1, yi 1)22.龙格一库塔方法在单步法中，应用最广泛的是龙格库塔（ Runge-kutta ）法，简称R K法。下面直接给出一种四阶的龙格一库塔法的计算公式（ 3-9）1yi 丄(Ki 2K2 2K3 6h f(Xi,yJ它也称为标准（古典）龙格-库塔法。例3-1 研究下列微分方程的初值问题2y1y 21 Xy(

6、0) 0解:这是一个特殊的微分方程，其解的解析式可以给出，为应用龙格库塔法，取 h=,根据式（3-9）编写一段程序，由零开始自左相右逐步算出 y（x）在所有节点Xi上的近似值yi。计算结果见表3-1。计算结果表明，四阶龙格库塔方法的精 度是较高的。表3-1Xnyny(Xn ) yn0.0.0.0.实际上，MATLAB为常微分方程提供了很好的解题指令，使得求解常微分方程变得很容易，并且能将问题及解答表现在图形上。因此，我们可以不用根据式（ 3-9 ）编写较复杂的程序，而只需应用MATLAB供的常微分方程解题器来解决问题。 下面给出用MATLAB编写的解题程序。fun cti on dy = my

7、fun( x,y)dy = zeros(1,1);dy=1/(1+xA2)-2*yA2;再编写利用解题器指令求解 y的程序。clearx0=0;for i=1:4xm=2*i;y0=0;x,y = ode45( myfun ,x0 xm,yO);format longy(le ngth(y)endplot(x,y,-)运行上述程序，在得到几个点的函数值的同时，也得到函数 y的曲线，如图3-1所示。x图3-1 根据运算结果画出y的曲线二 简单电力系统的暂态稳定性(一)物理过程分析某简单电力系统如图 3-2(a)所示，正常运行时发电机经过变压器和双回线路向无限大系统供电。发电机用电势 E作为其等值

8、电势，则电势 E与无限大系统间的电抗为这时发电机发出的电磁功率可表示为势E与无限大系统之间的联系电抗为Xl (Xd Xti )( XT2)2在故障情况下发电机输出的电磁功率为在短路故障发生之后，线路继电保护装置将迅速断开故障线路两端的断路器，如图T2 U= cE jXd jXT1-jXL jXT2 U jXLjx(b)LUG T1发电机输出的功率为E jXd jXT1 jXL jXT2(c)图3-2 简单电力系统及其等值电路(a)正常运行方式及其等值电路； (b)故障情况及其等值电路；(c)故障切除后及其等值电路Po。因此，可于Po。假定不计故障后几秒种之内调速器的作用，即认为机械功率始终保持

9、以得到此简单电力系统正常运行、故障期间及故障切除后的功率特性曲线如图 3-3所示。图3-3 简单系统正常运行、故障期间及故障切除后的功率特性曲线对于上述简单电力系统，我们可以根据等面积定则求得极限切除角。 但是，实际工作需要知道在多少时间之内切除故障线路， 也就是要知道与极限切除角对应的极限切除时间。 要解决这个问题，必须求解发电机的转子运动方程。（二）求解发电机的转子运动方程求解发电机转子运动方程可以得出3 -t和3 -t的关系曲线。其中3 -t曲线一般称为摇摆曲线。在上述简单电力系统中故障期间的转子运动方程为1)速度，即1=2 f =，其单位为弧度/秒；Tj 发电机的惯性时间常数，其单位为

10、秒； Pt、Pm 分别为机械和电磁功率，标幺值。这是两个一阶的非线性常微分方程，它的起始条件是已知的，即故障切除后，由于系统参数改变， 以致发电机功率特性发生变化， 必须开始求解另一组 微分方程:-(1) 1d 1 (3-17 )(Pt p m sin )dt Tj式中变量含义冋前述，其中P M也为标幺值。这组方程的起始条件为t =tc ； = c ； = c其中tc为给定的切除时间； c、 c为与tc时刻对应的 和，它们可由故障期间的3 -t和3 -t的关系曲线求得（ 和 都是不突变的）。一般来说，在计算故障发生后几秒种的过程 中，如果3始终不超过 180 0 ,而且振荡幅值越来越小，则系统

11、是暂态稳定的。当发电机与无限大系统之间发生振荡或失去同步时， 在发电机的转子回路中， 特别是阻尼绕组中将有感应电流而形成阻尼转矩（也称为异步转矩） 。当作微小振荡时，阻尼功率可表达为：Pd = D =D（ 1） （3-18）式中，D称为阻尼功率系数； 为转子角速度的偏移量，标幺值； 为转子角速度，标幺值。阻尼功率系数 d除了与发电机的参数有关外，还和原始功角、 的振荡频率有关。在般情况下它是正数。 在原始功角较小，或者定子回路中有串联电容使定子回路总电阻相对dt d1)(3-19 )dtTJPtD(1) P m sin 于总电抗较大时，d可能为负数。如果考虑阻尼功率的影响，则故障后的转子运动方

12、程又可 表达为电力系统暂态稳定计算包括两类问题，一类是应用数值计算法得出故障期间的曲线后,根据曲线找到与极限切除角对应的极限切除时间，此时只需要求解微分方程（ 3-16 ）;另类是已知故障切除时间，需要求出摇摆曲线来判断系统的稳定性， 此时需要分段分别求解微分方程（3-16 ）和（3-17 ）。如果考虑阻尼转矩的影响，则此时需要分段分别求解微分方程（3-16 ）和（3-19 ）。三、例题例3-2 某简单电力系统如图3-4所示，取基准值Sb=220MVA Ub=209KV。换算后的参数已经标在图中，其中一回线的电抗 Xl=, Tj =秒。设电力线路某一回的始端发生两相接地短路。假定E =常数。(

13、1)计算保持暂态稳定而要求的极限切除角。 (2)计算极限切除时间,并且作出在秒切除故障时的3 -t曲线。图3-4某简单电力系统的接线图解：计算系统正常运行方式，决定 E和0。由3-3(a)的正序网络可得，此时系统的总电抗为x =+=(1.0 0.798)21.00.798=01.0 0.2 0.7980.2 0.798)21.01发电机的暂态电势为:(2)故障后的功率特性又由3-3(b)的负序、零序网络可得故障点的负序、零序等值电抗为=(0.432 0.138) (0.243 0.122)X2 = (0.432 0.138) (0.243 0.122)0.138(0.972 0.122) _X

14、o = =0.138 (0.972 0.122)所以在正序网络故障点上的附加电抗为:x于是故障时等值电路如图222 123 0.0790.222 0.1233-3(c)所示，则0.433 0.365 x0.4330.3652.800.079因此，故障期间发电机的最大功率为:(3 )故障切除后的功率特性此时最大功率为0.2950.138 0.486 0.1221.041EU竺卫1.351.0411800 sin1 型1.351322j 0.138 j0.243 4 j 0.122I(c)计算极限切除角f(0)(b) U=1.01.41 j 0.295 j 0.138 j 0.243 j 0.12

15、2n_/vv% j 0.079U=1.01.41j 0.295 j 0.138 j 0.486 j 0.122:- rr-r (d)图3-5 例題7-12的等值电路(a)正常运行等值电路;(b)负序和零序等值电路；(C)故障时等值电路;(d)故障切除后等值电路COS cmPT( 0 h) P M COS h P M COS 0Pm P m00.504 cos34.530 面(1322 3453)35cos132.21.35 0.504cm62.740(5)找出极限切除时间tcm根据(3-16)，首先计算初值34.530.6027,1.0令y(1)= , y(2)=。编写描述故障期间转子运动方程

16、的 ODE文件，文件名为/ myequ。fun cti on dy = myequ(t,y)dy = zeros(2,1);f=50;w1=2*pi*f;dy(1) = (y(2)-1)*w1;dy(2) = (1/* 再编写利用解题器指令求解 y的程序。cleart0=0;tm=;d0=180)*pi;w0=1;T,Y = ode45( myequ ,t0 tm,dO wO);plot(T,(Y(:,1)/pi)*180, - ,*)text,60, delta_cmax=circ , FontSize ,10)text,56, t_cmax= , FontSize ,10)图3-6 例題7

17、-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 =

18、zeros(2,1);dy(1) = (y(2)-1)*w1;dy(2) = (1/TJ)*(Pt-P2m*sin(y(1);再编写描述故障切除后转子运动方程的 ODE文件,function dy = 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的小程序。cleart0=0; tc=; tm=;d0=180)*pi; w0=;T1,Y1 = ode45( myfun01 ,t0 tc,d0 w

19、0);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, text,50, itt_c= , FontSize ,8)text,43, itdelta_c=circ , FontSize ,8) xlabel( itt )ylabel( itdelta )计算结果表明， 功角 沿着故障切除后的功角特性曲线根据等面积定则作等幅振荡，最终运行在图3-7所示。实际上，由于阻尼转矩的影响，振荡的幅度是

20、逐渐衰减的，功角图3-7不考虑阻尼转矩影响，当秒切除故障时发电机的 t曲线k=o。因此，发电机能够保持暂态稳定。(7)考虑阻尼转矩的影响，当故障切除时间为秒时通过计算得出3 -t曲线描述故障期间转子运动方程的 ODE文件与(6)相同，文件名也为” myfunOI ”。重新编写描述故障切除后转子运动方程的 ODE文件，文件名为” myfun03”，阻尼功率系数D取为15。fun cti on dy = myfunO 3(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

21、)*(Pt-D*(y (2)-1)-P3m*si n(y(1);再编写利用解题器指令求解 y的小程序。cleart0=0; tc=; tm=4;d0=*180; w0=1;T1,Y1 = ode45( myfu n01 ,t0 tc,d0 w0);dc=Y1(le ngth(Y1),1);wc=Y1(le ngth(Y1),2);T2,Y2 = ode45( myfu n03 ,tc tm,dc wc);plot(T1,(Y1(:,1)/pi)*180, - ,T2,(Y2(:,1)/pi)*180, - ,tc,(dc/pi)*180, *)text,50, itt_c= ，Fo ntSize ,8)text,43, itdelta_c=circ , FontSize ,8)xlabel( itt )ylabel( itdelta )t图3-8 不考虑阻尼转矩影响，当秒切除故障时发电机的 t曲线计算结果表明，功角 沿着故障切除后的功角特性曲线作减幅振荡， 如图3-8所示。功 角 最终运行在 k=o。因此，发电机能够保持暂态稳定。