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

电力系统短路故障分析的MATLAB辅助程序设计电力系统短路故障可分为三相对称短路故障(three-phasebalancedfaults)和不对称短路故障(unbalancedfaults)。不对称短路故障又分为单相接地短路故障(singleline-to-groundfault)、两相短路故障(line-to-linefault)以及两相接地短路故障(doubleline-to-groundfault)。根据故障分析结果可以对继电保护装置、自动装置进行整定计算，我们可以建立算法来形成节点阻抗矩阵,利用节点阻抗矩阵来计算短路故障情况下的节点电压和线路电流。一、三相对称短路故障进行三相短路计算需要两个程序：zbuild/zbuildpi和symfault程序，zbuild、zbuildpi用来在MATLAB中形成节点阻抗矩阵，symfault用来计算三相对称故障。Zbus=zbuild(zdata)这里的参数zdata是一个(e×4)阶矩阵，e是拓扑图的总支路数目。第一列和第二列为元素两端的节点编号，第三列和第四列分别是线路的电阻、电抗的标幺值。连接在0节点和发电机节点之间的发电机阻抗可能是次暂态电抗、暂态电抗或同步电抗，而且这个矩阵中还包含并联电抗器和负荷阻抗。Zbus=zbuildpi(linedata,gendata,yload)这个函数与潮流计算程序是相容的，第一个参数linedata与潮流计算程序中的文件是一致的。第一列和第二列为节点编号；第三列到第五列分别是线路的电阻、电抗以及1/2线路电纳值，这三项都为在统一基准容量下的标幺值；最后一列是变压器分接头位置，对线路来说，必须输入1；线路无输入顺序。发电机参数不包含在Linedata参数中，而是包含在第二个参数gendata中，gendata是一个g×4阶矩阵,g是发电机总数。第一列和第二列为0节点、发电机节点编号，第三列和第四列为发电机的暂态电阻和暂态电抗。最后一个参数yload是可选择的，这个矩阵有两列，第一列为节点编号，第二列为复数导纳值，yload可以由潮流程序lfguss，lfnewton或者decouple自动生成。Zbuild和zbuildpi两个函数可以通过建立算法求出节点阻抗矩阵。首先添加所有与参考节点相连的树支，然后添加其余的树支，最后添加共轭连支。程序symfault(zdata,Zbus,V)用来计算三相对称故障，程序要求输入zdata和Zbus两个矩阵，第三个参数V是可选的。如果V不存在，程序将默认故障前所有的节点电压标幺值为1.0，如果变量V存在，那么V包括节点编号和复数电压值。电压向量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);nbr=length(linedata(:,1));nbus=max(max(nl),max(nr));fork=1:nbrifR(k)==inf|X(k)==infR(k)=999999999;X(k)=999999999;%无穷else,endendZB=R+j*X;Zbus=zeros(nbus,nbus);tree=0;%从参考总线0上添加一个分支forI=1:nbrntree(I)=1;ifnl(I)==0|nr(I)==0ifnl(I)==0n=nr(I);elseifnr(I)==0n=nl(I);endifabs(Zbus(n,n))==0Zbus(n,n)=ZB(I);tree=tree+1;%%newelseZbus(n,n)=Zbus(n,n)*ZB(I)/(Zbus(n,n)+ZB(I));endntree(I)=2;else,endend%添加一个新总线分支到现有总线上whiletreenbusforn=1:nbusnadd=1;ifabs(Zbus(n,n))==0forI=1:nbrifnadd==1;ifnl(I)==n|nr(I)==nifnl(I)==nk=nr(I);elseifnr(I)==nk=nl(I);endifabs(Zbus(k,k))~=0form=1:nbusifm~=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%增加两个原有总线间的支路阻抗forn=1:nbusforI=1:nbrifntree(I)==1ifnl(I)==n|nr(I)==nifnl(I)==nk=nr(I);elseifnr(I)==nk=nl(I);endDM=Zbus(n,n)+Zbus(k,k)+ZB(I)-2*Zbus(n,k);forjj=1:nbusAP=Zbus(jj,n)-Zbus(jj,k);forkk=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程序代码：%Thisprogramformsthecomplexbusimpedancematrixbythemethod%ofbuildingalgorithm.Buszeroistakenasreference.%Thisprogramiscompatiblewithpowerflowdata.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(:,1);nr=linedata(:,2);R=linedata(:,3);X=linedata(:,4);nbr=length(linedata(:,1));nbus=max(max(nl),max(nr));nc=length(linedata(1,:));fork=1:nbrifR(k)==inf|X(k)==infR(k)=99999999;X(k)=99999999;else,endendifnc4BC=linedata(:,5);forn=1:nbusyc(n)=0;nlc(n)=0;nrc(n)=n;fork=1:nbrifnl(k)==n|nr(k)==nyc(n)=yc(n)+j*BC(k);else,endendendelseifnc==4yc=zeros(1,nbr);endnlc=nlc';nrc=nrc';yc=yc.';ZB=R+j*X;ifexist('yload')==1yload=yload.';yc=yc+yload;else,endm=0;forn=1:nbusifabs(yc(n))~=0m=m+1;nlcc(m)=nlc(n);nrcc(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=[nlnrreal(ZB)imag(ZB)];nbr=length(nl);Zbus=zeros(nbus,nbus);tree=0;%%从参考总线0上添加一个分支forI=1:nbrntree(I)=1;ifnl(I)==0|nr(I)==0ifnl(I)==0n=nr(I);elseifnr(I)==0n=nl(I);endifabs(Zbus(n,n))==0Zbus(n,n)=ZB(I);tree=tree+1;elseZbus(n,n)=Zbus(n,n)*ZB(I)/(Zbus(n,n)+ZB(I));endntree(I)=2;else,endend%%添加一个新总线分支到现有总线上whiletreenbusforn=1:nbusnadd=1;ifabs(Zbus(n,n))==0forI=1:nbrifnadd==1ifnl(I)==n|nr(I)==nifnl(I)==nk=nr(I);elseifnr(I)==nk=nl(I);endifabs(Zbus(k,k))~=0form=1:nbusifm~=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%增加两个原有总线间的支路阻抗forn=1:nbusforI=1:nbrifntree(I)==1ifnl(I)==n|nr(I)==nifnl(I)==nk=nr(I);elseifnr(I)==nk=nl(I);endDM=Zbus(n,n)+Zbus(k,k)+ZB(I)-2*Zbus(n,k);forjj=1:nbusAP=Zbus(jj,n)-Zbus(jj,k);forkk=1:nbusAT=Zbus(n,kk)-Zbus(k,kk);DELZ(jj,kk)=AP*AT/DM;endendZbus=Zbus-DELZ;ntree(I)=2;else,endelse,endendend（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,:));ifnc4BC=zdata(:,5);elseifnc==4,BC=zeros(length(zdata(:,1)),1);endZB=R+j*X;nbr=length(zdata(:,1));nbus=max(max(nl),max(nr));ifexist('V')==1iflength(V)==nbusV0=V;else,endelse,V0=ones(nbus,1)+j*zeros(nbus,1);endfprintf('\Three-phasebalancedfaultanalysis\n')ff=999;whileff0nf=input('EnterFaultedBusNo.-');rtn=isempty(nf);ifrtn==1;nf=-1;endwhilenf=0|nfnbusfprintf('FaultedbusNo.mustbebetween1&%g\n',nbus)nf=input('EnterFaul