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

《MATLAB实验 电力系统暂态稳定分析.doc》由会员分享，可在线阅读，更多相关《MATLAB实验 电力系统暂态稳定分析.doc（12页珍藏版）》请在得力文库 - 分享文档赚钱的网站上搜索。

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-9

4、） （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，=209KV

11、。换算后的参数已经标在图中，其中一回线的电抗=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)*w1

13、;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/TJ

15、)*(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

18、(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 tm,dc wc);plot(T1,(Y1(:,1)/pi)*180,-,T2,(Y2(:,1)/pi)*180,-,tc,(dc/pi)*180,*)text

19、(0.3,50,itt_c=0.15s,FontSize,8)text(0.28,43,itdelta_c=51.71circ,FontSize,8)xlabel(itt)ylabel(itdelta)图3-8 不考虑阻尼转矩影响，当0.15秒切除故障时发电机的曲线计算结果表明，功角沿着故障切除后的功角特性曲线作减幅振荡，如图3-8所示。功角最终运行在=47.8。因此，发电机能够保持暂态稳定。- 30 - 定态持够，因. 在角。所如振减性功除障沿明曲机发故切.响转阻 ( ( , ,. , ,.= , ,* , ), ,-0 /(,() , 0 ( ,),) ),) )00 0 ( , 0 ;

20、. . = 0 序小解指解用)( - ( ( *-( )( = 0 . ,0 为取功，0 件，文 动转切述”0 ”名，）与 动子期曲-得计秒.0切故响转虑曲机时除切0响影阻 定稳保够，。 行角的渐是荡振的尼于上。 ，幅作积据曲功除切角， ( , ( = 0 .0)/ , ), -0* ( ) , ), ( ) )0 0, 0 *0/ 0 序程的令题) - )/() ( ) ( = .= 0= ;. 0 0 ”0 ”为件 的运后障描) *-( ) * )(), = 0= * ,0 ” 为， 的运间故描曲-过通 间除故响矩虑析析计在状线摇电故前之限极态去，大断角电， 0时除如见 。 0时切到角除根

21、曲 秒.0后短给曲-的- 00 . _ 0, , . _ ,.( 0,*)/, ) , = *0 . 0 序的解令解) *.0. =) *)( ) *=0) , 名， 动转障描写)(=初初先 -时切出 除切计电电切故 电障 ；等和负 值运正电等 -例率功所)(如值后除特率除故为率最电期，示(-图等时为电的点网 . .抗电零负点得络零) 特特障 . . 为势态 .= . 抗电系，络正)和定方常系线线系简 线 障秒.出且间限极（除极求稳暂保 数定短相发端回力设 .电回其图经参后 ， 0准取 统单 例） 方微别要时响矩阻果） 和 -分微段要，稳的判曲求要间切是另 -方解需时，切应对限与据后的障得计用

22、类，问包计态 达可方运后障，的虑如数可时较于阻路子电串路子小功始正它般在率频角原，有数电与率尼幺度子为值，偏角；系尼称 为为可，振作当步称矩阻而电将绕尼别中转发时去或发间大与的的态则，值荡且 过始如中的几障故，一）突都得线关间由它和对为间切给；件始程方幺也中同量 -（ 程方一解须化生率功致改参于除=0 ，的件起的分常非一值值功电机位，间惯机秒为位 =即角的幺标角转度位其 程方运的中统简在线摇一 线关 得方转机方运子电程程转电求必个决间除的角限与要，障切之多道需实但。限求定根们，系单曲性功除故障行运单 示 曲性后障故障、统力单以，持保功为用器内种几不假等械出原为功送统大向常电电及障 （等及障（；

23、其式运电值统系 - 率功 -（ 抗电间大与电时此 -如路的路开迅置护路，发 率磁的机下 - 抗系的系大势机障示) 图路对端线回一 示可率的电00（ 抗电间无势则等作电电供大向回和过经发运所( 如电分程定定的曲 果算 示示图线曲到得的函点得， ( )( ) , 0 : 0 序序求器解 - ), 它调解， 名，文 程微写序序写编 下问解解分的 用只序杂写 根用可此。现答题且易得变微解令解很提分为 - .- .- .- .-的较是方阶，表计 见计。的有在逐左自，序编-（据.=塔，出式解解方殊个解问问方列究 法库龙）标 -公计库龙四一接下 简） （库龙泛应中方塔 （ 格欧为法正的来样正之称算，） 代替

24、用，预之记似的求拉先采式隐是- -（-，0= 如解梯 值得公梯分右 式） 为 - - = 式的-（得，点左形分的 （。）（得法积采-题值拉拉解的初般法积值叫造公这处种作的对，公构果出以 -，积， 可 法多则值个前不，公括还右因公为称 （。公为表由的-公在为称，步前用，）（）值似节在步相，进采发（式公是，解 程方称程一满造值值值解节在长长下步在步节节区求方化而数了为的-问 - 问初方阶一 如方方公及问问初的法值题初方分绍此。方性非一程动转发线关- 出从问初动子发求上计态分分态力-0 分上发初出-动一。绍题值问公 方初-的了方节下长节造程程 解，（进相值）前，在的公 称右括不则法 积 构的作公积般

25、的值采（分形（式=- 分公 =（-隐采的记用 算样的法 塔龙（）简四库 ）龙究问殊解式塔（序自在计 ，是的 为令变得现此 杂用的解 编序 名，解 -器求 点得线示果曲程如所经回大电则间 （电可 回端图障势系 机磁 发护路如此与大 （功 -统电；障（电电送原出不内功，单障后 单故曲，们。需多障要角个必程运机 线摇统中 转幺角 为秒，电值常件 0参功生解方 中始；间和由线都）几如始值则的间时中绕将称作可 称；偏子尼电数频般它小电子较可，运达 计问计障与限时解另间求的，分 果矩要 ）例统准， 图回 发短保求除限.线简线和络电 .势为.特零得零 电时-期率除除值所率-电负；电切电 切-初)描转名 , =*) 0)解的0 * , , , . -的后秒曲除时。见0，断去限故线析虑除间-故的为 ” * = () *)后 ”0 0 .= ( )令序 0 0 ( (* ,)0 ( ( 角除据，上的荡角 够保影切时虑故秒得期动） 述动， 0 0 =) * 用解 0 , ,(,)( ,), = ( .故明沿性振。因，定 - 30 -