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

1、.实验三 电力系统暂态稳定分析电力系统暂态稳定计算实际上就是求解发电机转子运动方程的初值问题，从而得出-t和-t的关系曲线。每台发电机的转子运动方程是两个一阶非线性的常微分方程。因此，首先介绍常微分方程的初值问题的数值解法。一、 常微分方程的初值问题（一）问题及求解公式的构造方法我们讨论形如式（3-1）的一阶微分方程的初值问题 （3-1）设初值问题（3-1）的解为，为了求其数值解而采取离散化方法，在求解区间上取一组节点称（）为步长。在等步长的情况下，步长为用表示在节点处解的准确值的近似值。设法构造序列所满足的一个方程（称为差分方程） （3-2）作为求解公式，这是一个递推公式，从（，）出发，采用

2、步进方式，自左相右逐步算出在所有节点上的近似值（）。在公式（3-2）中，为求只用到前面一步的值，这种方法称为单步法。在公式（3-2）中的由明显表示出，称为显式公式。而形如（3-3） （3-3）的公式称为隐式公式，因为其右端中还包括。如果由公式求时，不止用到前一个节点的值，则称为多步法。由式（3-1）可得= （3-4）两边在，上积分，得 （3-5）由此可以看出，如果想构造求解公式，就要对右端的积分项作某种数值处理。这种求解公式的构造方法叫做数值积分法。（二）一般的初值问题的解法1 欧拉法和改进欧拉法对于初值问题（3-1），采用数值积分法，从而得到（3-5）。对于（3-5）右端的积分用矩形公式(取

3、左端点)，则得到进而得到（3-1）的求解公式（3-2） （=0，1，2，n-1） （3-6）此公式称为欧拉（Euler）格式。如果对式（3-5）右端的积分用梯形公式则可以得到初值问题（3-1）的梯形求解公式如式（3-7） （=0，1，2，n-1） （3-7）式（3-7）是个隐式公式。可以采取先用欧拉格式求一个的初步近似值，记作，称之为预报值，然后用预报值替代式（3-7）右端的，再计算得到，称之为校正值，这样建立起来的预报校正方法称为改进欧拉格式 （3-8）2 龙格库塔方法在单步法中，应用最广泛的是龙格库塔（Runge-kutta）法，简称RK法。下面直接给出一种四阶的龙格库塔法的计算公式（3-

4、9） （3-9）它也称为标准（古典）龙格库塔法。例3-1 研究下列微分方程的初值问题解：这是一个特殊的微分方程，其解的解析式可以给出，为应用龙格库塔法，取=0.25，根据式（3-9）编写一段程序，由零开始自左相右逐步算出在所有节点上的近似值。计算结果见表3-1。计算结果表明，四阶龙格库塔方法的精度是较高的。表3-12.00.399956994.3e-54.00.235291592.5e-66.00.162161793.7e-78.00.123076839.2e-8实际上，MATLAB为常微分方程提供了很好的解题指令，使得求解常微分方程变得很容易，并且能将问题及解答表现在图形上。因此，我们可以不

5、用根据式（3-9）编写较复杂的程序，而只需应用MATLAB提供的常微分方程解题器来解决问题。下面给出用MATLAB编写的解题程序。首先编写描述常微分方程的ODE文件，文件名为myfun，便于解题器调用它。function dy = myfun(x,y)dy = zeros(1,1);dy=1/(1+x2)-2*y2;再编写利用解题器指令求解y的程序。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的曲

6、线，如图3-1所示。图3-1 根据运算结果画出y的曲线二、 简单电力系统的暂态稳定性（一）物理过程分析某简单电力系统如图3-2(a)所示，正常运行时发电机经过变压器和双回线路向无限大系统供电。发电机用电势作为其等值电势，则电势与无限大系统间的电抗为 （3-10）这时发电机发出的电磁功率可表示为 （3-11）如果突然在一回输电线路始端发生不对称短路，如图3-2(b)所示。故障期间发电机电势与无限大系统之间的联系电抗为 （3-12）在故障情况下发电机输出的电磁功率为 （3-13）在短路故障发生之后，线路继电保护装置将迅速断开故障线路两端的断路器，如图3-2(c)所示。此时发电机电势与无限大系统间的

7、联系电抗为 （3-14）发电机输出的功率为 （3-15）图3-2 简单电力系统及其等值电路（a）正常运行方式及其等值电路；（b）故障情况及其等值电路；（c）故障切除后及其等值电路如果正常时发电机向无限大系统输送的有功功率为，则原动机输出的机械功率等于。假定不计故障后几秒种之内调速器的作用，即认为机械功率始终保持。因此，可以得到此简单电力系统正常运行、故障期间及故障切除后的功率特性曲线如图3-3所示。图3-3 简单系统正常运行、故障期间及故障切除后的功率特性曲线对于上述简单电力系统，我们可以根据等面积定则求得极限切除角。但是，实际工作需要知道在多少时间之内切除故障线路，也就是要知道与极限切除角对

