function buduichen  
%數據來源于教材《電力系統暫態分析》P129例（5-1）
%給學生提供的不對稱短路故障的部分程序。
%要求的工作（1）自動形成導納陣；（2）完善短路類型，兩相相間短路和兩相短路接地；（3）求短路后發電機的電壓
 clear  %清除所有變量的值
N1=3; %輸入網絡的節點數
N2=5;  %輸入包括所有發電機節點的網絡的節點數
k=3;%輸入短路點的節點號
fault=1;  %輸入短路類型 ;f(1)=1單相短路接地;f(2)=2兩相相間短路;f(1,1)=3兩相短路接地
%輸入支路各序阻抗，ZZ1表示支路正序阻抗，ZZ2表示支路負序阻抗，ZZ0表示支路零序阻抗

ZZ1=[j*0.15 j*0.1   j*0.1
     j*0.1  j*0.075 j*0.1
     j*0.1  j*0.1   0    ];
 ZZ0=[j*0.05 j*0.2  j*0.2
      j*0.2  j*0.025  j*0.2
      j*0.2  j*0.2    0   ];
  ZG1=[0      j*0.1   j*0.1   j*0.05    0
      j*0.1    0     j*0.1     0     j*0.025
      j*0.1  j*0.1     0       0       0
      j*0.05   0       0     j*0.1     0
      0      j*0.025   0       0     j*0.05];
  ZZ2=zeros(3,3);
  ZZ2=ZZ1;
  

 %不計及發電機節點的節點導納矩陣
   %求得正序節點導納矩陣
 Y1=zeros(3,3);
 for i=1:N1
     for p=1:N1
         if ZZ1(i,p)~=0
         Y1(i,i)=Y1(i,i)+1/ZZ1(i,p);
            if  i~=p
        Y1(i,p)=-1/ZZ1(i,p);
            end
         end
 end
end
 Y1
   %求得負序節點導納矩陣
 Y2=zeros(3,3);
 Y2=Y1
   %求得零序節點導納矩陣
  Y0=zeros(3,3);
 for i=1:N1
     for p=1:N1
         if ZZ0(i,p)~=0
         Y0(i,i)=Y0(i,i)+1/ZZ0(i,p);  
            if  i~=p
        Y0(i,p)=-1/ZZ0(i,p);
            end
         end
 end
end
Y0
%計及發電機節點的節點導納矩陣
  %YY1和YY2分別是此時的正序、負序節點導納矩陣
 YY1=zeros(5,5);
 for i=1:N2
     for p=1:N2
         if ZG1(i,p)~=0
         YY1(i,i)=YY1(i,i)+1/ZG1(i,p);
            if  i~=p
        YY1(i,p)=-1/ZG1(i,p);
            end
         end
 end
end
 YY1
 YY2=zeros(5,5);
 YY2=YY1


%第一部分：計算所有節點的a,b,c三相電壓

for p=1:N1
    if p==k
        I(p)=1;
    else
        I(p)=0;
    end
end

%正序
Z1=eye(N1)/Y1;% eye(N1)是3階的單位矩陣,求節點導納矩陣Y1的逆陣,即節點阻抗矩陣Z
Z1(:,k)=Z1*I';%短路點k的自阻抗和互阻抗     
Zk1=Z1(:,k);

%負序
Z2=eye(N1)/Y2;
Z2(:,k)=Z2*I';
Zk2=Z2(:,k);

%零序
Z0=eye(N1)/Y0;
Z0(:,k)=Z0*I';     
Zk0=Z0(:,k);


%計根據故障類型選擇不同的計算公式，算故障點各序電流
if fault==1  
    Ik1=1/(Z1(k,k)+Z2(k,k)+Z0(k,k));
    Ik2=Ik1;                               %單相短路接地時的各序電流分量
    Ik0=Ik1;
