电力系统暂态稳定计算

运用Matlab变成，对单机无穷大系统三相短路后的暂态行为进行计算，分析不同切除时间对暂态稳定的影响，比较不同算法，如欧拉法（显隐式）、龙格库塔法的对计算性能的影响。

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）+

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件     244564  2014-10-16 15:17  陈文韬-暂稳大作业\Report_.docx

     文件       3471  2014-01-05 15:58  陈文韬-暂稳大作业\runga_kutta.m

     文件       7105  2014-01-05 15:58  陈文韬-暂稳大作业\transient.m

     文件       7027  2014-01-05 16:41  陈文韬-暂稳大作业\附件1_主程序transient 代码.txt

     文件       3475  2014-01-04 16:07  陈文韬-暂稳大作业\附件2_函数文件Runga_kutta 代码.txt

     目录          0  2014-10-16 15:17  陈文韬-暂稳大作业

----------- ---------  ---------- -----  ----

               265642                    6