NF=input('请输入短路点的数目：NF=');   %输入短路点的数目NF

for i1=1:NF     %给i1赋值短路点的数目，从1到NF

clear; %清除

n1=input('请输入正序节点数：n1=');  %输入正序节点数n1

n2=input('请输入负序节点数：n2=');  %输入负序节点数n2

n0=input('请输入零序节点数：n0=');  %输入零序节点数n0

nl1=input('请输入正序支路数：nl1='); %输入正序支路数nl1

nl2=input('请输入负序支路数：nl2='); %输入负序支路数nl2

nl0=input('请输入零序支路数：nl0='); %输入零序支路数nl0

Lf=input('请输入短路类型：Lf=');  %输入短路类型Lf

f=input('请输入短路节点号：f=');  %输入短路节点号f

If=input('请输入纵横故障标志：If=');  %输入纵横故障标志If

zf=input('请输入接地阻抗：zf=');  %输入接地阻抗zf

zg=input('请输入接地阻抗：zg=');  %输入接地阻抗zg

B1=input('请输入正序支路参数：B1=');  %输入正序支路参数B1

B2=input('请输入负序支路参数：B2=');  %输入负序支路参数B2

B0=input('请输入零序支路参数：B0=');  %输入零序支路参数B0

A1=input('请输入不对称系数：A1=');  %输入不对称系数A1

Vcs=input('请输入由各节点的初始电压标幺值形成的列矩阵：Vcs=');  %输入由各节点的初始电压标幺值形成的列矩阵Vcs

Z1=zeros(n1);Z2=zeros(n2);Z0=zeros(n0); %正序节点阻抗矩阵Z1中正序节点数组n1清零，负序节点阻抗矩阵Z2中负序节点数组n2清零，零序节点阻抗矩阵Z0中零序节点数组n0清零

Y1=zeros(n1);Y2=zeros(n2);Y0=zeros(n0); %正序节点导纳矩阵Z1中正序节点数组n1清零，负序节点导纳矩阵Z2中负序节点数组n2清零，零序节点导纳矩阵Z0中零序节点数组n0清零

V1=zeros(n1,1);V2=zeros(n2,1);V0=zeros(n0,1); %正序节点电压矩阵V1中电压清零，负序节点电压矩阵V2中电压清零，零序节点电压矩阵V0中电压清零

I1=zeros(nl1,1);I2=zeros(nl2,1);I0=zeros(nl0,1); %正序节点电流矩阵I1中电流清零，负序节点电流矩阵I2中电流清零，零序节点电流矩阵I0中电流清零

for m1=1:3 %网络化简为正负零序等值网络

m=0;

if m1==1 %在正序等值网络下参数

nl=nl1;B=B1;Z=Z1;n=n1; %正序支路数nl1、正序支路参数B1、正序阻抗矩阵Z1、正序节点数n1

elseif m1==2 %在负序等值网络下参数

nl=nl2;B=B2;Z=Z2;n=n2; %负序支路数nl2、负序支路参数B2、负序导纳矩阵Z2、负序节点数n2

elseif m1==3 %在零序等值网络下参数

nl=nl0;B=B0;Z=Z0;n=n0; %零序支路数nl0、零序支路参数B0、零序导纳矩阵Z0、零序节点数n0

end

m=0;Z=zeros(n);V=zeros(n);I=zeros(nl); %循环次数m清零，节点阻抗矩阵Z中节点数n清零，节点电压矩阵V中节点数n清零，节点电流矩阵I中支路数nl清零

for k1=1:nl %B矩阵的第一行支路到nl行支路循环

p=B(k1,1);q=B(k1,2); %取支路的第1、2列的值为首末端号

if B(k1,6)==0 %支路首端p未处于高压侧

k=1./B(k1,5);  %支路变比取反

else

k=B(k1,5); %首端p处于高压侧则支路变比不变

end

if p==0 %首端号为0，表示为接地支路

if q>m                   %追加接地树支

Z(q,q)=B(k1,3);m=m+1; %树支支路阻抗为输入阻抗

else                      %追加接地连支

for i=1:m %令i从1到m次循环

Z(i,m+1)=-Z(i,q);Z(m+1,i)=-Z(q,i); %得到非对角元素暂增行与列的阻抗值

end

Z(m+1,m+1)=Z(q,q)+B(k1,3); %得到对角元素上暂增元素的阻抗值

for i=1:m   %令i从1到m次循环

for j=1:m  %令j从1到m+1次循环

Z(i,j)=Z(i,j)-Z(i,m+1)*Z(m+1,j)./Z(m+1,m+1); %第51行运用公式：Z'ij=Zij-(Zi(m+1)*Z(m+1)j/Z(m+1)(m+1))用高斯消去法消去矩阵暂增行与列后得到原矩阵的元素阻抗值

end

