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

1电力系统稳定性分析作业一1euler.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.277038.023438.8763龙格35.161535.221935.401635.698936.111636.637637.274738.020738.8731时间0．090.100.110.120.130.140.150.160.17欧拉39.364640.383541.504342.636643.777344.923046.070937.642238.4500改进39.833240．891842.005143.125644.250145.375746.499547.618748.7305龙格39.829540.887542.000243.120044.244045.369046.492347.611048.7224（1）欧拉法在matlaB中输入命令[t,x,y,z]=euler('doty','doty2','doty3',0,5,0.1,0.01)可得t-w曲线，t-δ曲线分别如下图所示。具体功角，角速度数据分别见文件1.mat和2.mat2（2）欧拉改进法在matlab命令窗口输入[t,x,y,z]=reuler('doty','doty2','doty3',0,5,0.1,0.01)t-w曲线，t-δ曲线分别如下图所示。具体功角，角速度数据分别见文件3.mat和4.mat（3）龙格库塔法在matlab命令窗口输入[t,x,y,z]=kunta('doty','doty2','doty3',0,5,0.1,0.01)t-w曲线，t-δ曲线分别如下图所示。具体功角，角速度数据分别见文件5.mat和6.mat32运行Runge-Kutta,将参数阻尼D设置为0.05，不断更改参数切除时间t的值，当t=0.2728和t=0.2730时，运行程序分别得到如下两图：则当阻尼D=0.05时，临界切除时间CCT=0.2729类似可以求得：阻尼D=0.2时，临界切除时间为CCT=0.5729由以上数据我们可以看出：阻尼增大时，临界切除时间也增大了。即伴随阻尼的增大，功角和角速度振荡衰减更明显，系统更容易回到平衡状态，系统的稳定性更好。3接地阻抗X=0.05时，临界切除时间CCT=0.2462接地阻抗X=0.1时，临界切除时间CCT=0.3112由以上数据我们可以看出：接地阻抗增大时，系统临界切除时间也增大了，系统稳定性变好。附注：以下为详细的程序清单。【Euler.m】欧拉法主程序function[t,x,y,z]=euler(fun1,fun2,fun3,t0,xfinal,tm,h)n=(xfinal-t0)/h;n1=(tm-t0)/h;4globalKwp0pp2d1wpp1f=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=w^2/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;forii=1:n1t(ii+1)=t(ii)+h;x(ii+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;endforii=n1+1:nt(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;globalKwp0pp2d1wpp1f=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=w^2/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;5forii=1:n1t(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);endforii=n1+1:nt(ii+1)=t(ii)+h;k1=feval(fun1,y(ii));g1=feval(fun3,x(ii),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】龙格库塔法主程序function[t,x,y,z]=kunta(fun1,fun2,fun3,t0,xfinal,tm,h)n=(xfinal-t0)/h;n1=(tm-t0)/h;globalKwp0pp2d1wpp1f=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=w^2/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;6pp2=E0*V0/x3;t(1)=t0;x(1)=asin(p0*x1/E0/V0);y(1)=2*pi*f;z(1)=x(1)*180/pi;forii=1:n1t(ii+1)=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)+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);endforii=n1+1:nt(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);g3=feval(fun3,x22,y22);x33=x(ii)+h*k2;y33=y(ii)+h*g2;k4=feval(fun1,y33);7g4=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】functionfun1=doty(y)globalwfun1=(y-w);【doty2.m】functionfun2=doty2(x,y)globalwp0pp1d1Kwfun2=Kw*(p0-pp1*sin(x)-d1*(y-w))/y【doty3.m】functionfun3=doty3(x,y)globalKwp0pp2d1wfun3=Kw*(p0-pp2*sin(x)-d1*(y-w))/y;以下为四阶龙格库塔法求临界切除时间程序：functionjj%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%初始值symsE0xdTjxt1xt2xlDxzU0P0Q0;E0=1.4239;xd=0.29;Tj=11;xt1=0.13;xt2=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.2729det_c_lim=det3(2729)%xdet=0.07149;D=0.2修改t可得到t=CCT=0.5729det_c_lim=det3(5729)%xdet=0.05;D=0.05修改t可得到t=CCT=0.2462det_c_lim=det3(2462)%xdet=0.1;D=0.05修改t可得到t=CCT=0.3111det_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=wn^2/Tj;Xdnum1=xd+xt1+0.5*xl+xt2;Peli1=E0*U0/Xdnum1;8Xdnum2=Xdnum1+(xd+xt1)*(0.5*x