电力系统短路故障分析的MATLAB辅助程序设计短路计算程序

1、电力系统短路故障分析的MATLABffi助程序设计电力系统短路故障可分为三相对称短路故障(three-phase balanced faults)和不对称短路故障(unbalanced faults )。不对称短路故障又分为单相接地短路故障 (single line-to-ground fault) 、两相短路故障(line-to-line fault)以及两相接地短路故障(double line-to-ground fault) 。根据故障分析结果可以对继电保护装置、自动装置进行整定计算，我们可以建立算法来形成节点阻抗矩阵 , 利用节点阻抗矩阵来计算短路故障情况下的节点电压和线路电流。一、三

2、相对称短路故障进行三相短路计算需要两个程序：zbuild / zbuildpi 和 symfault 程序 ， zbuild 、zbuildpi用来在MATLA冲形成节点阻抗矩阵，symfault用来计算三相对称故障。Zbus=zbuild(zdata) 这里的参数zdata是一个(eX4)阶矩阵，e是拓扑图的总支路数目。第一列和第二列为元素两端的节点编号，第三列和第四列分别是线路的电阻、电抗的标幺值。连接在 0 节点和发电机节点之间的发电机阻抗可能是次暂态电抗、暂态电抗或同步电抗，而且这个矩阵中还包含并联电抗器和负荷阻抗。Zbus=zbuildpi(linedata,gendata,yloa

3、d ) 这个函数与潮流计算程序是相容的，第一个参数linedat a 与潮流计算程序中的文件是一致的。第一列和第二列为节点编号；第三列到第五列分别是线路的电阻、电抗以及1/2 线路电纳值，这三项都为在统一基准容量下的标幺值；最后一列是变压器分接头位置，对线路来说，必须输入1；线路无输入顺序。发电机参数不包含在Linedata 参数中，而是包含在第二个参数gendata中，gendata是一个gX4阶矩阵,g是发电机总数。 第一列和第二列为0 节点、 发电机节点编号，第三列和第四列为发电机的暂态电阻和暂态电抗。最后一个参数yload 是可选择的，这个矩阵有两列，第一列为节点编号，第二列为复数导纳

4、值，yload 可以由潮流程序lfguss ， lfnewton或者 decouple 自动生成。Zbuild 和 zbuildpi 两个函数可以通过建立算法求出节点阻抗矩阵。首先添加所有与参考节点相连的树支，然后添加其余的树支，最后添加共轭连支。程序 symfault(zdata,Zbus,V )用来计算三相对称故障，程序要求输入zdata和Zbus两个矩阵，第三个参数 V是可选的。如果V不存在，程序将默认故障前 所有的节点电压标幺值为1.0,如果变量V存在，那么V包括节点编号和复数电 压值。电压向量V也可以由潮流计算程序自动生成。当 symfault程序运行时， 用户要输入故障节点编号和故