Z(i,m+1)=0; %暂增列数值归零

end

for i=1:m+1 %令i从1到m+1次循环

Z(m+1,i)=0; %暂增行数值归零

end

end

else if q>m                  %追加不接地树支

for i=1:m  %令i从1到m次循环

Z(i,q)=Z(i,p)*k;Z(q,i)=Z(p,i)*k; %得到非对角元素行与列的阻抗值

end

Z(q,q)=k^2*Z(p,p)+k^2*B(k1,3);%得到对角元素上的阻抗值

m=m+1;

else

for i=1:m   %令i从1到m次循环       %追加不接地连支

Z(i,m+1)=k*Z(i,p)-Z(i,q); %得到非对角元素暂增列的阻抗值

Z(m+1,i)=k*Z(p,i)-Z(q,i); %得到非对角元素暂增行的阻抗值

end

Z(m+1,m+1)=k^2*Z(p,p)+Z(q,q)-2*k*Z(p,q)+k^2*B(k1,3); %得到对角元素上暂增元素的阻抗值

for i=1:m %令i从1到m次循环

for j=1:m  %令j从1到m次循环

Z(i,j)=Z(i,j)-Z(i,m+1)*Z(m+1,j)./Z(m+1,m+1);%用高斯消去法消去矩阵暂增行与列后得到原矩阵的元素阻抗值

end

Z(i,m+1)=0; %暂增列数值归零

end

for i=1:m+1 %令i从1到m+1次循环

Z(m+1,i)=0; %暂增行数值归零

end

end

end

end

if m1==1  %正序等值网络下

Z1=Z;  %形成正序阻抗矩阵

elseif m1==2 %负序等值网络下

Z2=Z; %形成负序阻抗矩阵

elseif m1==3 %零序等值网络下

Z0=Z; %形成零序阻抗矩阵

end

end %以上程序求出正，负，零序对应的阻抗阵

if If==0 %求出横向故障的阻抗阵

Z1(f,f)=Z1(f,f); %第92行运用公式：Zq(F,F)=Zq(f,f) (q=1,2,0)

Z2(f,f)=Z2(f,f);

Z0(f,f)=Z0(f,f);

else  %求出纵向故障的阻抗阵