8、应的极限切除时间。要解决这个问题，必须求解发电机的转子运动方程。（二）求解发电机的转子运动方程求解发电机转子运动方程可以得出-t和-t的关系曲线。其中-t曲线一般称为摇摆曲线。在上述简单电力系统中故障期间的转子运动方程为 （3-16）式中，功率角，其单位为弧度；转子角速度，标幺值；转子的同步角速度，即=314.16，其单位为弧度/秒；发电机的惯性时间常数，其单位为秒；、分别为机械和电磁功率，标幺值。这是两个一阶的非线性常微分方程，它的起始条件是已知的，即=0； =1.0；=故障切除后，由于系统参数改变，以致发电机功率特性发生变化，必须开始求解另一组微分方程： （3-17）式中变量含义同前述，其

9、中也为标幺值。这组方程的起始条件为=；=；=其中为给定的切除时间；、为与时刻对应的和，它们可由故障期间的-t和-t的关系曲线求得（和都是不突变的）。一般来说，在计算故障发生后几秒种的过程中，如果始终不超过180，而且振荡幅值越来越小，则系统是暂态稳定的。当发电机与无限大系统之间发生振荡或失去同步时，在发电机的转子回路中，特别是阻尼绕组中将有感应电流而形成阻尼转矩（也称为异步转矩）。当作微小振荡时，阻尼功率可表达为：= （3-18）式中，称为阻尼功率系数；为转子角速度的偏移量，标幺值；为转子角速度，标幺值。阻尼功率系数除了与发电机的参数有关外，还和原始功角、的振荡频率有关。在一般情况下它是正数。

10、在原始功角较小，或者定子回路中有串联电容使定子回路总电阻相对于总电抗较大时，D可能为负数。如果考虑阻尼功率的影响，则故障后的转子运动方程又可表达为 （3-19）电力系统暂态稳定计算包括两类问题，一类是应用数值计算法得出故障期间的曲线后，根据曲线找到与极限切除角对应的极限切除时间，此时只需要求解微分方程（3-16）；另一类是已知故障切除时间，需要求出摇摆曲线来判断系统的稳定性，此时需要分段分别求解微分方程（3-16）和（3-17）。如果考虑阻尼转矩的影响，则此时需要分段分别求解微分方程（3-16）和（3-19）。三、 例题例3-2 某简单电力系统如图3-4所示，取基准值=220MVA，=209K

11、V。换算后的参数已经标在图中，其中一回线的电抗=0.486，=8.18秒。设电力线路某一回的始端发生两相接地短路。假定=常数。（1）计算保持暂态稳定而要求的极限切除角。（2）计算极限切除时间，并且作出在0.15秒切除故障时的-t曲线。图3-4 某简单电力系统的接线图解：计算系统正常运行方式，决定和。由3-3(a)的正序网络可得，此时系统的总电抗为=0.295+0.138+0.243+0.122=0.798发电机的暂态电势为：=1.41=34.53（2）故障后的功率特性又由3-3(b)的负序、零序网络可得故障点的负序、零序等值电抗为=0.222=0.123所以在正序网络故障点上的附加电抗为：于是

12、故障时等值电路如图3-3(c)所示，则因此，故障期间发电机的最大功率为：（3）故障切除后的功率特性故障切除后的等值电路如图3-3(d)所示此时最大功率为图3-5 例題7-12的等值电路（a）正常运行等值电路；（b）负序和零序等值电路；（c）故障时等值电路；（d）故障切除后等值电路（4）计算极限切除角=0.458（5）找出极限切除时间根据（3-16），首先计算初值，令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)*w

13、1;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(0.194,56,t_cmax=0.194s,FontSize,10)图3-6 例題7-12的-t曲线图3-6给出短路发生后0秒到0.25秒期间的-t计算曲线，根

14、据最大切除角（）找到极限切除时间为0.194秒。由图3-6可见，如果故障切除时间大于0.194秒，则发电机的功角将不断地增大，最终失去暂态稳定。在极限切除时间之前切除故障，发电机的摇摆曲线的状况将在下面作计算、分析。（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/T

15、J)*(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的小程序。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

16、=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)text(0.28,43,itdelta_c=51.71circ,FontSize,8)xlabel(itt)ylabel(itdelta)计算结果表明，功角沿着故障切除后的功角特性曲线根据等面积定则作等幅振荡，如图3-7所示。实际上，由于阻尼转矩的

17、影响，振荡的幅度是逐渐衰减的，功角最终运行在=47.8。因此，发电机能够保持暂态稳定。图3-7 不考虑阻尼转矩影响，当0.15秒切除故障时发电机的曲线（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的小程序。cleart0=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