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

1、电力系统短路故障分析的 MATLAB助程序设计电力系统短路故障可分为三相对称短路故障 (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用来在MATLAB中形成节点阻抗矩阵，symfault用来计算三相对称故障。Zbus=zbuild(zdata) 这里的参数zdata是一个(e x 4)阶矩阵，e是拓扑图 的总支路数目。 第一列和第二列为元素两端的节点编号， 第三列和第四列分别是 线路的电阻、 电抗的标幺值。 连接在 0 节点和发电机节点之间的发电机阻抗可能 是次暂态电抗、 暂态电抗或同步电抗， 而且这个矩阵中还包含并联电抗器和负荷 阻抗。Zbus=zbuildpi(li

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

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

5、 压值。电压向量V也可以由潮流计算程序自动生成。当 symfault程序运行时， 用户要输入故障节点编号和故障阻抗， 运行可得到总的故障电流， 节点电压幅值 以及故障情况下的线路电流。在三相短路计算中， zbuild 和 symfault 程序， zbuildpi 和 symfault 程序 都可 以进行计算， 下面是三相短路计算使用的程序代码：(1)Zbuild.m 程序代码：function Zbus = zbuild(linedata)nl = linedata(:,1); nr = linedata(:,2); R = linedata(:,3); X = linedata(:,4);

6、nbr=length(linedata(:,1); nbus = max(max(nl), max(nr);for k=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

7、);end%newif 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 4 BC = linedata(:,5);for n = 1:nbusyc(n) = 0; nlc(n) = 0; nrc(n) = n;for k = 1:nbrif nl(k) = n | nr(k) = n yc(n) = yc(n) + j*BC(

8、k);else , end endendelseifnc=4 yc= zeros(1, nbr);endnlc=nlc; nrc=nrc; yc=yc.;ZB = R + j*X;if exist( yload ) = 1yload = yload.;yc =yc + yload; else , end m = 0; for n = 1:nbusif abs(yc(n) =0 m=m+1; nlcc(m) = nlc(n); nrcc(m) = nrc(n); zc(m) = 1/yc(n);else , end end nlcc=nlcc; nrcc=nrcc; zc=zc.; nl=nlg

9、; 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 = nr(I);elseif nr(I) = 0 n = nl(I);endif abs(Zbus(n, n) = 0 Zbus(n,n) = ZB(I); tree=tre

10、e+1;else Zbus(n,n) = Zbus(n,n)*ZB(I)/(Zbus(n,n) + ZB(I); endntree(I) = 2;else , endend% 添加一个新总线分支到现有总线上while tree 4BC = zdata(:,5);eIseif nc =4, BC = zeros(Iength(zdata(:,1), 1);endZB = R + j*X;nbr=Iength(zdata(:,1); nbus = max(max(nI), max(nr);if exist( V ) = 1if Iength(V) = nbusV0 = V;else , endel

11、se , V0 = ones(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 nbus, nbus)fprintf(Faulted bus No. must be between 1 & %g nnf = input(Enter Faulted Bus No. - );rtn=isempt

12、y(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(Balanced three-phase fault at bus No. %gn, nf)If = V0(nf)/(Zf + Zbus(nf, nf);Ifm = abs(If); Ifmang=angle(If)*

13、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 Anglen )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;els

14、e , Vf(n) = V0(n) - V0(n)*Zbus(n,nf)/(Zf + Zbus(nf,nf);Vfm = abs(Vf(n); angv=angle(Vf(n)*180/pi; end%13.4fnfprintf( %4g , n), fprintf(%13.4f , Vfm),fprintf( angv) endfprintf(n )fprintf( Line currents for fault at bus No. %gnn , nf) fprintf( From To Current Anglen ) fprintf( Bus Bus Magnitude degrees

15、n ) for n= 1:nbusfor I = 1:nbrif nl(I) = n | nr(I) = nif nl(I) =n k = nr(I); elseifnr(I) = n k = nl(I);endif k=0Ink = (V0(n) - Vf(n)/ZB(I); Inkm = abs(Ink); th=angle(Ink);if real(Ink) 0fprintf( G ), fprintf(%7g,n),fprintf(%12.4fInkm)fprintf(%12.4fn, th*180/pi)elseifreal(Ink) =0 &imag(Ink) 0fprintf(%

16、7g ,n), fprintf(%10g , k),fprintf(%12.4f, Inkm),fprintf( %12.4fn, th*180/pi)elseifreal(Ink) =0 & imag(Ink) );if strcmp(resp, n )=1 & strcmp(resp,N)=1 & strcmp(resp, y )=1& strcmp(resp, Y )=1fprintf(n Incorrect reply, try again nn), endendif resp = y | resp =Ynf = 999;else ff = 0; endend二、不对称短路故障不对称故

17、障包括单相接地短路 (single line-to-ground fault) ，两相短路 (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 和 dlgf

18、ault 需要输入正序、 零序、负序节点阻抗矩阵 Zbus1、 ZbusO和Zbus2,而llfault只需输入正序和负序节点阻抗矩阵 Zbusl、Zbus2,最后参数V是可以选择的。如果用户未输入 V,程序则默认所有故障前的电压为 I.Opu。如果用户输入参数 V,则需要输入节点编号以及节点电压的复数值。变 量 V 可由潮流计算程序获得。零序网络和正序网络的节点阻抗矩阵分别通过函数ZbusO=zbuild (zdataO )和 Zbus1=zbuild (zdata1 )获得。 参数 zdata1 包含正序网络阻抗, 参数 zdataO 包含零序网络阻抗。参数 zdatal、zdata2、z

19、dataO均为ex 4阶矩阵，e为元件 的个数。第一列和第二列是节点编号, 第三列和第四列分别是线路电阻和电抗 (均 为标幺值)。节点 O 为参考节点,发电机节点到节点 O 的阻抗为发电机阻抗,可 以是发电机的次暂态电抗、 暂态电抗或者同步电抗, 而且矩阵还包括并联电容和 负荷阻抗等。负序网络和正序网络有同样的拓扑结构。 线路和变压器的负序阻抗等于正序 阻抗,而发电机的负序阻抗不等于正序阻抗, 但在大型电力系统故障分析中通常 按相等来计算。 零序网络的拓扑结构不同于正序网络, 零序网络的构建取决于变 压器的绕组连接方式，除了两侧中性点都接地的丫-Y连接的变压器，对于其他连 接方式的变压器, 一

20、次侧和二次侧的零序网络都是隔离开的。 在这样的连接方式 的零序网络中，用inf填写对应的电阻和电抗列。对于 Y-(中性点接地)连 接的变压器，丫侧的阻抗到节点0的数据需要填写。当中性点接地阻抗为Xn时，需要填写为 3Xn。在不对称故障计算时， 还需要一个函数用来求节点阻抗矩阵， Zbus=zbuildpi(linedata ，gendata ，yload )， 函数与潮流计算程序是相容的。第一个参数 linedata 与潮流程序所需要的数据相同。第一列和第二列是节点编号；第三列 到第五列分别是线路的电阻、电抗、 1/2 线路电纳(均为标幺值，以指定的 MVA 为基准功率)；最后一列为变压器分接

21、头位置，对于线路，在这一列中必须填写 1。在潮流计算程序中， 发电机电抗不包含在参数 linedata 中，而是包含在参数 gen data中。参数gen data是一个(g x 4)阶矩阵，其四列分别为参考节点0、发 电机节点编号、 发电机电阻和电抗。 最后一个参数 yload 是可选的， 这是一个两 列的矩阵，包括节点编号和负荷复数导纳值。这些数据可以从潮流计算程序 lfgauss ， lfnewton ， decouple 执行过程中自动生成。程序运行时会提示用户输 入故障节点编号和故障阻抗 Zf ，运行结果可得到故障电流，节点电压和线路电 流。不对称短路故障计算使用的程序如下：( 1)

22、 lgfault.m程序代码：% The program lgfault is designed for the single line-to-ground% fault 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

23、by the inversion of Ybus or it may be% determined either 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 reserve

24、d Vector V. The array V may be defined or% it is returned 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 th

25、e fault.function lgfault(zdata0, Zbus0, zdata1, Zbus1, zdata2, Zbus2, V)if exist( zdata2 ) = 1 zdata2=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(z

26、data0(:,1);R0 = zdata0(:,3); X0 = zdata0(:,4);R1 = zdata1(:,3); 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 ) = 1if length(V) = nbusV0 = V;else ,

27、endelse , V0 = ones(nbus, 1) + j*zeros(nbus, 1);endfprintf( 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 nbusfprintf( Faulted bus No. must be between 1 & %g n , nbus) nf = input( Enter Faulted Bus No. - );rt

28、n=isempty(nf);if rtn=1; nf=-1; endendrtz=1;while rtz=1fprintf( nEnter Fault Impedance Zf = R + j*X in );Ia1=Ia0;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 , nf) a =cos(2*pi/3)+j*sin(2*pi/3);sctm = 1 1 1

29、; 1 aA2 a; 1 a aA2;Ia0 = V0(nf)/(Zbus1(nf,nf)+Zbus2(nf,nf)+ Zbus0(nf, nf)+3*Zf);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 during the fault in per unit nn)fprintf( Bus Voltage Magnitude n)fprintf(No

30、. 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);fprintf( %5g ,n), abs(Vb(n), nf)fprintf( %11.4f, abs(Va(n),fprintf( %11.4ffprintf( %11.4

31、fn, abs(Vc(n)endfprintf( n )fprintf(Line currents for fault at bus No. %gnnfprintf( From To Line Current Magnitude nfprintf( BusBus Phase a Phase b Phase c nfor 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);endif k = 0 Ink1(n, k) = (Vf1(n) - Vf1(k)/ZB

32、1(I); Ink2(n, k) = (Vf2(n) - Vf2(k)/ZB2(I);else , endelse , endendfor I = 1:nbr0if nl0(I) = n | nr0(I) = nif nl0(I) =n k = nr0(I);elseifnr0(I) = n k = nl0(I);endif k = 0 Ink0(n, k) = (Vf0(n) - Vf0(k)/ZB0(I);else , endelse , endendfor I = 1:nbrif nl(I) = n | nr(I) = nif nl(I) =n k = nr(I);elseif nr(I

33、) = n k = nl(I); endif k = 0Inkabc = sctm*Ink0(n, k); Ink1(n, k); Ink2(n, k); Inkabcm = abs(Inkabc); th=angle(Inkabc);if real(Inkabc(1) 0fprintf(%7g , n), fprintf( %10g ,k),fprintf( %11.4f, abs(Inkabc(1),fprintf( %11.4fabs(Inkabc(2)fprintf( %11.4fn, abs(Inkabc(3)elseifreal(Inkabc(1) =0 & imag(Inkabc

34、(1) );if strcmp(resp, n)=1 & strcmp(resp,N )=1& strcmp(resp,singley )=1& strcmp(resp,Y )=1), endfprintf( n Incorrect reply, try again nnendif resp = y | resp = Ynf = 999;else ff = 0; end end2) llfault.m 程序代码：% The program llfault is designed for the line-to-line% fault analysis of a power system net

35、work. 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 from the function Zbus = zbuild(zdata) % or the function Zbus = zbuild

36、pi(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 programs lfgauss, lfnewton,% decouple or perturb

37、. 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)if exist( zdata2 ) = 1 zdata2=zdata1;else , endif exist( Zbus2 ) = 1

38、Zbus2=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);if exist( V ) = 1if length(V) = nbusV0 = V;else , endelse , V0 = ones(nbus, 1) + j*zeros

39、(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 nbus, nbus)fprintf(Faulted bus No. must be between 1 & %g nnf = input(Enter Faulted Bus No. - );rtn=isempty(nf);if rtn=1; nf=-1; endendrtz=1;while

40、rtz=1fprintf(nEnter Fault Impedance Zf = R + j*X in Zf = input( complex form (for bolted fault enter 0). Zf = rtz=isempty(Zf);end fprintf( 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)+Z

41、f); Ia2=-Ia1;I012=Ia0; Ia1; Ia2;Ifabc = sctm*I012;Ifabcm = abs(Ifabc);fprintf( Total fault current = %9.4f per unitnn fprintf( Bus Voltages during the fault in per unit nn fprintf( Bus Voltage Magnitude nfprintf( No. Phase a Phase b Phase c nfor n = 1:nbusVf0(n)= 0;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); fprintf( %5g,n)fprintf( %11.4f, abs(Va(n),fprintf( %11.4ffprintf( %11.4fn, abs(Vc(n)endfprintf( n)fprintf( Line currents for fault at bus No. %gnn fprintf( FromToLine Current M