else
   if fault==2
      Ik1=1/(Z1(k,k)+Z2(k,k));
       Ik2=-1*Ik1;                          %兩相相間短路時的各序電流分量
       Ik0=0;
     else
        if fault==3
         Ik1=1/(Z1(k,k)+Z2(k,k)*Z0(k,k)/(Z2(k,k)+Z0(k,k)));   
         Ik2=-Ik1*Z0(k,k)/(Z2(k,k)+Z0(k,k));        %兩相短路接地時的各序電流分量
         Ik0=-Ik1*Z2(k,k)/(Z2(k,k)+Z0(k,k));
        end
    end
   end
end

%計算各序電流的故障分量
for p=1:N1  
    p;
    if p==k
        I1(p)=-1*Ik1;
        I2(p)=-1*Ik2;
        I0(p)=-1*Ik0;
    else
        I1(p)=0;
        I2(p)=0;
        I0(p)=0;
    end
end

%計算故障點各序電壓的故障分量 P96公式（4-13）
uu1(:,k)=Y1\I1.' %“\”表示左除 A\B=inv(A)*B  inv是求逆陣。 即AX=B
uu2(:,k)=Y2\I2.'
uu0(:,k)=Y0\I0.'

%假設故障前空載，各節點電壓正序分量為1 
for p=1:N1
    U1(p)=1;
end

%計算故障點發生故障后的電壓值:正常運行的分量+故障分量，P37公式（5-39）
u1=U1'+uu1(:,k)  %計算所有節點正序電壓
u2=uu2(:,k)      %計算所有節點負序電壓
u0=uu0(:,k)       %計算所有節點零序電壓

a=-0.5+j*sqrt(3)/2
T=[1    1      1        % T為對稱分量法的變換矩陣，見P94公式（4-4）
   a^2  a      1
   a    a^2    1]
   
for p=1:N1
    U=[u1(p)   u2(p)    u0(p)]
    p
    Uabc=T*U.'       %Uabc表示p節點的a,b,c三相電壓向量,見P94公式（4-4）
    UUabc=abs(Uabc)   %UUabc表示p節點的a,b,c三相電壓有效值
end

%第二部分：計算支路電流
for  m=1:N1
    for n=1:N1
        if m<n
            mn(1)=m;mn(2)=n;
            mn
            I1(m,n)=(u1(m)-u1(n))/ZZ1(m,n);   %支路正序電流
            I2(m,n)=(u2(m)-u2(n))/ZZ2(m,n);   %支路負序電流
            I0(m,n)=(u0(m)-u0(n))/ZZ0(m,n);   %支路零序電流
            Iabc=T*[I1(m,n)   I2(m,n)    I0(m,n)].';%由序分量合成相分量
            Iabc
            abs(Iabc)   
        end       
   end
end

%第三部分：計算發電機的端電壓，只寫了部分程序，余下的工作由同學們完成。要考慮變壓器繞組是11點接線的
for p=1:N2
    if p==k
        II(p)=-Ik1;
    else
        II(p)=0;
    end
end
vv1(:,k)=YY1\II.';
vv2(:,k)=YY2\II.';            
for p=1:N2
    V1(p)=1;
end
%計算發機障后的電壓值=正常運行分量+故障分量
v1=V1'+vv1(:,k)   %計算發電機節點的正序分量
v2=vv2(:,k)       %計算發電機節點的負序分量
%計及11點接線形式
c=sqrt(3)/2+j*0.5  %正序分量三角形側電壓較星型側超前30度（11點鐘）
d=sqrt(3)/2-j*0.5  %負序分量三角形側電壓較星型側落后30度（11點鐘）
v1(4)=v1(4)*c;
v2(4)=v2(4)*d;
v1(5)=v1(5)*c;
v2(5)=v2(5)*d;
VG1=[v1(4)  v2(4)   0 ]
VG1abc=T*VG1.'
VG1abc1=abs(VG1abc)    %VG1abc是發電機1的機端三相電壓
VG2=[v1(5)  v2(5)    0 ]
VG2abc=T*VG2.'
VG2abc1=abs(VG2abc)    %VG2abc是發電機1的機端三相電壓