电力系统稳定性分析matlab程序文档格式.docx

1、35.402335.699936.113036.639437.277038.023438.8763龙格35.221935.401635.698936.111636.637637.274738.020738.87310090.100.110.120.130.140.150.160.1739.364640.383541.504342.636643.777344.923046.070939.833240891842.005143.125644.250145.375746.499547.618748.730539.829540.887542.000243.120044.244045.369046.49

2、2347.611048.7224（1）欧拉法在matlaB中输入命令t,x,y,z=euler（doty,doty2doty3,0,5,0.1,0.01）可得t-w曲线，t-曲线分别如下图所示。具体功角，角速度数据分别见文件1.mat 和2.mat（2）欧拉改进法在matlab命令窗口输入t,x,y,z=reuler（,0,5,0.1,0.01）具体功角，角速度数据分别见文件3.mat 和4.mat（3）龙格库塔法在matlab命令窗口输入t,x,y,z=kunta（具体功角，角速度数据分别见文件5.mat 和6.mat2 运行Runge-Kutta,将参数阻尼D设置为0.05，不断更改参数切

3、除时间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】欧拉法主程序functiont

4、,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.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）=

5、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（ii+1）=y（ii）+h*feval（fun2,x（ii）,y（ii）; z（ii+1）=x（ii+1）*180/pi;endfor ii=n1+1:n y（ii+1）=y（ii）+h*feval（fun3,x（ii）,y（ii）;【reuler.m】改进欧拉法主程序：function t,x,y,z=reuler（fun1,fun2,fun3,t0,xf

6、inal,tm,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）; y（ii+1）=y（ii）+h/2*（g1+g2）; g1=feval（fun3,x（ii）,y（ii）; g2=feval（fun3,x0,y0）;subplot（1,2,1）plot（t,y）subplot（1,2,2）plot（t,z）;【kunta.m】龙格库塔法主程序func

7、tiont,x,y,z=kunta（fun1,fun2,fun3,t0,xfinal,tm,h） 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）

8、=x（ii）+h/6*（k1+2*k2+2*k3+k4）; y（ii+1）=y（ii）+h/6*（g1+2*g2+2*g3+g4）; g2=feval（fun3,x11,y11）; g3=feval（fun3,x22,y22）; g4=feval（fun3,x33,y33）;以下为子程序【doty.m】function fun1=doty（y）global 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 fun

9、3=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;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） D=0.2 修改t可

10、得到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）%初始公式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=E

11、0*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 K1w（i）=Kw*（P0-Pd-Peli3*sin（det3（i）*2*pi/360）/w3（i）; K2w（i）=Kw*（P0-Pd-Peli3*sin（det_c0*2*pi/360）/w_c0; K3w（i）=Kw*（P0-Pd-Peli3*sin（det_c1*2*pi/360）/w_c1; K4w（i）=Kw*（P0-Pd-Peli3*sin（det_c2*2*pi/360）/w_c2;%显示数据及图形t=0:0.0001:5; plot（t,det3（）,b）;%t=0: plot（t,w3（）,Welcome !欢迎您的下载，资料仅供参考！