function[]=runga_kutta（TDdeltaXLstep）
E0=1.4239;
Xd=0.29;
Tj=15;
V0=1.0;
P0=1.0;
Q0=0.2;
XT1=0.13;
XT2=0.11;
XL=0.58;
h=0.01;
N=5/h;
XI=Xd+XT1+0.5*XL+XT2;????????????????????????????????%pre-fault
XII=XI+（Xd+XT1）*（0.5*XL+XT2）/deltaXL;??????????????????%during-fault
XIII=XI+0.5*XL;??????????????????????????????????????%post-fault
wn=2*pi*50;
Kw=wn^2/Tj;
P1m=E0*V0/XI;
P2m=E0*V0/XII;
P3m=E0*V0/XIII;
delta0=asin（P0/P1m）;????????????????????????????????
deltaB=pi-asin（P0/P3m）;
%--------------------------------------------------------------------------
%the?4th?order?R-K?method
delta3=zeros（1N）;
w3=zeros（1N）;
for?a=1:1:N/step
????k3=1;
????delta3（1）=delta0;
????w3（1）=wn;
????Pe3=P2m;?
?while?k3<=（N-1）
???if?（（h*k3>T）&&（h*（k3-1）<=T））?????????????????????%change??P2m??to??P3m??near?at?fault?clearing?time
???????Kdelta1=w3（k3）-wn;
???????Kw1=Kw/w3（k3）*（P0-D*（w3（k3）-wn）-Pe3*sin（delta3（k3）））;
???????deltaT1=delta3（k3）+0.5*（T-h*（k3-1））*Kdelta1;
???????wT1=w3（k3）+0.5*（T-h*（k3-1））*Kw1;
???????
???????Kdelta2=wT1-wn;
???????Kw2=Kw/wT1*（P0-D*（wT1-wn）-Pe3*sin（deltaT1））;
???????deltaT2=delta3（k3）+0.5*（T-h*（k3-1））*Kdelta2;
???????wT2=w3（k3）+0.5*（T-h*（k3-1））*Kw2;
???????
???????Kdelta3=wT2-wn;
???????Kw3=Kw/wT2*（P0-D*（wT2-wn）-Pe3*sin（deltaT2））;
???????deltaT3=delta3（k3）+（T-h*（k3-1））*Kdelta3;
???????wT3=w3（k3）+（T-h*（k3-1））*Kw3;
???????
???????Kdelta4=wT3-wn;
???????Kw4=Kw/wT3*（P0-D*（wT3-wn）-Pe3*sin（deltaT3））;
???????deltaC3=delta3（k3）+1/6*（T-h*（k3-1））*（Kdelta1+2*Kdelta2+2*Kdelta3+Kdelta4）;
???????wC3=w3（k3）+1/6*（T-h*（k3-1））*（Kw1+2*Kw2+2*Kw3+Kw4）;
???????deltalim（a）=deltaC3;
??????????
???????Pe3=P3m;
???????Kdelta1=wC3-wn;
???????Kw1=Kw/wC3*（P0-D*（wC3-wn）-Pe3*sin（deltaC3））;
???????deltaT1=deltaC3+0.5*（h*k3-T）*Kdelta1;
???????wT1=wC3+0.5*（h*k3-T）*Kw1;
???????
???????Kdelta2=wT1-wn;
???????Kw2=Kw/wT1*（P0-D*（wT1-wn）-Pe3*sin（deltaT1））;
???????deltaT2=deltaC3+0.5*（h*k3-T）*Kdelta2;
???????wT2=wC3+0.5*（h*k3-T）*Kw2;
???????
???????Kdelta3=wT2-wn;
???????Kw3=Kw/wT2*（P0-D*（wT2-wn）-Pe3*sin（deltaT2））;
???????deltaT3=deltaC3+（h*k3-T）*Kdelta3;
???????wT3=wC3+（h*k3-T）*Kw3;
???????
???????Kdelta4=wT3-wn;
???????Kw4=Kw/wT3*（P0-D*（wT3-wn）-Pe3*sin（deltaT3））;
???????delta3（k3+1）=deltaC3+1/6*（h*k3-T）*（Kdelta1+2*Kdelta2+2*Kdelta3+Kdelta4）;
???????w3（k3+1）=wC3+1/6*（h*k3-T）*（Kw1+2*Kw2+2*Kw3+Kw4）;
???????k3=k3+1;
???????continue
???end
???????Kdelta1=w3（k3）-wn;
???????Kw1=Kw/w3（k3）*（P0-D*（w3（k3）-wn）-Pe3*sin（delta3（k3）））;
???????deltaT1=delta3（k3）+0.5*h*Kdelta1;
???????wT1=w3（k3）+0.5*h*Kw1;
???????
???????Kdelta2=wT1-wn;
???????Kw2=Kw/wT1*（P0-D*（wT1-wn）-Pe3*sin（deltaT1））;
???????deltaT2=delta3（k3）+0.5*h*Kdelta2;
???????wT2=w3（k3）+0.5*h*Kw2;
???????
???????Kdelta3=wT2-wn;
???????Kw3=Kw/wT2*（P0-D*（wT2-wn）-Pe3*sin（deltaT2））;
???????deltaT3=delta3（k3）+h*Kdelta3;
???????wT3=w3（k3）+h*Kw3;
???????
???????Kdelta4=wT3-wn;
???????Kw4=Kw/wT3*（P0-D*（wT3-wn）-Pe3*sin（deltaT3））;
???????delta3（k3+1）=delta3（k3）+

?屬性????????????大小?????日期????時(shí)間???名稱
-----------?---------??----------?-----??----

?????文件?????244564??2014-10-16?15:17??陳文韜-暫穩(wěn)大作業(yè)\Report_.docx

?????文件???????3471??2014-01-05?15:58??陳文韜-暫穩(wěn)大作業(yè)\runga_kutta.m

?????文件???????7105??2014-01-05?15:58??陳文韜-暫穩(wěn)大作業(yè)\transient.m

?????文件???????7027??2014-01-05?16:41??陳文韜-暫穩(wěn)大作業(yè)\附件1_主程序transient?代碼.txt

?????文件???????3475??2014-01-04?16:07??陳文韜-暫穩(wěn)大作業(yè)\附件2_函數(shù)文件Runga_kutta?代碼.txt

?????目錄??????????0??2014-10-16?15:17??陳文韜-暫穩(wěn)大作業(yè)

-----------?---------??----------?-----??----

???????????????265642????????????????????6