电力系统稳定性分析matlab程序.pdf

《电力系统稳定性分析matlab程序.pdf》由会员分享，可在线阅读，更多相关《电力系统稳定性分析matlab程序.pdf（6页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、.-电力系统稳定性分析作业一电力系统稳定性分析作业一1euler.m,reuler.m,kunta.m 分别为1中的欧拉法，改良欧拉法，龙格库塔法的主程序；doty.m,doty2.m,doty3.m 均为1中子函数程序。Runge-Kutta.m 为2和3的运行程序。下表为三种方法的局部运行结果功角数据：时时0 0间间欧欧拉拉改改良良龙龙格格35.161535.161535.16150.010.0135.161535.222235.22190.100.100.020.0235.282835.402335.40160.110.110.030.0335.523635.699935.69890.1

2、20.1242.636643.125643.12000.040.0435.882036.113036.11160.130.1343.777344.250144.24400.050.0536.356036.639436.63760.140.1444.923045.375745.36900.060.0636.943437.277037.27470.150.1546.070946.499546.49230.070.0737.642238.023438.02070.160.1637.642247.618747.61100.080.0838.450038.876338.87310.170.1738.45

3、0048.730548.7224时时0 00909间间欧欧拉拉改改良良龙龙格格39.364639.833239.829540.383541.504340891842.005140.887542.00021欧拉法在 matlaB 中输入命令t,x,y,z=euler(doty,doty2,doty3,0,5,0.1,0.01)可得t-w 曲线，t-曲线分别如以下列图所示。具体功角，角速度数据分别见文件 1.mat 和 2.mat2欧拉改良法在 matlab 命令窗口输入t,x,y,z=reuler(doty,doty2,doty3,0,5,0.1,0.01)t-w 曲线，t-曲线分别如以下列图所

4、示。具体功角，角速度数据分别见文件 3.mat 和 4.mat3龙格库塔法在 matlab 命令窗口输入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 p

6、0 pp2 d1 w pp1f=50;Tj=11;p0=1.0;d1=0.05;xd=0.29;xt1=0.13;xt2=0.11;*=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)/*;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(i

7、i+1)=x(ii)+h*feval(fun1,y(ii);y(ii+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=(x

8、final-t0)/h;n1=(tm-t0)/h;global Kw p0 pp2 d1 w pp1f=50;Tj=11;p0=1.0;d1=0.05;xd=0.29;xt1=0.13;xt2=0.11;*=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)/*;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

9、(1)*180/pi;for ii=1:n1 t(ii+1)=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),

10、y(ii);x0=x(ii)+h*k1;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;glob

11、al Kw p0 pp2 d1 w pp1f=50;Tj=11;p0=1.0;d1=0.05;xd=0.29;xt1=0.13;xt2=0.11;*=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)/*;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

12、)=t(ii)+h;k1=feval(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(fun1,y22);g

14、3=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 wfun1=(y-w);【

15、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 xdTjxt1xt2xlDxzU0P0Q0;E0=1.4239;xd=0.29;Tj=11;xt1=0.13;xt2=0.11;xl=0.58;U0=1;P0=1;Q

16、0=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=0.05修改t可得到t=CCT=0.3111 det_c_lim=det3(3111)

17、%初始公式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*U0/Xdnum3;%四阶龙格库塔.优选-.-当t0.1s后去除故障for i=T+1:50000 K1det

18、(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);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-Peli3*sin(det_c0*2*pi/360)/w_c0;det_c1=det3(i)+0.5*h*K2det(i);w_c1=w3(i)+0.5*h*K2w(i);K3det(i)=(w_c

19、1-wn)*180/pi;Pd=D*(w_c1-wn);K3w(i)=Kw*(P0-Pd-Peli3*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-Peli3*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:0.0001:5;plot(t,det3(),b);%t=0:0.0001:5;plot(t,w3(),b);end.优选-