电力系统稳定性分析matlab程序.docx
《电力系统稳定性分析matlab程序.docx》由会员分享，可在线阅读，更多相关《电力系统稳定性分析matlab程序.docx（14页珍藏版）》请在冰点文库上搜索。
电力系统稳定性分析matlab程序
电力系统稳定性分析作业一
1
euler.m,reuler.m,kunta.m分别为
（1）中的欧拉法，改进欧拉法，龙格库塔法的主程序；doty.m,doty2.m,doty3.m均为
（1）中子函数程序。
Runge-Kutta.m为
（2）和（3）的运行程序。
下表为三种方法的部分运行结果功角数据：
时间
0
0.01
0.02
0.03
0.04
0.05
0.06
0.07
0.08
欧拉
35.1615
35.1615
35.2828
35.5236
35.8820
36.3560
36.9434
37.6422
38.4500
改进
35.1615
35.2222
35.4023
35.6999
36.1130
36.6394
37.2770
38.0234
38.8763
龙格
35.1615
35.2219
35.4016
35.6989
36.1116
36.6376
37.2747
38.0207
38.8731
时间
0．09
0.10
0.11
0.12
0.13
0.14
0.15
0.16
0.17
欧拉
39.3646
40.3835
41.5043
42.6366
43.7773
44.9230
46.0709
37.6422
38.4500
改进
39.8332
40．8918
42.0051
43.1256
44.2501
45.3757
46.4995
47.6187
48.7305
龙格
39.8295
40.8875
42.0002
43.1200
44.2440
45.3690
46.4923
47.6110
48.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）龙格库塔法
在matlab命令窗口输入[t,x,y,z]=kunta（'doty','doty2','doty3',0,5,0.1,0.01）
t-w曲线，t-δ曲线分别如下图所示。
具体功角，角速度数据分别见文件5.mat和6.mat
2
运行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;
globalKwp0pp2d1wpp1
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=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:
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;
end
forii=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;
globalKwp0pp2d1wpp1
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=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:
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）;
end
forii=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;
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）;
end
subplot（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;
globalKwp0pp2d1wpp1
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=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:
n1
t（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）;
end
forii=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）;
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）;
end
subplot（1,2,1）
plot（t,y）
subplot（1,2,2）
plot（t,z）;
以下为子程序【doty.m】
functionfun1=doty（y）
globalw
fun1=（y-w）;
【doty2.m】
functionfun2=doty2（x,y）
globalwp0pp1d1Kw
fun2=Kw*（p0-pp1*sin（x）-d1*（y-w））/y
【doty3.m】
functionfun3=doty3（x,y）
globalKwp0pp2d1w
fun3=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;
Xdnum2=Xdnum1+（xd+xt1）*（0.5*xl+xt2）/xdet;Peli2=E0*U0/Xdnum2;
Xdnum3=Xdnum1+0.5*xl;Peli3=E0*U0/Xdnum3;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%四阶龙格库塔
当t<0.1s时，故障
fori=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（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后清除故障
fori=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）;
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_c1-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
升级会员 