forn=1:
nbus
nadd=1;
ifabs（Zbus（n,n））==0
forI=1:
nbr
ifnadd==1
ifnl（I）==n|nr（I）==n
ifnl（I）==nk=nr（I）;
elseifnr（I）==nk=nl（I）;
end
ifabs（Zbus（k,k））~=0
form=1:
nbus
ifm~=n
Zbus（m,n）=Zbus（m,k）;
Zbus（n,m）=Zbus（m,k）;
else,end
end
Zbus（n,n）=Zbus（k,k）+ZB（I）;tree=tree+1;
nadd=2;ntree（I）=2;
else,end
else,end
else,end
end
else,end
end
end
%增加两个原有总线间的支路阻抗
forn=1:
nbus
forI=1:
nbr
ifntree（I）==1
ifnl（I）==n|nr（I）==n
ifnl（I）==nk=nr（I）;
elseifnr（I）==nk=nl（I）;
end
DM=Zbus（n,n）+Zbus（k,k）+ZB（I）-2*Zbus（n,k）;
forjj=1:
nbus
AP=Zbus（jj,n）-Zbus（jj,k）;
forkk=1:
nbus
AT=Zbus（n,kk）-Zbus（k,kk）;
DELZ（jj,kk）=AP*AT/DM;
end
end
Zbus=Zbus-DELZ;
ntree（I）=2;
else,end
else,end
end
end
（3）symfault.m程序代码：
%此程序用来计算电网的三相对称故障
%计算前需要用户生成节点阻抗矩阵
%节点阻抗矩阵可由函数Zbus=zbuild（zdata）生成
%此程序需要用户按提示输入短路节点编号和过度电阻Zf
%向量V是可选参数，包含节点编号和复数电压
%V可由潮流程序lfgauss,lfnewton,decouple自动生成.
%如果V不存在，程序将默认故障前所有的节点电压标幺值为1.0
%程序可得到总的故障电流,节点电压幅值及输电线路的电流
functionsymfaul（zdata,Zbus,V）
nl=zdata（:
1）;nr=zdata（:
2）;R=zdata（:
3）;
X=zdata（:
4）;
nc=length（zdata（1,:
））;
ifnc>4
BC=zdata（:
5）;
elseifnc==4,BC=zeros（length（zdata（:
1））,1）;
end
ZB=R+j*X;
nbr=length（zdata（:
1））;nbus=max（max（nl）,max（nr））;
ifexist（'V'）==1
iflength（V）==nbus
V0=V;
else,end
else,V0=ones（nbus,1）+j*zeros（nbus,1）;
end
fprintf（'\Three-phasebalancedfaultanalysis\n'）
ff=999;
whileff>0
nf=input（'EnterFaultedBusNo.->'）;
rtn=isempty（nf）;
ifrtn==1;nf=-1;end
whilenf<=0|nf>nbus
fprintf（'FaultedbusNo.mustbebetween1&%g\n',nbus）
nf=input（'EnterFaultedBusNo.->'）;
rtn=isempty（nf）;
ifrtn==1;nf=-1;end
end
rtz=1;
whilertz==1
fprintf（'\nEnterFaultImpedanceZf=R+j*Xin'）
Zf=input（'complexform（forboltedfaultenter0）.Zf='）;
rtz=isempty（Zf）;
end
fprintf（'\n'）
fprintf（'Balancedthree-phasefaultatbusNo.%g\n',nf）
If=V0（nf）/（Zf+Zbus（nf,nf））;
Ifm=abs（If）;Ifmang=angle（If）*180/pi;
fprintf（'Totalfaultcurrent=%8.4fperunit\n\n',Ifm）
%fprintf（'p.u.\n\n',Ifm）
fprintf（'BusVoltagesduringfaultinperunit\n\n'）
fprintf（'BusVoltageAngle\n'）
fprintf（'No.Magnitudedegrees\n'）
forn=1:
nbus
ifn==nf
Vf（nf）=V0（nf）*Zf/（Zf+Zbus（nf,nf））;Vfm=abs（Vf（nf））;angv=angle（Vf（nf））*180/pi;
else,Vf（n）=V0（n）-V0（n）*Zbus（n,nf）/（Zf+Zbus（nf,nf））;
Vfm=abs（Vf（n））;angv=angle（Vf（n））*180/pi;
end
fprintf（'%4g',n）,fprintf（'%13.4f',Vfm）,fprintf（'%13.4f\n',angv）
end
fprintf（'\n'）
fprintf（'LinecurrentsforfaultatbusNo.%g\n\n',nf）
fprintf（'FromToCurrentAngle\n'）
fprintf（'BusBusMagnitudedegrees\n'）
forn=1:
nbus
forI=1:
nbr
ifnl（I）==n|nr（I）==n
ifnl（I）==nk=nr（I）;
elseifnr（I）==nk=nl（I）;
end
ifk==0
Ink=（V0（n）-Vf（n））/ZB（I）;
Inkm=abs（Ink）;th=angle（Ink）;
ifreal（Ink）>0
fprintf（'G'）,fprintf（'%7g',n）,fprintf（'%12.4f',Inkm）
fprintf（'%12.4f\n',th*180/pi）
elseifreal（Ink）==0&imag（Ink）<0
fprintf（'G'）,fprintf（'%7g',n）,fprintf（'%12.4f',Inkm）
fprintf（'%12.4f\n',th*180/pi）
else,end
Ign=Ink;
elseifk~=0
Ink=（Vf（n）-Vf（k））/ZB（I）+BC（I）*Vf（n）;
Inkm=abs（Ink）;th=angle（Ink）;
ifreal（Ink）>0
fprintf（'%7g',n）,fprintf（'%10g',k）,
fprintf（'%12.4f',Inkm）,fprintf（'%12.4f\n',th*180/pi）
elseifreal（Ink）==0&imag（Ink）<0
fprintf（'%7g',n）,fprintf（'%10g',k）,
fprintf（'%12.4f',Inkm）,fprintf（'%12.4f\n',th*180/pi）
else,end
else,end
else,end
end
ifn==nf
fprintf（'%7g',n）,fprintf（'F'）,fprintf（'%12.4f',Ifm）
fprintf（'%12.4f\n',Ifmang）
else,end
end
resp=0;
whilestrcmp（resp,'n'）~=1&strcmp（resp,'N'）~=1&strcmp（resp,'y'）~=1&strcmp（resp,'Y'）~=1
resp=input（'Anotherfaultlocation?
Enter''y''or''n''withinsinglequote->'）;
ifstrcmp（resp,'n'）~=1&strcmp（resp,'N'）~=1&strcmp（resp,'y'）~=1&strcmp（resp,'Y'）~=1
fprintf（'\nIncorrectreply,tryagain\n\n'）,end
end
ifresp=='y'|resp=='Y'
nf=999;
elseff=0;end
end
二、不对称短路故障
不对称故障包括单相接地短路（singleline-to-groundfault），两相短路（line-to-linefault），和两相接地短路（doubleline-to-groundfault）。
这里介绍三个程序来分析不对称短路故障：
单相接地lgfault（zdata0,zbus0,zdata1,zbus1,zdata2,zbus2,V）
两相短路llfault（zdata1,zbus1,zdata2,zbus2,V）
两相接地短路dlgfault（zdata0,zbus0,zdata1,zbus1,zdata2,zbus2,V）
对于lgfault和dlgfault需要输入正序、零序、负序节点阻抗矩阵Zbus1、Zbus0和Zbus2，而llfault只需输入正序和负序节点阻抗矩阵Zbus1、Zbus2，最后参数V是可以选择的。
如果用户未输入V，程序则默认所有故障前的电压为1.0pu。
如果用户输入参数V，则需要输入节点编号以及节点电压的复数值。
变量V可由潮流计算程序获得。
零序网络和正序网络的节点阻抗矩阵分别通过函数Zbus0=zbuild（zdata0）和Zbus1=zbuild（zdata1）获得。
参数zdata1包含正序网络阻抗，参数zdata0包含零序网络阻抗。
参数zdata1、zdata2、zdata0均为e×4阶矩阵，e为元件的个数。
第一列和第二列是节点编号，第三列和第四列分别是线路电阻和电抗（均为标幺值）。
节点0为参考节点，发电机节点到节点0的阻抗为发电机阻抗，可以是发电机的次暂态电抗、暂态电抗或者同步电抗，而且矩阵还包括并联电容和负荷阻抗等。
负序网络和正序网络有同样的拓扑结构。
线路和变压器的负序阻抗等于正序阻抗，而发电机的负序阻抗不等于正序阻抗，但在大型电力系统故障分析中通常按相等来计算。
零序网络的拓扑结构不同于正序网络，零序网络的构建取决于变压器的绕组连接方式，除了两侧中性点都接地的Y-Y连接的变压器，对于其他连接方式的变压器，一次侧和二次侧的零序网络都是隔离开的。
在这样的连接方式的零序网络中,，用inf填写对应的电阻和电抗列。
对于Y-△（中性点接地）连接的变压器，Y侧的阻抗到节点0的数据需要填写。
当中性点接地阻抗为Xn时，需要填写为3Xn。
在不对称故障计算时，还需要一个函数用来求节点阻抗矩阵，Zbus=zbuildpi（linedata，gendata，yload），函数与潮流计算程序是相容的。
第一个参数linedata与潮流程序所需要的数据相同。
第一列和第二列是节点编号；第三列到第五列分别是线路的电阻、电抗、1/2线路电纳（均为标幺值，以指定的MVA为基准功率）；最后一列为变压器分接头位置，对于线路，在这一列中必须填写1。
在潮流计算程序中，发电机电抗不包含在参数linedata中，而是包含在参数gendata中。
参数ge
升级会员 