1、電力系統穩定性分析作業一1 euler.m ,reuler.m, kunta.m分別為（1）中的歐拉法，改進歐拉法，龍格庫塔法的主程序；doty.m,doty2.m,doty3.m均為（1）中子函數程序。Runge-Kutta.m為（2）和（3）的運行程序。下表為三種方法的部分運行結果功角數據：時間00.010.020.030.040.050.060.070.08歐拉35.161535.161535.282835.523635.882036.356036.943437.642238.4500改進35.161535.222235.402335.699936.113036.639437.27703

2、8.023438.8763龍格35.161535.221935.401635.698936.111636.637637.274738.020738.8731時間0090.100.110.120.130.140.150.160.17歐拉39.364640.383541.504342.636643.777344.923046.070937.642238.4500改進39.833240891842.005143.125644.250145.375746.499547.618748.7305龍格39.829540.887542.000243.120044.244045.369046.492347.61

3、1048.7224（1）歐拉法在matlaB中輸入命令t,x,y,z=euler('doty','doty2','doty3',0,5,0.1,0.01)可得t-w曲線，t-曲線分別如下圖所示。具體功角，角速度數據分別見文件1.mat 和2.mat（2）歐拉改進法在matlab命令窗口輸入t,x,y,z=reuler('doty','doty2','doty3',0,5,0.1,0.01)t-w曲線，t-曲線分別如下圖所示。具體功角，角速度數據分別見文件3.mat 和4.mat（3）龍格庫塔法在ma

4、tlab命令窗口輸入t,x,y,z=kunta('doty','doty2','doty3',0,5,0.1,0.01)t-w曲線，t-曲線分別如下圖所示。具體功角，角速度數據分別見文件5.mat 和6.mat2 運行Runge-Kutta,將參數阻尼D設置為0.05，不斷更改參數切除時間t的值，當t=0.2728和t=0.2730時，運行程序分別得到如下兩圖： 則當阻尼D=0.05 時，臨界切除時間 CCT=0.2729類似可以求得：阻尼D=0.2時，臨界切除時間為CCT=0.5729由以上數據我們可以看出：阻尼增大時，臨界切除時間也增大了。即

5、伴隨阻尼的增大，功角和角速度振蕩衰減更明顯，系統更容易回到平衡狀態，系統的穩定性更好。3接地阻抗X=0.05時，臨界切除時間CCT=0.2462接地阻抗X=0.1時，臨界切除時間CCT=0.3112由以上數據我們可以看出：接地阻抗增大時，系統臨界切除時間也增大了，系統穩定性變好。附注：以下為詳細的程序清單。【Euler.m】歐拉法主程序functiont,x,y,z=euler(fun1,fun2,fun3,t0,xfinal,tm,h)n=(xfinal-t0)/h;n1=(tm-t0)/h; global Kw p0 pp2 d1 w pp1 f=50;Tj=11;p0=1.0;d1=0.

