三机九节点的暂态稳定性仿真程序

王锡凡版三机九节点的暂态稳定性仿真程序，包括计凸极效应还是不计凸极效应的仿真

clc;clear;
n1=9;n2=9;
B1=load（‘b1.txt‘‘ascii‘）;
B2=load（‘b2.txt‘‘ascii‘）;
m=0;r=0;
for i=1:n1
    if B1（i8）==1
        m=m+1;
    end
end
r=n1-m-1;
e=B1（1:n12）;f=B1（1:n13）;
ps=B1（1:n14）-B1（1:n16）;
qs=B1（1:n15）-B1（1:n17）;

Y=zeros（n1）;
for j=1:n2
      p=B2（j2）;q=B2（j3）;J=zeros（2*n1）;
      Y（pq）=-1./（（B2（j4）+1i*B2（j5））*B2（j7））;
      Y（qp）=Y（pq）;
      Y（pp）=Y（pp）+1/（B2（j4）+1i*B2（j5））+1i*B2（j6）;
      Y（qq）=Y（qq）+1/（（B2（j4）+1i*B2（j5））*B2（j7）^2）+1i*B2（j6）;
end
%导纳矩阵Y
H1=zeros（n1）;N1=zeros（n1）;J1=zeros（n1）;L1=zeros（n1）;
G=real（Y）;B=imag（Y）;
dp（1:（n1-1））=1;dq（1:m）=1;dx（1:（2*m+r））=1;
while （abs（max（dp）>=1.0e-5）|abs（max（dq）>=1.0e-5）|abs（max（dx（1:（m+r）））>=1.0e-5）|abs（max（dx（（m+r+1）:（2*m+r））.*e（1:m））>=1.0e-5））
for i=1:n1
    for j=1:n1
        si（ij）=sin（f（i）-f（j））;
        co（ij）=cos（f（i）-f（j））;
        if i==j
            H1（ij）=e（i）^2*B（ii）+qs（i）;
            N1（ij）=-e（i）^2*G（ii）-ps（i）;
            J1（ij）=e（i）^2*G（ii）-ps（i）;
            L1（ij）=e（i）^2*B（ii）-qs（i）;
        else
            H1（ij）=-e（i）*e（j）*（G（ij）*si（ij）-B（ij）*co（ij））; 
            N1（ij）=-e（i）*e（j）*（G（ij）*co（ij）+B（ij）*si（ij））; 
            J1（ij）=e（i）*e（j）*（G（ij）*co（ij）+B（ij）*si（ij））;
            L1（ij）=-e（i）*e（j）*（G（ij）*si（ij）-B（ij）*co（ij））;
        end      
    end
end
H=H1（2:n12:n1）;
N=N1（2:n1r+2:n1）;
J=J1（r+2:n12:n1）;
L=L1（r+2:n1r+2:n1）;
J=[H N;J L];
%雅克比矩阵J
 for i=2:n1
    a（i）=0;b（i）=0;
    for j=1:n1
        a（i）=a（i）+e（j）*（G（ij）*co（ij）+B（ij）*si（ij））;
        if i>=（r+2）
        b（i）=b（i）+e（j）*（G（ij）*si（ij）-B（ij）*co（ij））;
        end
     end
    dp（i-1）=ps（i）-e（i）*a（i）;
    if i>=（r+2）
       dq（i-（r+1））=qs（i）-e（i）*b（i）;
    end
end
F=[dp‘;dq‘];
dx=-J\F;
f（2:n1）=f（2:n1）+dx（1:（m+r））;
e1=e（（r+2）:n1）+dx（（m+r+1）:（2*m+r））.*e（（r+2）:n1）;
e2=e（1:（r+1））;
e=[e2;e1];
end
%disp（‘e=‘）;disp（e‘）;
%disp（‘单位：度 f=‘）;disp（f‘*180/3.1415926）
%求出各节点电压幅值及相角
    a=0;
 for i=1:n1
    a=a+conj（Y（1i））*（e（i）*cos（f（i））-1i*e（i）*sin（f（i）））;
 end
    a=a*（e（1）*cos（f（1））+1i*e（1）*sin（f（1）））;
%disp（‘平衡节点P+jQ=‘）;disp（a）;

%初值计算
 for j=1:n1
     s（j）=0;
     for i=1:n1
     s（j）=s（j）+conj（Y（ji））*（e（i）*cos（f（i））-1i*e（i）*sin（f（i）））;
     end
     s（j）=s（j）*（e（j）*cos（f（j））+1i*e（j）*sin（f（j）））;
 end
 Y1=conj（-s）./（e‘.^2）;
 %disp（‘S=‘）;disp（s）;
 
 Ra=[0 0 0];X_d=[0.0608 0.1198 0.1813];Xq=[0.0969 0.8645 1.2578];
 %Xq=X_d;%不计凸极效应
 v=e.*cos（f）+1i*e.*sin（f）;
 for i=1:3
 I（i）=conj（s（i）/（e（i）*cos（f（i））+1i*e（i）*sin（f（i））））;
 EQ（i）=e（i）*cos（f（i））+1i*e（i）*sin（f（i））+（Ra（i）+1i*Xq（i））*I（i）;
 EQx（i）=real（EQ（i））;EQy（i）=imag（EQ（i））;
 delta（i）=atan（EQy（i）/EQx（i））;
 end
 %disp（‘delta=‘）;disp（delta*180/3.1415926）;
 for i=1:3
     Vx（i）=e（i）*cos（f（i））; Vy（i）=e（i）*sin（f（i））;
     Vdq=[sin（delta（i）） -cos（delta（i））;cos（delta（i）） sin（delta（i））]*[Vx（i）;Vy（i）];
     Vd（i）=Vdq（1）;Vq（i）=Vdq（2）;
     Ix（i）=real（I（i））;Iy（i）=imag（I（i））;
     Idq=[sin（delta（i）） -cos（delta（

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

     文件        814  2014-01-03 16:43  暂态\b1.txt

     文件        611  2014-01-03 16:51  暂态\b2.txt

     文件       6478  2014-12-07 21:39  暂态\zantai.m

     目录          0  2014-12-07 21:39  暂态

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

                 7903                    4