三相短路电流计算

clear all;
format long e
Y1=[-j*23.619417  0         0            j*4.444444       j*9.523810    j*4.651163
0            -j*35.238095  0             j*28.571429      0            0
0            0            -j*27.809524   0                j*23.80524    0
j*4.444444    j*28.571429   0            -j*41.015873      0             j*8.000000
j*9.523810    0             j*23.809524   0              -j*39.047619   j*5.714286
j*4.651163    0           0             j*8.000000       j*5.714286    -j*18.365449] %输入正序负序导纳矩阵

Y0=[-j*31.3736      0          0          j*2.5        j*5.000000    j*2.564103 
   0             -j*33.8346     0          j*28.5714      0           0
   0              0           -j*9.2083       0       j*5.2083        0
   j*2.5          j*28.5714      0         -j*35.0714     0            j*4  
   j*5            0           j*5.2083     0         -j*5.2083      j*2.857182   
    0           0           0          j*4      j*2.8571      -j*3.2784]        %输入零序导纳矩阵、
disp（‘正序导纳矩阵‘）;
Y1
disp（‘负序导纳矩阵‘）;
Y0
N1=input（‘please input 正负序网络节点数 ：‘）;
N2=input（‘please input 零序网络节点数 ：‘）;
N3=input（‘please input 正负序短路编号（正负序导纳矩阵中短路点对应的编号）  ：‘）;
N4=input（‘please input 零序短路编号（零序导纳矩阵中短路点对应的编号）  ：‘）;
disp（‘fault表示输入短路类型 1表示单相a短路 2表示两相bc短路 3表示两相bc短路接地‘）
fault=input（‘please input fault的值  ：‘）;



%求正序和负序因子表
YY1=zeros（N1N1）;
for m=1:N1
    for n=m :N1 
        if n==m
            YY1（mm）=Y1（mm）;
            for k=1:m-1 
                YY1（mm）=YY1（mm）-YY1（km）*YY1（km）*YY1（kk）;
            end
        else
            YY1（mn）=Y1（mn）;
            for k=1:m-1
                YY1（mn）=YY1（mn）-YY1（km）*YY1（kn）*YY1（kk）; 
            end
            YY1（mn）=YY1（mn）/YY1（mm）;
        end
    end
end
for m=1:N1
    YY1（mm）=1/YY1（mm）;
end
YY2=YY1;
disp（‘正序因子表‘）;
YY1 
disp（‘负序因子表‘）;
YY2

%求零序因子表
YY0=zeros（N2N2）;
for m=1:N2
    for n=m :N2 
        if n==m
            YY0（mm）=Y0（mm）;
            for k=1:m-1 
                YY0（mm）=YY0（mm）-YY0（km）*YY0（km）*YY0（kk）;
            end
        else
            YY0（mn）=Y0（mn）;
            for k=1:m-1
                YY0（mn）=YY0（mn）-YY0（km）*YY0（kn）*YY0（kk）; 
            end
            YY0（mn）=YY0（mn）/YY0（mm）;
        end
    end
end
for m=1:N2
    YY0（mm）=1/YY0（mm）;
end
disp（‘零序因子表‘）;
YY0

%正负序阻抗