6、05;xd=0.29;xt1=0.13;xt2=0.11;xx=0.07149;xl=0.58;E0=1.4239;V0=1.0;w=2*pi*f;Kw=w2/Tj;x1=xd+xt1+0.5*xl+xt2;x2=x1+(xd+xt1)*(0.5*xl+xt2)/xx;x3=x1+0.5*xl;pp1=E0*V0/x2;pp2=E0*V0/x3;t(1)=t0;x(1)=asin(p0*x1/E0/V0);y(1)=2*pi*f;z(1)=x(1)*180/pi;for ii=1:n1 t(ii+1)=t(ii)+h; x(ii+1)=x(ii)+h*feval(fun1,y(ii); y(i

7、i+1)=y(ii)+h*feval(fun2,x(ii),y(ii); z(ii+1)=x(ii+1)*180/pi;endfor ii=n1+1:n t(ii+1)=t(ii)+h; x(ii+1)=x(ii)+h*feval(fun1,y(ii); y(ii+1)=y(ii)+h*feval(fun3,x(ii),y(ii); z(ii+1)=x(ii+1)*180/pi;end【reuler.m】改進歐拉法主程序：function t,x,y,z=reuler(fun1,fun2,fun3,t0,xfinal,tm,h)n=(xfinal-t0)/h;n1=(tm-t0)/h;glob

8、al Kw p0 pp2 d1 w pp1 f=50;Tj=11;p0=1.0;d1=0.05;xd=0.29;xt1=0.13;xt2=0.11;xx=0.07149;xl=0.58;E0=1.4239;V0=1.0;w=2*pi*f;Kw=w2/Tj;x1=xd+xt1+0.5*xl+xt2;x2=x1+(xd+xt1)*(0.5*xl+xt2)/xx;x3=x1+0.5*xl;pp1=E0*V0/x2;pp2=E0*V0/x3;t(1)=t0;x(1)=asin(p0*x1/E0/V0);y(1)=2*pi*f;z(1)=x(1)*180/pi; for ii=1:n1 t(ii+1)=

9、t(ii)+h; k1=feval(fun1,y(ii); g1=feval(fun2,x(ii),y(ii); x0=x(ii)+h*k1; y0=y(ii)+h*g1; k2=feval(fun1,y0); g2=feval(fun2,x0,y0); x(ii+1)=x(ii)+h/2*(k1+k2); z(ii+1)=x(ii+1)*180/pi; y(ii+1)=y(ii)+h/2*(g1+g2);endfor ii=n1+1:n t(ii+1)=t(ii)+h; k1=feval(fun1,y(ii); g1=feval(fun3,x(ii),y(ii); x0=x(ii)+h*k1

10、; y0=y(ii)+h*g1; k2=feval(fun1,y0); g2=feval(fun3,x0,y0); x(ii+1)=x(ii)+h/2*(k1+k2); z(ii+1)=x(ii+1)*180/pi; y(ii+1)=y(ii)+h/2*(g1+g2);endsubplot(1,2,1)plot(t,y)subplot(1,2,2)plot(t,z);【kunta.m】龍格庫塔法主程序functiont,x,y,z=kunta(fun1,fun2,fun3,t0,xfinal,tm,h)n=(xfinal-t0)/h;n1=(tm-t0)/h;global Kw p0 pp2

11、d1 w pp1 f=50;Tj=11;p0=1.0;d1=0.05;xd=0.29;xt1=0.13;xt2=0.11;xx=0.07149;xl=0.58;E0=1.4239;V0=1.0;w=2*pi*f;Kw=w2/Tj;x1=xd+xt1+0.5*xl+xt2;x2=x1+(xd+xt1)*(0.5*xl+xt2)/xx;x3=x1+0.5*xl;pp1=E0*V0/x2;pp2=E0*V0/x3;t(1)=t0;x(1)=asin(p0*x1/E0/V0);y(1)=2*pi*f;z(1)=x(1)*180/pi; for ii=1:n1 t(ii+1)=t(ii)+h; k1=f

12、eval(fun1,y(ii); g1=feval(fun2,x(ii),y(ii); x11=x(ii)+0.5*h*k1; y11=y(ii)+0.5*h*g1; k2=feval(fun1,y11); g2=feval(fun2,x11,y11); x22=x(ii)+0.5*h*k2; y22=y(ii)+0.5*h*g2; k3=feval(fun1,y22); g3=feval(fun2,x22,y22); x33=x(ii)+h*k2; y33=y(ii)+h*g2; k4=feval(fun1,y33); g4=feval(fun2,x33,y33); x(ii+1)=x(ii

13、)+h/6*(k1+2*k2+2*k3+k4); z(ii+1)=x(ii+1)*180/pi; y(ii+1)=y(ii)+h/6*(g1+2*g2+2*g3+g4);endfor ii=n1+1:n t(ii+1)=t(ii)+h; k1=feval(fun1,y(ii); g1=feval(fun3,x(ii),y(ii); x11=x(ii)+0.5*h*k1; y11=y(ii)+0.5*h*g1; k2=feval(fun1,y11); g2=feval(fun3,x11,y11); x22=x(ii)+0.5*h*k2; y22=y(ii)+0.5*h*g2; k3=feval(

14、fun1,y22); g3=feval(fun3,x22,y22); x33=x(ii)+h*k2; y33=y(ii)+h*g2; k4=feval(fun1,y33); g4=feval(fun3,x33,y33); x(ii+1)=x(ii)+h/6*(k1+2*k2+2*k3+k4); z(ii+1)=x(ii+1)*180/pi; y(ii+1)=y(ii)+h/6*(g1+2*g2+2*g3+g4);endsubplot(1,2,1)plot(t,y)subplot(1,2,2)plot(t,z);以下為子程序【doty.m】function fun1=doty(y)global

15、wfun1=(y-w);【doty2.m】function fun2=doty2(x,y)global w p0 pp1 d1 Kwfun2=Kw*(p0-pp1*sin(x)-d1*(y-w)/y【doty3.m】function fun3=doty3(x,y)global Kw p0 pp2 d1 wfun3=Kw*(p0-pp2*sin(x)-d1*(y-w)/y;以下為四階龍格庫塔法求臨界切除時間程序：function jj % 初始值syms E0 xd Tj xt1 xt2 xl D xz U0 P0 Q0;E0=1.4239;xd=0.29;Tj=11;xt1=0.13;xt2=

16、0.11;xl=0.58;U0=1;P0=1;Q0=0.2; h=0.0001; %²參數調試xdet=0.07149; %xdet=0.07149;D=0.05; %D=0.05;t=0.2730;%xdet=0.07149; D=0.05 修改t可得到t= CCT=0.2729 det_c_lim=det3(2729)%xdet=0.07149; D=0.2 修改t可得到t= CCT=0.5729 det_c_lim=det3(5729)%xdet=0.05; D=0.05 修改t可得到t= CCT=0.2462 det_c_lim=det3(2462)%xdet=0.1; D=

17、0.05 修改t可得到t= CCT=0.3111 det_c_lim=det3(3111) %初始公式wn=2*pi*50; w1(1)=wn;w2(1)=wn;w3(1)=wn; T=t*10000;det1(1)=35.1615;det2(1)=35.1615;det3(1)=35.1615; Kw=wn2/Tj; Xdnum1=xd+xt1+0.5*xl+xt2; Peli1=E0*U0/Xdnum1; Xdnum2=Xdnum1+(xd+xt1)*(0.5*xl+xt2)/xdet; Peli2=E0*U0/Xdnum2; Xdnum3=Xdnum1+0.5*xl; Peli3=E0*

18、U0/Xdnum3; %四階龍格庫塔當t<0.1s時，故障for i=1:T K1det(i)=(w3(i)-wn)*180/pi; Pd=D*(w3(i)-wn); K1w(i)=Kw*(P0-Pd-Peli2*sin(det3(i)*2*pi/360)/w3(i); det_c0=det3(i)+0.5*h*K1det(i); w_c0=w3(i)+0.5*h*K1w(i); K2det(i)=(w_c0-wn)*180/pi; Pd=D*(w_c0-wn); K2w(i)=Kw*(P0-Pd-Peli2*sin(det_c0*2*pi/360)/w_c0; det_c1=det3(

19、i)+0.5*h*K2det(i); w_c1=w3(i)+0.5*h*K2w(i); K3det(i)=(w_c1-wn)*180/pi; Pd=D*(w_c1-wn); K3w(i)=Kw*(P0-Pd-Peli2*sin(det_c1*2*pi/360)/w_c1; det_c2=det3(i)+h*K3det(i); w_c2=w3(i)+h*K3w(i); K4det(i)=(w_c2-wn)*180/pi; Pd=D*(w_c2-wn); K4w(i)=Kw*(P0-Pd-Peli2*sin(det_c2*2*pi/360)/w_c2; det3(i+1)=det3(i)+h*(K1det(i)+ 2*K2det(i)+ 2*K3det(i)+ K4det(i)/6; w3(i+1)=w3(i)+h*(K1w(i)+2*K2w(i)+2*K3w(i)+K4w(i)/6;end% t>0.1s后清除故障for i=T+1:50000 K1det(i)=(w3(i)-wn)*180/pi; Pd=D*(w3(i)-wn); K1w(i)=Kw*(P0-Pd-Peli3*sin(det3(i)*2*pi/360)/w3(i);