5、障阻抗，运行可得到总的故障电流，节点电压幅值以及故障情况下的线路电流。在三相短路计算中，zbuild 和 symfault 程序， zbuildpi 和 symfault 程序 都可以进行计算，下面是三相短路计算使用的程序代码：( 1) Zbuild.m 程序代码：function Zbus = zbuild(linedata)nl = linedata(:,1); nr = linedata(:,2); R = linedata(:,3);X = linedata(:,4);nbr=length(linedata(:,1); nbus = max(max(nl), max(nr);for k

6、=1:nbrif R(k) = inf | X(k) =infR(k) = 999999999; X(k) = 999999999;%无穷else , endendZB = R + j*X;Zbus = zeros(nbus, nbus);tree=0;%从参考总线0 上添加一个分支for I = 1:nbrntree(I) = 1;if nl(I) = 0 | nr(I) = 0if nl(I) = 0 n = nr(I);elseif nr(I) = 0 n = nl(I);endif abs(Zbus(n, n) = 0 Zbus(n,n) = ZB(I);tree=tree+1;%ne

7、welse Zbus(n,n) = Zbus(n,n)*ZB(I)/(Zbus(n,n) + ZB(I);endntree(I) = 2;else , endend% 添加一个新总线分支到现有总线上while tree < nbusfor n = 1:nbusnadd = 1;if abs(Zbus(n,n) = 0for I = 1:nbrif nadd = 1;if nl(I) = n | nr(I) = nif nl(I) = n k = nr(I);elseif nr(I) = n k = nl(I);endif abs(Zbus(k,k) = 0for m = 1:nbusif

8、 m = nZbus(m,n) = Zbus(m,k);Zbus(n,m) = Zbus(m,k);else , endendZbus(n,n) = Zbus(k,k) + ZB(I); tree=tree+1;nadd = 2; ntree(I) = 2;else , end28 / 28elseendelse , endendelse , endendend%增加两个原有总线间的支路阻抗for n = 1:nbusfor I = 1:nbrif ntree(I) = 1if nl(I) = n | nr(I) = nif nl(I) = n k = nr(I);elseif nr(I) =

9、 n k = nl(I);endDM = Zbus(n,n) + Zbus(k,k) + ZB(I) - 2*Zbus(n,k);for jj = 1:nbusAP = Zbus(jj,n) - Zbus(jj,k);for kk = 1:nbusAT = Zbus(n,kk) - Zbus(k, kk);DELZ(jj,kk) = AP*AT/DM;endendZbus = Zbus - DELZ;ntree(I) = 2;else , endelse , endendend(2) Zbuildpi.m 程序代码：% This program forms the complex bus im

10、pedance matrix by the method % of building algorithm. Bus zero is taken as reference.% This program is compatible with power flow data.function Zbus, linedata = zbuildpi(linedata, gendata, yload)ng=length(gendata(:,1);nlg=gendata(:,1);nrg=gendata(:,2);zg= gendata(:,2) + j*gendata(:,3);nl = linedata(

11、:,1); nr = linedata(:,2); R = linedata(:,3);X = linedata(:,4);nbr=length(linedata(:,1); nbus = max(max(nl), max(nr);nc = length(linedata(1,:);for k=1:nbrif R(k) = inf | X(k) = infR(k) = 99999999; X(k) = 99999999;else , endendif nc > 4BC = linedata(:,5);for n = 1:nbusyc(n) = 0;nlc(n) = 0; nrc(n) =

12、 n;for k = 1:nbrif nl(k) = n | nr(k) = n yc(n) = yc(n) + j*BC(k);else , endendendelseif nc=4 yc= zeros(1, nbr);endnlc=nlc' nrc=nrc' yc=yc.'ZB = R + j*X;if exist( 'yload' ) = 1yload = yload.'yc =yc + yload;else , endm = 0;for n = 1:nbusif abs(yc(n) =0m=m+1;nlcc(m) = nlc(n);nrc

13、c(m) = nrc(n);zc(m) = 1/yc(n);else , endendnlcc=nlcc' nrcc=nrcc' zc=zc.'nl=nlg; nlcc; nl; nr = nrg; nrcc; nr; ZB = zg; zc; ZB;linedata=nl nr real(ZB) imag(ZB);nbr= length(nl);Zbus = zeros(nbus, nbus);tree=0;%从参考总线0上添加一个分支for I = 1:nbrntree(I) = 1;if nl(I) = 0 | nr(I) = 0if nl(I) = 0 n =

14、nr(I);elseif nr(I) = 0 n = nl(I);endif abs(Zbus(n, n) = 0 Zbus(n,n) = ZB(I); tree=tree+1;else Zbus(n,n) = Zbus(n,n)*ZB(I)/(Zbus(n,n) + ZB(I); endntree(I) = 2;else , endend%添加一个新总线分支到现有总线上while tree < nbusfor n = 1:nbusnadd = 1;if abs(Zbus(n,n) = 0for I = 1:nbrif nadd = 1if nl(I) = n | nr(I) = nif

15、 nl(I) = n k = nr(I);elseif nr(I) = n k = nl(I);endif abs(Zbus(k,k) = 0for m = 1:nbusif m = nZbus(m,n) = Zbus(m,k);Zbus(n,m) = Zbus(m,k);else , endendZbus(n,n) = Zbus(k,k) + ZB(I); tree=tree+1;nadd = 2; ntree(I) = 2;else , endelse , endelse , endendelse , endendend%增加两个原有总线间的支路阻抗for n = 1:nbusfor I

16、= 1:nbrif ntree(I) = 1if nl(I) = n | nr(I) = nif nl(I) = n k = nr(I);elseif nr(I) = n k = nl(I); endDM = Zbus(n,n) + Zbus(k,k) + ZB(I) - 2*Zbus(n,k);for jj = 1:nbusAP = Zbus(jj,n) - Zbus(jj,k);for kk = 1:nbusAT = Zbus(n,kk) - Zbus(k, kk);DELZ(jj,kk) = AP*AT/DM;endendZbus = Zbus - DELZ;ntree(I) = 2;e

17、lse , endelse , end endend(3)symfault.m 程序代码：% 此程序用来计算电网的三相对称故障% 计算前需要用户生成节点阻抗矩阵% 节点阻抗矩阵可由函数Zbus = zbuild(zdata)生成%此程序需要用户按提示输入短路节点编号和过度电阻Zf%向量V是可选参数，包含节点编号和复数电压%VW由潮流程序 Ifgauss, lfnewton,decouple自动生成.%如果V不存在，程序将默认故障前所有的节点电压标幺值为1.0% 程序可得到总的故障电流, 节点电压幅值及输电线路的电流function symfaul(zdata, Zbus, V)nl = zda

18、ta(:,1); nr = zdata(:,2); R = zdata(:,3);X = zdata(:,4);nc = length(zdata(1,:);if nc > 4BC = zdata(:,5);elseif nc =4, BC = zeros(length(zdata(:,1), 1);endZB = R + j*X;nbr=length(zdata(:,1); nbus = max(max(nl), max(nr);if exist( 'V' ) = 1if length(V) = nbusV0 = V;else , endelse , V0 = ones

19、(nbus, 1) + j*zeros(nbus, 1); endfprintf( 'Three-phase balanced fault analysis n')ff = 999; while ff > 0nf = input( 'Enter Faulted Bus No. -> ');rtn=isempty(nf);if rtn=1; nf=-1;endwhile nf <= 0 | nf > nbusfprintf( 'Faulted bus No. must be between 1 & %g n', nb

20、us)nf = input('Enter Faulted Bus No. -> ');rtn=isempty(nf); if rtn=1; nf=-1;endendrtz=1; while rtz=1fprintf( 'nEnter Fault Impedance Zf = R + j*X in ');Zf = input('complex form (for bolted fault enter 0). Zf = 'rtz=isempty(Zf);endfprintf( ' n')fprintf( 'Balance

21、d three-phase fault at bus No. %gn', nf)If = V0(nf)/(Zf + Zbus(nf, nf);Ifm = abs(If); Ifmang=angle(If)*180/pi;fprintf( 'Total fault current = %8.4f per unit nn', Ifm)%fprintf(' p.u. nn', Ifm)fprintf( 'Bus Voltages during fault in per unit nn')fprintf( ' Bus Voltage An

22、glen')fprintf( ' No. Magnitude degreesn')for n = 1:nbus if n=nfVf(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;endfprintf( ' %4g', n), fprintf('%

23、13.4f', Vfm),fprintf('%13.4fn'angv) endfprintf(n')fprintf( 'Line currents for fault at bus No. %gnn, nf)fprintf(From To Current Anglen'fprintf( ' Bus Bus Magnitude degreesn for n= 1:nbusfor I = 1:nbrif nl(I) = n | nr(I) = nif nl(I) =n k = nr(I);elseif nr(I) = n k = nl(I);

24、endif k=0Ink = (V0(n) - Vf(n)/ZB(I);Inkm = abs(Ink); th=angle(Ink);if real(Ink) > 0'%7g',n),fprintf('%12.4f''%7g',n),fprintf('%12.4f'fprintf(' G '), fprintf(Inkm)fprintf('%12.4fn', th*180/pi)elseif real(Ink) =0 & imag(Ink) < 0fprintf(' G

25、'), fprintf(Inkm)fprintf('%12.4fn', th*180/pi)else , endIgn=Ink;elseif k = 0Ink = (Vf(n) - Vf(k)/ZB(I)+BC(I)*Vf(n);Inkm = abs(Ink); th=angle(Ink);if real(Ink) > 0fprintf('%7g' , n), fprintf('%10g', k),fprintf('%12.4f', Inkm),fprintf('%12.4fn', th*180/pi

26、)elseif real(Ink) =0 & imag(Ink) <0fprintf('%7g' , n), fprintf('%10g', k),fprintf('%12.4f', Inkm),fprintf('%12.4fn', th*180/pi)else , endelse, endelse , endendif n=nfF' ), fprintf('%12.4f', Ifm)fprintf( '%7g' ,n), fprintf(fprintf( '%12.4

27、fn', Ifmang)else , endendresp=0;while strcmp(resp,n' )=1 & strcmp(resp,N' )=1 & strcmp(resp,'y' )=1 & strcmp(resp,'Y' )=1resp = input( 'Another fault location? quote -> ');if strcmp(resp, 'n' )=1 & strcmp(resp, & strcmp(resp, 'Y&

28、#39; )=1Enter ''y'' or ''n'' within single'N' )=1 & strcmp(resp, 'y' )=1), endfprintf( 'n Incorrect reply, try again nnendif resp ='y'| resp ='Y'nf = 999;else ff = 0; endend二、不对称短路故障不对称故障包括单相接地短路(single line-to-ground fault) ，两相

29、短路(line-to-line fault ) ，和两相接地短路(double line-to-ground fault) 。这里介绍三个程序来分析不对称短路故障：单相接地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， 而 llfaul

30、t 只需输入正序和负序节点阻抗矩阵Zbus1、 Zbus2，最后参数V是可以选择的。如果用户未输入 V,程序则默认所有故障前的电压为 1.0puo如果用户输入参数 V,则需要输入节点编号以及节点电压的复数值。变 量 V 可由潮流计算程序获得。零序网络和正序网络的节点阻抗矩阵分别通过函数Zbus0=zbuild ( zdata0 )和 Zbus1=zbuild ( zdata1 ) 获得。 参数 zdata1 包含正序网络阻抗，参数 zdata0包含零序网络阻抗。参数 zdatal、zdata2、zdataO均为eX4阶矩阵，e为元件 的个数。 第一列和第二列是节点编号，第三列和第四列分别是线路

31、电阻和电抗(均为标幺值)。节点 0 为参考节点，发电机节点到节点0 的阻抗为发电机阻抗，可以是发电机的次暂态电抗、暂态电抗或者同步电抗，而且矩阵还包括并联电容和负荷阻抗等。负序网络和正序网络有同样的拓扑结构。线路和变压器的负序阻抗等于正序阻抗， 而发电机的负序阻抗不等于正序阻抗，但在大型电力系统故障分析中通常按相等来计算。零序网络的拓扑结构不同于正序网络，零序网络的构建取决于变压器的绕组连接方式，除了两侧中性点都接地的 Y-Y连接的变压器，对于其他连接方式的变压器，一次侧和二次侧的零序网络都是隔离开的。在这样的连接方式的零序网络中用inf填写对应的电阻和电抗列。对于 Y-A (中性点接地)连接

32、的变压器，Y 侧的阻抗到节点0 的数据需要填写。当中性点接地阻抗为Xn 时，需要填写为3Xn。在不对称故障计算时，还需要一个函数用来求节点阻抗矩阵，Zbus=zbuildpi( linedata ， gendata ， yload ) ， 函数与潮流计算程序是相容的。第一个参数linedata 与潮流程序所需要的数据相同。第一列和第二列是节点编号；第三列到第五列分别是线路的电阻、电抗、1/2 线路电纳(均为标幺值，以指定的MVA为基准功率)；最后一列为变压器分接头位置，对于线路，在这一列中必须填写1。 在潮流计算程序中，发电机电抗不包含在参数linedata 中， 而是包含在参数gendata

33、中。参数gendata是一个(g x4)阶矩阵，其四列分别为参考节点0、发电机节点编号、发电机电阻和电抗。最后一个参数yload 是可选的，这是一个两列的矩阵，包括节点编号和负荷复数导纳值。这些数据可以从潮流计算程序lfgauss ， lfnewton ， decouple 执行过程中自动生成。程序运行时会提示用户输入故障节点编号和故障阻抗Zf,运行结果可得到故障电流，节点电压和线路电流。不对称短路故障计算使用的程序如下：( 1) lgfault.m 程序代码：% The program lgfault is designed for the single line-to-ground% fa

34、ult analysis of a power system network. The program requires% the positive-, negative- and zero-sequence bus impedance matrices,% Zbus1 Zbus2,and Zbus0.The bus impedances matrices may be defined% by the user, obtained by the inversion of Ybus or it may be% determined either from the function Zbus =

35、zbuild(zdata)% or the function Zbus = zbuildpi(linedata, gendata, yload).% The program prompts the user to enter the faulted bus number% and the fault impedance Zf. The prefault bus voltages are% defined by the reserved Vector V. The array V may be defined or% it is returned from the power flow prog

36、rams lfgauss, lfnewton,% decouple or perturb. If V does not exist the prefault bus voltages% are automatically set to 1.0 per unit. The program obtains the% total fault current, bus voltages and line currents during the fault.function lgfault(zdata0, Zbus0, zdata1, Zbus1, zdata2, Zbus2, V) if exist(

37、 'zdata2' ) = 1zdata2=zdata1;else , endif exist( 'Zbus2' ) = 1Zbus2=Zbus1;else , endnl = zdata1(:,1); nr = zdata1(:,2);nl0 = zdata0(:,1); nr0 = zdata0(:,2);nbr=length(zdata1(:,1); nbus = max(max(nl), max(nr);nbr0=length(zdata0(:,1);R0 = zdata0(:,3); X0 = zdata0(:,4);R1 = zdata1(:,3);

38、 X1 = zdata1(:,4);R2 = zdata1(:,3); X2 = zdata1(:,4);for k=1:nbr0if R0(k)=inf | X0(k) =infR0(k) = 99999999; X0(k) = 99999999; else , endendZB1 = R1 + j*X1; ZB0 = R0 + j*X0;ZB2 = R2 + j*X2;if exist( 'V' ) = 1 if length(V) = nbusV0 = V; else , endelse , V0 = ones(nbus, 1) + j*zeros(nbus, 1); e

39、ndfprintf( 'nLine-to-ground fault analysis n')ff = 999; while ff > 0nf = input( 'Enter Faulted Bus No. -> ');rtn=isempty(nf);if rtn=1; nf=-1;endwhile nf <= 0 | nf > nbusfprintf( 'Faulted bus No. must be between 1 & %g n', nbus)nf = input('Enter Faulted Bus

40、 No. -> ');rtn=isempty(nf);if rtn=1; nf=-1;endend rtz=1;while rtz=1fprintf( 'nEnter Fault Impedance Zf = R + j*X in ')Zf = input('complex form (for bolted fault enter 0). Zf = ');rtz=isempty(Zf);endfprintf( ' n')fprintf( 'Single line to-ground fault at bus No. %gn&

41、#39;, nf)a =cos(2*pi/3)+j*sin(2*pi/3);sctm = 1 1 1; 1 aA2 a; 1 a aA2;Ia0 = V0(nf)/(Zbus1(nf,nf)+Zbus2(nf,nf)+ Zbus0(nf, nf)+3*Zf); Ia1=Ia0;Ia2=Ia0;I012=Ia0; Ia1; Ia2;Ifabc = sctm*I012;Ifabcm = abs(Ifabc);fprintf( 'Total fault current = %9.4f per unitnn', Ifabcm(1)fprintf( 'Bus Voltages d

42、uring the fault in per unit nn')fprintf( ' Bus Voltage Magnituden')fprintf(No. Phase a Phase b Phase c n for n = 1:nbusVf0(n)= 0 - Zbus0(n, nf)*Ia0;Vf1(n)= V0(n) - Zbus1(n, nf)*Ia1;Vf2(n)= 0 - Zbus2(n, nf)*Ia2;Vabc = sctm*Vf0(n); Vf1(n); Vf2(n);Va(n)=Vabc(1); Vb(n)=Vabc(2); Vc(n)=Vabc(3)

43、; fprintf(' %5g',n)fprintf( ' %11.4f', abs(Va(n),fprintf(' %11.4f', abs(Vb(n)fprintf(' %11.4fn', abs(Vc(n)end fprintf( ' n')fprintf( 'Line currents for fault at bus No. %gnn', nf)fprintf( ' From To Line Current Magnituden')fprintf( ' Bus Bu

44、s Phase a Phase b Phase c n')for n= 1:nbusfor I = 1:nbrif nl(I) = n | nr(I) = n if nl(I) =n k = nr(I);elseif nr(I) = n k = nl(I); endif k = 0Ink1(n, k) = (Vf1(n) - Vf1(k)/ZB1(I);Ink2(n, k) = (Vf2(n) - Vf2(k)/ZB2(I); else , end else , end end for I = 1:nbr0if nl0(I) = n | nr0(I) = n if nl0(I) =n

45、k = nr0(I);elseif nr0(I) = n k = nl0(I); end if k = 0Ink0(n, k) = (Vf0(n) - Vf0(k)/ZB0(I); else , end else , end end for I = 1:nbr if nl(I) = n | nr(I) = n if nl(I) =n k = nr(I);elseif nr(I) = n k = nl(I); end if k = 0Inkabc = sctm*Ink0(n, k); Ink1(n, k); Ink2(n, k);Inkabcm = abs(Inkabc); th=angle(I

46、nkabc);if real(Inkabc(1) > 0fprintf('%7g' , n), fprintf('%10g' , k),fprintf(' %11.4f', abs(Inkabc(1),fprintf(' %11.4f'abs(Inkabc(2)fprintf(' %11.4fn', abs(Inkabc(3)elseifreal(Inkabc(1) =0 & imag(Inkabc(1) < 0fprintf('%7g' , n), fprintf('%

47、10g' , k),fprintf(' %11.4f', abs(Inkabc(1),fprintf(' %11.4f'abs(Inkabc(2)fprintf(' %11.4fn', abs(Inkabc(3)else ,endelse , endelse , endendif n=nffprintf(' %11.4f', Ifabcm(1),fprintf(fprintf(' %11.4fn', Ifabcm(3)else ,endendfprintf( '%7g' ,n), fprin

48、tf(' F'),' %11.4f', Ifabcm(2)N' )=1 & strcmp(resp,resp=0;while strcmp(resp,'y' )=1 & strcmp(resp,n' )=1 & strcmp(resp, 'Y' )=1Enter ''y''or ''n''withinsingle'N' )=1& strcmp(resp,'y' )=1), endresp

49、= input( 'Another fault location? quote -> ');if strcmp(resp, 'n' )=1 & strcmp(resp, & strcmp(resp, 'Y' )=1fprintf( 'n Incorrect reply, try again nnendif resp ='y'| resp ='Y'nf = 999;else ff = 0; endend(2) llfault.m 程序代码：% The program llfault is

50、 designed for the line-to-line% fault analysis of a power system network. The program requires % the positive- and negative-sequence bus impedance matrices, % Zbus1, and Zbus2.The bus impedances matrices may be defined % by the user, obtained by the inversion of Ybus or it may be% determined either

51、from the function Zbus = zbuild(zdata)% or the function Zbus = zbuildpi(linedata, gendata, yload).% The program prompts the user to enter the faulted bus number% and the fault impedance Zf. The prefault bus voltages are% defined by the reserved Vector V. The array V may be defined or% it is returned

52、 from the power flow programs lfgauss, lfnewton,% decouple or perturb. If V does not exist the prefault bus voltages% are automatically set to 1.0 per unit. The program obtains the% total fault current, bus voltages and line currents during the fault.function llfault(zdata1, Zbus1, zdata2, Zbus2, V)

53、if exist( 'zdata2' ) = 1zdata2=zdata1;else , endif exist( 'Zbus2' ) = 1Zbus2=Zbus1; else , endnl = zdata1(:,1); nr = zdata1(:,2);R1 = zdata1(:,3); X1 = zdata1(:,4);R2 = zdata2(:,3); X2 = zdata2(:,4);ZB1 = R1 + j*X1; ZB2 = R2 + j*X2;nbr=length(zdata1(:,1); nbus = max(max(nl), max(nr);

54、if exist( 'V' ) = 1if length(V) = nbusV0 = V; else , endelse , V0 = ones(nbus, 1) + j*zeros(nbus, 1); endfprintf( 'nLine-to-line fault analysis n')ff = 999; while ff > 0nf = input( 'Enter Faulted Bus No. -> ');rtn=isempty(nf); if rtn=1; nf=-1;endwhile nf <= 0 | nf &g

55、t; nbusfprintf( 'Faulted bus No. must be between 1 & %g n', nbus)nf = input('Enter Faulted Bus No. -> ');rtn=isempty(nf); if rtn=1; nf=-1;endendrtz=1;while rtz=1fprintf( 'nEnter Fault Impedance Zf = R + j*X in ')Zf = input( rtz=isempty(Zf);'complex form (for bolted

56、 fault enter 0). Zf = ');endfprintf( ' n')fprintf( 'Line-to-line fault at bus No. %gn', nf)a =cos(2*pi/3)+j*sin(2*pi/3);sctm = 1 1 1; 1 aA2 a; 1 a aA2;Ia0=0;Ia1 = V0(nf)/(Zbus1(nf,nf)+Zbus2(nf, nf)+Zf); Ia2=-Ia1;I012=Ia0; Ia1; Ia2;Ifabc = sctm*I012;Ifabcm = abs(Ifabc);fprintf( &#

57、39;Total fault current = %9.4f per unitnn', Ifabcm(2)fprintf( 'Bus Voltages during the fault in per unit nn')fprintf( ' Bus Voltage Magnituden')fprintf( ' No. Phase a Phase b Phase c n') for n = 1:nbusVf0(n)= 0;Vf1(n)= V0(n) - Zbus1(n, nf)*Ia1;Vf2(n)= 0 - Zbus2(n, nf)*Ia2

58、;Vabc = sctm*Vf0(n); Vf1(n); Vf2(n);Va(n)=Vabc(1); Vb(n)=Vabc(2); Vc(n)=Vabc(3);fprintf(' %5g',n)fprintf( ' %11.4f', abs(Va(n),fprintf(' %11.4f', abs(Vb(n)fprintf(' %11.4fn', abs(Vc(n)end fprintf( ' n')fprintf( 'Line currents for fault at bus No. %gnn', nf)fprintf( ' From To Line Current Magnituden')fprintf( ' Bus Bus Phase a Phas