Z1(f,f)=Z1(f,f)+Z1(1,1)-2*Z1(f,1); %第96行运用公式：Zq(F,F)=Zq(f,f)+Zq(f',f')-2*Zq(f,f') (q=1,2,0)

Z2(f,f)=Z2(f,f)+Z2(1,1)-2*Z2(f,1);

Z0(f,f)=Z0(f,f)+Z0(1,1)-2*Z0(f,1);

end

if Lf==1 %单相短路时

Z4=Z2(f,f)+Z0(f,f)+3*zf;k2=1;k0=1;I1(f)=Vcs(f)./(Z1(f,f)+Z4);%第101行运用公式：附加阻抗Z4=Z2(F,F)+Z0(F,F)+3*Zf,系数k2=1,k0=1,I1(F)=UF(0)/(Z1(F,F)+Z4)

I2(f)=k2*I1(f);I0(f)=k0*I1(f);a1=I1(f);a2=I2(f);a0=I0(f); %第101行运用公式：I2(F)=k2*I1(F),I0(F)=k0*I1(F),将求得的故障口的正负零序电流的值赋给a1、a2、a0

elseif Lf==2 %两相短路时

Z4=Z2(f,f)+2*zf;k2=-1;k0=0;I1(f)=Vcs(f)./(Z1(f,f)+Z4); %第104行运用公式：附加阻抗Z4=Z2(F,F)+2*Zf，系数k2=-1,k0=0,I1(F)=UF(0)/(Z1(F,F)+Z4)

I2(f)=k2*I1(f);I0(f)=k0*I1(f);a1=I1(f);a2=I2(f);a0=I0(f);%第105行运用公式：I2(F)=k2*I1(F),I0(F)=k0*I1(F),将求得的故障口的正负零序电流的值赋给a1、a2、a0

elseif Lf==3 %两相短路接地时

Z4=zf+(Z2(f,f)+zf)*(Z0(f,f)+zf+3*zg)./(Z2(f,f)+Z0(f,f)+2*zf+3*zg);  %第107行运用公式：附加阻抗Z4=Zf+(Z2(F,F)+Zf)*(Z0(F,F)+Zf+3*Zg)/(Z2(F,F)+Z0(F,F)+2*Zf+3*Zg),

k2=-(Z0(f,f)+zf+3*zg)./(Z2(f,f)+Z0(f,f)+2*zf+3*zg); %第108行运用公式：系数k2=(Z0(F,F)+Zf+3*Zg)/(Z2(F,F)+Z0(F,F)+2*Zf+3*Zg)

k0=-(Z2(f,f)+zf)./(Z2(f,f)+Z0(f,f)+2*zf+3*zg); %第109行运用公式：系数k0=(Z2(F,F)+Zf)/(Z2(F,F)+Z0(F,F)+2*Zf+3*Zg)

I1(f)=Vcs(f)./(Z1(f,f)+Z4); %第110行运用公式：I1(F)=UF(0)/(Z1(F,F)+Z4)

I2(f)=k2*I1(f);I0(f)=k0*I1(f);a1=I1(f);a2=I2(f);a0=I0(f);%第111行运用公式:I2(F)=k2*I1(F),I0(F)=k0*I1(F),将求得的故障口的正负零序电流的值赋给a1、a2、a0

elseif Lf==4 %单相断线时

Z4=zf+(Z2(f,f)+zf)*(Z0(f,f)+zf)./(Z2(f,f)+Z0(f,f)+2*zf); %第113行运用公式：附加阻抗Z4=Zf+(Z2(F,F)+Zf)*(Z0(F,F)+Zf)/(Z2(F,F)+Z0(F,F)+2*Zf)

k2=-(Z0(f,f)+zf)./(Z2(f,f)+Z0(f,f)+2*zf); %第114行运用公式：系数k2=-(Z0(F,F)+Zf)/(Z2(F,F)+Z0(F,F)+2*Zf)

k0=-(Z2(f,f)+zf)./(Z2(f,f)+Z0(f,f)+2*zf); %第115行运用公式：系数k0=-((Z2(F,F)+Zf)/(Z2(F,F)+Z0(F,F)+2*Zf)

I1(f)=Vcs(f)./(Z1(f,f)+Z4);%第116行运用公式：I1(F)=UF(0)/(Z1(F,F)+Z4)

I2(f)=k2*I1(f);I0(f)=k0*I1(f);a1=I1(f);a2=I2(f);a0=I0(f); %第117行运用公式：I2(F)=k2*I1(F),I0(F)=k0*I1(F),将求得的故障口的正负零序电流的值赋给a1、a2、a0

elseif Lf==5 %两相断线时

Z4=Z2(f,f)+Z0(f,f)+3*zf;k2=1;k0=1;I1(f)=Vcs(f)./(Z1(f,f)+Z4); %第119行运用公式：附加阻抗Z4=Z2(F,F)+Z0(F,F)+3*Zf，系数k2=1,k0=1,I1(F)=UF(0)/(Z1(F,F)+Z4)

I2(f)=k2*I1(f);I0(f)=k0*I1(f);a1=I1(f);a2=I2(f);a0=I0(f); %第120行运用公式：I2(F)=k2*I1(F),I0(F)=k0*I1(F),将求得的故障口的正负零序电流的值赋给a1、a2、a0

end

disp('故障点处的正序电流标幺值为：');

disp(a1); %输出故障点处的正序电流标幺值

disp('故障点处的负序电流标幺值为：');

disp(a2); %输出故障点处的负序电流标幺值

disp('故障点处的零序电流标幺值为：');

disp(a0); %输出故障点处的零序电流标幺值

for i=1:n1 %各节点的正序电压标幺值

V1(i)=Vcs(i)-I1(f)*Z1(i,f); %第129行运用公式：Ui=Ui(0)-ZifIf，i=1,2,3,…,n

end

for i=1:n2 %各节点的负序电压标幺值

V2(i)=-I2(f)*Z2(i,f); %第132行运用公式：Ui=Ui(0)-ZifIf，Ui(0)=0,i=1,2,3,…,n

end

for i=1:n0 %各节点的零序电压标幺值

V0(i)=-I0(f)*Z0(i,f); %第135行运用公式：Ui=Ui(0)-ZifIf，Ui(0)=0,i=1,2,3,…,n

end

disp('各节点的正序电压标幺值(节点号从小到大)依次为：');

disp(V1); %输出各节点的正序电压标幺值

disp('各节点的负序电压标幺值(节点号从小到大)依次为：');

disp(V2); %输出各节点的负序电压标幺值

disp('各节点的零序电压标幺值(节点号从小到大)依次为：');

disp(V0);   %输出各节点的零序电压标幺值

for i=1:3   %求故障点A，B，C三相各相的总电流，各个节点A，B，C三相各相的总电压，各支路A，B，C三相各相的总电流

if i==1  %在正序等值网络下

B=B1;V=V1;nl=nl1; %列出节点正序导纳矩阵，节点的正序电压标幺值矩阵，正序支路数

elseif i==2  %在负序等值网络下

B=B2;V=V2;nl=nl2; %列出节点负序导纳矩阵，节点的负序电压标幺值矩阵，负序支路数

elseif i==3  %在零序等值网络下

B=B0;V=V0;nl=nl0; %列出节点零序导纳矩阵，节点的零序电压标幺值矩阵，零序支路数

end

for j=1:nl %求出各支路电流的标幺值

if B(j,6)==0 %确定变比

k=B(j,5); %支路变比不变

else

k=1./B(j,5); %支路变比取反

end

p=B(j,1);q=B(j,2); %取取支路的第1、2列的值为首末端号

if p==0  %求接地支路电流标幺值

e=0;b=B(j,3); %支路首端号为0的节点处的电压为0；

I(j)=(e-V(q)./k)./b; %第160行运用公式：Ipq=((0-Uq)/k)/Zpq

else

I(j)=(V(p)-V(q)./k)./B(j,3); %第160行运用公式: Ipq=(Up-Uq)/Zpq

end

end

if i==1  %在正序等值网络下时

I1=I;  %列出各支路的正序电流标幺值

elseif i==2  %在负序等值网络下时

I2=I;   %列出各支路的负序电流标幺值

elseif i==3  %在零序等值网络下时

I0=I;  %列出各支路的零序电流标幺值

end

end

disp('各支路的正序电流标幺值(顺序同输入B时一样)依次为：');

for i=1:nl %从B1矩阵的第一行支路到nl行支路循环

disp(I1(i)); %输出各支路的正序电流标幺值

end

disp('各支路的负序电流标幺值(顺序同输入B时一样)依次为：');

for i=1:nl %从B2矩阵的第一行支路到nl行支路循环

disp(I2(i)); %输出各支路的负序电流标幺值

end

disp('各支路的零序电流标幺值(顺序同输入B时一样)依次为：');

for i=1:nl  %从B0矩阵的第一行支路到nl行支路循环

disp(I0(i)); %输出各支路的零序电流标幺值

end

for i=1:3   %给i赋值1到3

if i==1 %汇总故障口的正、负、零序电流

s0=a0;s1=a1;s2=a2; %列出已求得故障点的正、负、零序电流

elseif i==2 %汇总各节点的正，负，零序电压标幺值

s0=V0;s1=V1;s2=V2; %列出已求得各节点的正、负、零序电压标幺值

else

s0=I0;s1=I1;s2=I2;  %列出已求得各支路的正、负、零序电流标幺值

end

for j=1:3   %给j赋值1到3

CA=s0*A1(j,1)+s1*A1(j,2)+s2*A1(j,3); %计算故障点A，B，C三相各相的总电流，各个节点A，B，C三相各相的总电压，各支路A，B，C三相各相的总电流

if j==1   %如果j==1

A=CA; %求得故障点A相的总电流，各个节点A相的总电压，各支路A相的总电流

elseif j==2   %如果j==2

B=CA;  %求得故障点B相的总电流，各个节点B相的总电压，各支路B相的总电流

elseif j==3   %如果j==3

C=CA;   %求得故障点C相的总电流，各个节点C相的总电压，各支路C相的总电流

end

end

if i==1  %汇总故障点处三相各相的总电流

Iadl=A;Ibdl=B;Icdl=C; %故障点处A相的总电流，故障点处B相的总电流，故障点处C相的总电流

elseif i==2   %汇总各个节点A，B，C三相各相的总电压

Vajy=A;Vbjy=B;Vcjy=C;  %各节点A相的总电压，各节点B相的总电压，各节点C相的总电压

elseif i==3   %汇总各支路A，B，C三相各相的总电流

Iazl=A;Ibzl=B;Iczl=C;  %各支路A相的总电流，各支路B相的总电流，各支路C相的总电流，

end

end

disp('故障点处A相电流Iadl为：');

disp(Iadl); %输出故障点处A相电流

disp('故障点处B相电流Iadl为：');

disp(Ibdl); %输出故障点处B相电流

disp('故障点处C相电流Iadl为：');

disp(Icdl); %输出故障点处C相电流

disp('各节点处A相电压Vajy(节点号从小到大)依次为：');

disp(Vajy);%输出各节点处A相电压

disp('各节点处B相电压Vbjy(节点号从小到大)依次为：');

disp(Vbjy); %输出各节点处B相电压

disp('各节点处C相电压Vcjy(节点号从小到大)依次为：');

disp(Vcjy); %输出各节点处C相电压

disp('各支路A相电流Iazl(顺序同输入B时一样)依次为：');

for i=1:nl  %从B1、B2、B0矩阵的第一行支路到nl行支路循环

disp(Iazl(i)); %输出各支路A相电流Iazl

end

disp('各支路B相电流Ibzl(顺序同输入B时一样)依次为：');

for i=1:nl   %从从B1、B2、B0矩阵的第一行支路到nl行支路循环

disp(Ibzl(i)); %输出各支路B相电流

end

disp('各支路C相电流Iczl(顺序同输入B时一样)依次为：');

for i=1:nl   %从从B1、B2、B0矩阵的第一行支路到nl行支路循环

disp(Iczl(i)); %输出各支路C相电流

end

end