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

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

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

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

4、积分用 矩形公式 （取左端点 ） ，则得到xi 1 f（x,y（x）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 ）右端的积分用梯形公式1xi xh(x,y(x)dx h2 (f(xi,y(xi) f(xi 1,y(xi 1)则可以得到初值问题（ 3-1 ）的梯形求解公式如式（ 3-7 ）yi 1 yi h2 f（xi,yi） f（xi 1,yi1） （i =0，1， 2，n-1 ） （3-7）式（ 3-7 ）是

5、个隐式公式。可以采取先用欧拉格式求一个 y（xi 1） 的初步近似值，记作 yi 1 ，称之为预报值，然后用预报值 yi 1替代式（ 3-7 ）右端的 yi 1 ，再计算得到 yi 1 ，称之为校正 值，这样建立起来的预报校正方法称为改进欧拉格式3-8)yi 1 yi h f（xi, yi）yi 1 yi h f（xi,yi） f（xi 1,yi 1） 22 龙格库塔方法在单步法中，应用最广泛的是龙格库塔（ Runge-kutta ）法，简称 RK 法。下面直接给出一种四阶的龙格库塔法的计算公式（ 3-9 ）1yi 1 yi (K1 2K2 2K 3 K4)6K1 h f (xi ,yi)3-

6、9)h1K2 h f(xi , yi K1)22h1K3 h f (xi ,yi K2)22K4 h f(xi h,yi K3)它也称为标准(古典)龙格库塔法。例 3-1 研究下列微分方程的初值问题12 y 2 2y1 x2 y(0) 0解：这是一个特殊的微分方程，其解的解析式可以给出，为1x应用龙格库塔法，取 h =0.25 ，根据式( 3-9 )编写一段程序，由零开始自左相右逐步算出y(x)在所有节点 xi 上的近似值 yi 。计算结果见表 3-1 。计算结果表明， 四阶龙格库塔方法 的精度是较高的。表 3-1xnyny(xn ) yn2.00.399956994.3e-54.00.235

7、291592.5e-66.00.162161793.7e-78.00.123076839.2e-8实际上， MATLAB为常微分方程提供了很好的解题指令，使得求解常微分方程变得很容易，并且能将问题及解答表现在图形上。因此，我们可以不用根据式( 3-9 )编写较复杂的程序，而只需应用 MATLAB提供的常微分方程解题器来解决问题。 下面给出用 MATLAB编写的解题程序。首先编写描述常微分方程的 ODE文件，文件名为 myfun，便于解题器调用它。 function dy = myfun(x,y)dy = zeros(1,1);dy=1/(1+x2)-2*y2;再编写利用解题器指令求解 y 的程

8、序。clearx0=0;for i=1:4xm=2*i;y0=0;x,y = ode45( myfun ,x0 xm,y0); format longy(length(y)endplot(x,y, - )运行上述程序，在得到几个点的函数值的同时，也得到函数 y 的曲线，如图 3-1 所示。简单电力系统的暂态稳定性一)物理过程分析这时发电机发出的电磁功率可表示为势E 与无限大系统之间的联系电抗为在故障情况下发电机输出的电磁功率为在短路故障发生之后，线路继电保护装置将迅速断开故障线路两端的断路器，如图3-2(c) 所示。此时发电机电势 E 与无限大系统间的联系电抗为发电机输出的功率为(c)图 3-

9、2 简单电力系统及其等值电路P0 。因此，可3-3 所示。于 P0 。假定不计故障后几秒种之内调速器的作用，即认为机械功率始终保持 以得到此简单电力系统正常运行、故障期间及故障切除后的功率特性曲线如图PP图 3-3 简单系统正常运行、故障期间及故障切除后的功率特性曲线PPT P00 k c m但是，实际工作需对于上述简单电力系统， 我们可以根据等面积定则求得极限切除角。要知道在多少时间之内切除故障线路， 也就是要知道与极限切除角对应的极限切除时间。解决这个问题，必须求解发电机的转子运动方程。二）求解发电机的转子运动方程求解发电机转子运动方程可以得出 -t 和 -t 的关系曲线。 其中 -t 曲

10、线一般称为摇摆曲线。在上述简单电力系统中故障期间的转子运动方程为d( 1) 1ddt 1 (P P sin ) (3-16 )ddt T1J (PT PM sin )式中， 功率角，其单位为弧度； 转子角速度，标幺值； 1转子的同步角速度，即 1=2 f =314.16 ，其单位为弧度 /秒；TJ 发电机的惯性时间常数， 其单位为秒；PT 、 P M 分别为机械和电磁功率，标幺值。这是两个一阶的非线性常微分方程，它的起始条件是已知的，即故障切除后， 由于系统参数改变， 以致发电机功率特性发生变化， 必须开始求解另一组 微分方程：d ( 1) 1dt (3-17 )d1dt T (PT P M

11、sin )式中变量含义同前述，其中 P M 也为标幺值。这组方程的起始条件为t =tc； = c ； = c其中 tc 为给定的切除时间； c、 c 为与 t c时刻对应的 和 ，它们可由故障期间的 -t 和 -t 的关系曲线求得( 和 都是不突变的) 。一般来说，在计算故障发生后几秒种的过程 中，如果始终不超过 180o ，而且振荡幅值越来越小，则系统是暂态稳定的。当发电机与无限大系统之间发生振荡或失去同步时， 在发电机的转子回路中， 特别是阻 尼绕组中将有感应电流而形成阻尼转矩(也称为异步转矩) 。当作微小振荡时，阻尼功率可 表达为：PD = D = D( 1) (3-18 )式中， D

12、称为阻尼功率系数； 为转子角速度的偏移量，标幺值； 为转子角速度，标幺值。阻尼功率系数 D 除了与发电机的参数有关外，还和原始功角、 的振荡频率有关。在一般情况下它是正数。 在原始功角较小， 或者定子回路中有串联电容使定子回路总电阻相对 于总电抗较大时， D可能为负数。如果考虑阻尼功率的影响，则故障后的转子运动方程又可 表达为ddt ( 1) 1d 1 PT D( 1) P M sin dt TJ3-19 )电力系统暂态稳定计算包括两类问题，一类是应用数值计算法得出故障期间的曲线后， 根据曲线找到与极限切除角对应的极限切除时间，此时只需要求解微分方程（ 3-16 ）；另一 类是已知故障切除时间

13、， 需要求出摇摆曲线来判断系统的稳定性， 此时需要分段分别求解微 分方程（ 3-16 ）和（ 3-17 ）。如果考虑阻尼转矩的影响，则此时需要分段分别求解微分方程 （3-16 ）和（ 3-19 ）。三、 例题例 3-2 某简单电力系统如图 3-4 所示，取基准值 SB =220MVA，U B =209KV。换算后的参数 已经标在图中，其中一回线的电抗 xL =0.486 ， TJ =8.18 秒。设电力线路某一回的始端发生 两相接地短路。假定 E =常数。（ 1）计算保持暂态稳定而要求的极限切除角。 （ 2）计算极限 切除时间，并且作出在 0.15 秒切除故障时的 -t 曲线。U =1.0图

14、3-4 某简单电力系统的接线图解：计算系统正常运行方式，决定 E 和 0 。由 3-3（a） 的正序网络可得，此时系统的总电抗为x =0.295+0.138+0.243+0.122=0.798发电机的暂态电势为：E = (1.0 0.2 0.798)2 (1.0 0.798)2 =1.412）故障后的功率特性又由 3-3（b） 的负序、零序网络可得故障点的负序、零序等值电抗为(0.432 0.138) (0.243 0.122)= =0.222(0.432 0.138) (0.243 0.122)x = 0.138(0.972 0.122) =0.123 x0 = 0.138 (0.972 0

15、.122)所以在正序网络故障点上的附加电抗为：0.079x 0.222 0.123x 0.222 0.123于是故障时等值电路如图 3-3（c） 所示，则因此，故障期间发电机的最大功率为：3）故障切除后的功率特性故障切除后的等值电路如图 3-3（d） 所示x 0.295 0.138 0.486 0.122 1.041此时最大功率为j 0.138 j 0 .243 4 j 0.122(b)f(0)U =1.0E 1.41 j 0.295 j 0.138 j 0.243 j 0.122U =1.0E 1.41j 0.295 j 0.138 j 0.486 j 0.122(d)故障切除后等值电路图

16、3-5 例題 7-12 的等值电路（a）正常运行等值电路； （ b）负序和零序等值电路； （ c）故障时等值电路；4）计算极限切除角cos cmPT ( 0 h ) P M cos h P M cos 0001.0 (132.2 34.53) 1.35cos132.20 0.504 cos34.53 0 1801.35 0.504=0.458cm 62.740(5)找出极限切除时间 t cm根据( 3-16 )，首先计算初值34.530 0.6027 ， 0 1.00 180 0令 y(1)= ， y(2)= 。编写描述故障期间转子运动方程的 ODE文件，文件名为 myequ。function

17、 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 的程序。cleart0=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.76circ , FontSize ,10) text(

18、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”。function d

19、y = 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

20、);编写利用解题器指令求解 y 的小程序。cleart0=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, itt_c=0.15s , FontSize ,8)

21、text(0.28,43, itdelta_c=51.71circ , FontSize ,8)xlabel( itt )ylabel( itdelta )计算结果表明， 功角 沿着故障切除后的功角特性曲线根据等面积定则作等幅振荡， 如 图 3-7 所示。实际上，由于阻尼转矩的影响，振荡的幅度是逐渐衰减的，功角 最终运行在 k =47.8 o 。 因此，发电机能够保持暂态稳定。7）考虑阻尼转矩的影响，当故障切除时间为 0.15 秒时通过计算得出 -t 曲线描述故障期间转子运动方程的 ODE文件与（ 6）相同，文件名也为” myfun01 ”。 重新编写描述故障切除后转子运动方程的 ODE文件，

22、文件名为” 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);text(0.28,43, itdelta_c=51.71circ xlabel( itt ) ylabel( itdelta )t图 3-8 不考虑阻尼转矩影响，当 0.15 秒切除故障时发电机的 t 曲线 计算结果表明， 功角 沿着故障切除后的功角特性曲线作减幅振荡， 如图 3-8 所示。 功 角 最终运行在 k =47.8 o 。因此，发电机能够保持暂态稳定。