電力系統(tǒng)短路故障分析的MATLAB輔助程序設計-短路計算程序

電力系統(tǒng)短路故障分析的MATLAB輔助程序設計-短路計算程序電力系統(tǒng)短路故障分析的MATLAB輔助程序設計-短路計算程序電力系統(tǒng)短路故障分析的MATLAB輔助程序設計-短路計算程序xxx公司電力系統(tǒng)短路故障分析的MATLAB輔助程序設計-短路計算程序文件編號：文件日期：修訂次數(shù)：第1.0次更改批準審核制定方案設計，管理制度電力系統(tǒng)短路故障分析的MATLAB輔助程序設計電力系統(tǒng)短路故障可分為三相對稱短路故障(three-phasebalancedfaults)和不對稱短路故障(unbalancedfaults)。不對稱短路故障又分為單相接地短路故障(singleline-to-groundfault)、兩相短路故障(line-to-linefault)以及兩相接地短路故障(doubleline-to-groundfault)。根據(jù)故障分析結(jié)果可以對繼電保護裝置、自動裝置進行整定計算，我們可以建立算法來形成節(jié)點阻抗矩陣,利用節(jié)點阻抗矩陣來計算短路故障情況下的節(jié)點電壓和線路電流。一、三相對稱短路故障進行三相短路計算需要兩個程序：zbuild/zbuildpi和symfault程序，zbuild、zbuildpi用來在MATLAB中形成節(jié)點阻抗矩陣，symfault用來計算三相對稱故障。Zbus=zbuild(zdata)這里的參數(shù)zdata是一個(e×4)階矩陣，e是拓撲圖的總支路數(shù)目。第一列和第二列為元素兩端的節(jié)點編號，第三列和第四列分別是線路的電阻、電抗的標幺值。連接在0節(jié)點和發(fā)電機節(jié)點之間的發(fā)電機阻抗可能是次暫態(tài)電抗、暫態(tài)電抗或同步電抗，而且這個矩陣中還包含并聯(lián)電抗器和負荷阻抗。Zbus=zbuildpi(linedata,gendata,yload)這個函數(shù)與潮流計算程序是相容的，第一個參數(shù)linedata與潮流計算程序中的文件是一致的。第一列和第二列為節(jié)點編號；第三列到第五列分別是線路的電阻、電抗以及1/2線路電納值，這三項都為在統(tǒng)一基準容量下的標幺值；最后一列是變壓器分接頭位置，對線路來說，必須輸入1；線路無輸入順序。發(fā)電機參數(shù)不包含在Linedata參數(shù)中，而是包含在第二個參數(shù)gendata中，gendata是一個g×4階矩陣,g是發(fā)電機總數(shù)。第一列和第二列為0節(jié)點、發(fā)電機節(jié)點編號，第三列和第四列為發(fā)電機的暫態(tài)電阻和暫態(tài)電抗。最后一個參數(shù)yload是可選擇的，這個矩陣有兩列，第一列為節(jié)點編號，第二列為復數(shù)導納值，yload可以由潮流程序lfguss，lfnewton或者decouple自動生成。Zbuild和zbuildpi兩個函數(shù)可以通過建立算法求出節(jié)點阻抗矩陣。首先添加所有與參考節(jié)點相連的樹支，然后添加其余的樹支，最后添加共軛連支。程序symfault(zdata,Zbus,V)用來計算三相對稱故障，程序要求輸入zdata和Zbus兩個矩陣，第三個參數(shù)V是可選的。如果V不存在，程序?qū)⒛J故障前所有的節(jié)點電壓標幺值為，如果變量V存在，那么V包括節(jié)點編號和復數(shù)電壓值。電壓向量V也可以由潮流計算程序自動生成。當symfault程序運行時，用戶要輸入故障節(jié)點編號和故障阻抗，運行可得到總的故障電流，節(jié)點電壓幅值以及故障情況下的線路電流。在三相短路計算中，zbuild和symfault程序，zbuildpi和symfault程序都可以進行計算，下面是三相短路計算使用的程序代碼：（1）程序代碼：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)=9;X(k)=9;%無窮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%添加一個新總線分支到現(xiàn)有總線上whiletree<nbusforn=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）程序代碼：%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)=;X(k)=;else,endendifnc>4BC=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%%添加一個新總線分支到現(xiàn)有總線上whiletree<nbusforn=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）程序代碼：%此程序用來計算電網(wǎng)的三相對稱故障%計算前需要用戶生成節(jié)點阻抗矩陣%節(jié)點阻抗矩陣可由函數(shù)Zbus=zbuild(zdata)生成%此程序需要用戶按提示輸入短路節(jié)點編號和過度電阻Zf%向量V是可選參數(shù)，包含節(jié)點編號和復數(shù)電壓%V可由潮流程序lfgauss,lfnewton,decouple自動生成.%如果V不存在，程序?qū)⒛J故障前所有的節(jié)點電壓標幺值為%程序可得到總的故障電流,節(jié)點電壓幅值及輸電線路的電流functionsymfaul(zdata,Zbus,V)nl=zdata(:,1);nr=zdata(:,2);R=zdata(:,3);X=zdata(:,4);nc=length(zdata(1,:));ifnc>4BC=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;whileff>0nf=input('EnterFaultedBusNo.->');rtn=isempty(nf);ifrtn==1;nf=-1;endwhilenf<=0|nf>nbusfprintf('FaultedbusNo.mustbebetween1&%g\n',nbus)nf=input('EnterFaultedBusNo.->');rtn=isempty(nf);ifrtn==1;nf=-1;endendrtz=1;whilertz==1fprintf('\nEnterFaultImpedanceZf=R+j*Xin')Zf=input('complexform(forboltedfaultenter0).Zf=');rtz=isempty(Zf);endfprintf('\n')fprintf('Balancedthree-phasefaultatbusNo.%g\n',nf)If=V0(nf)/(Zf+Zbus(nf,nf));Ifm=abs(If);Ifmang=angle(If)*180/pi;fprintf('Totalfaultcurrent=%perunit\n\n',Ifm)%fprintf('.\n\n',Ifm)fprintf('BusVoltagesduringfaultinperunit\n\n')fprintf('BusVoltageAngle\n')fprintf('No.Magnitudedegrees\n')forn=1:nbusifn==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('%',Vfm),fprintf('%\n',angv)endfprintf('\n')fprintf('LinecurrentsforfaultatbusNo.%g\n\n',nf)fprintf('FromToCurrentAngle\n')fprintf('BusBusMagnitudedegrees\n')forn=1:nbusforI=1:nbrifnl(I)==n|nr(I)==nifnl(I)==nk=nr(I);elseifnr(I)==nk=nl(I);endifk==0Ink=(V0(n)-Vf(n))/ZB(I);Inkm=abs(Ink);th=angle(Ink);ifreal(Ink)>0fprintf('G'),fprintf('%7g',n),fprintf('%',Inkm)fprintf('%\n',th*180/pi)elseifreal(Ink)==0&imag(Ink)<0fprintf('G'),fprintf('%7g',n),fprintf('%',Inkm)fprintf('%\n',th*180/pi)else,endIgn=Ink;elseifk~=0Ink=(Vf(n)-Vf(k))/ZB(I)+BC(I)*Vf(n);Inkm=abs(Ink);th=angle(Ink);ifreal(Ink)>0fprintf('%7g',n),fprintf('%10g',k),fprintf('%',Inkm),fprintf('%\n',th*180/pi)elseifreal(Ink)==0&imag(Ink)<0fprintf('%7g',n),fprintf('%10g',k),fprintf('%',Inkm),fprintf('%\n',th*180/pi)else,endelse,endelse,endendifn==nffprintf('%7g',n),fprintf('F'),fprintf('%',Ifm)fprintf('%\n',Ifmang)else,endendresp=0;whilestrcmp(resp,'n')~=1&strcmp(resp,'N')~=1&strcmp(resp,'y')~=1&strcmp(resp,'Y')~=1resp=input('AnotherfaultlocationEnter''y''or''n''withinsinglequote->');ifstrcmp(resp,'n')~=1&strcmp(resp,'N')~=1&strcmp(resp,'y')~=1&strcmp(resp,'Y')~=1fprintf('\nIncorrectreply,tryagain\n\n'),endendifresp=='y'|resp=='Y'nf=999;elseff=0;endend二、不對稱短路故障不對稱故障包括單相接地短路(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需要輸入正序、零序、負序節(jié)點阻抗矩陣Zbus1、Zbus0和Zbus2，而llfault只需輸入正序和負序節(jié)點阻抗矩陣Zbus1、Zbus2，最后參數(shù)V是可以選擇的。如果用戶未輸入V，程序則默認所有故障前的電壓為。如果用戶輸入?yún)?shù)V，則需要輸入節(jié)點編號以及節(jié)點電壓的復數(shù)值。變量V可由潮流計算程序獲得。零序網(wǎng)絡和正序網(wǎng)絡的節(jié)點阻抗矩陣分別通過函數(shù)Zbus0=zbuild（zdata0）和Zbus1=zbuild（zdata1）獲得。參數(shù)zdata1包含正序網(wǎng)絡阻抗，參數(shù)zdata0包含零序網(wǎng)絡阻抗。參數(shù)zdata1、zdata2、zdata0均為e×4階矩陣，e為元件的個數(shù)。第一列和第二列是節(jié)點編號，第三列和第四列分別是線路電阻和電抗（均為標幺值）。節(jié)點0為參考節(jié)點，發(fā)電機節(jié)點到節(jié)點0的阻抗為發(fā)電機阻抗，可以是發(fā)電機的次暫態(tài)電抗、暫態(tài)電抗或者同步電抗，而且矩陣還包括并聯(lián)電容和負荷阻抗等。負序網(wǎng)絡和正序網(wǎng)絡有同樣的拓撲結(jié)構(gòu)。線路和變壓器的負序阻抗等于正序阻抗，而發(fā)電機的負序阻抗不等于正序阻抗，但在大型電力系統(tǒng)故障分析中通常按相等來計算。零序網(wǎng)絡的拓撲結(jié)構(gòu)不同于正序網(wǎng)絡，零序網(wǎng)絡的構(gòu)建取決于變壓器的繞組連接方式，除了兩側(cè)中性點都接地的Y-Y連接的變壓器，對于其他連接方式的變壓器，一次側(cè)和二次側(cè)的零序網(wǎng)絡都是隔離開的。在這樣的連接方式的零序網(wǎng)絡中,，用inf填寫對應的電阻和電抗列。對于Y-△（中性點接地）連接的變壓器，Y側(cè)的阻抗到節(jié)點0的數(shù)據(jù)需要填寫。當中性點接地阻抗為Xn時，需要填寫為3Xn。在不對稱故障計算時，還需要一個函數(shù)用來求節(jié)點阻抗矩陣，Zbus=zbuildpi（linedata，gendata，yload），函數(shù)與潮流計算程序是相容的。第一個參數(shù)linedata與潮流程序所需要的數(shù)據(jù)相同。第一列和第二列是節(jié)點編號；第三列到第五列分別是線路的電阻、電抗、1/2線路電納（均為標幺值，以指定的MVA為基準功率）；最后一列為變壓器分接頭位置，對于線路，在這一列中必須填寫1。在潮流計算程序中，發(fā)電機電抗不包含在參數(shù)linedata中，而是包含在參數(shù)gendata中。參數(shù)gendata是一個(g×4)階矩陣，其四列分別為參考節(jié)點0、發(fā)電機節(jié)點編號、發(fā)電機電阻和電抗。最后一個參數(shù)yload是可選的，這是一個兩列的矩陣，包括節(jié)點編號和負荷復數(shù)導納值。這些數(shù)據(jù)可以從潮流計算程序lfgauss，lfnewton，decouple執(zhí)行過程中自動生成。程序運行時會提示用戶輸入故障節(jié)點編號和故障阻抗Zf，運行結(jié)果可得到故障電流，節(jié)點電壓和線路電流。不對稱短路故障計算使用的程序如下：（1）程序代碼：%Theprogramlgfaultisdesignedforthesingleline-to-ground%faultanalysisofapowersystemnetwork.Theprogramrequires%thepositive-,negative-andzero-sequencebusimpedancematrices,%Zbus1Zbus2,andbusimpedancesmatricesmaybedefined%bytheuser,obtainedbytheinversionofYbusoritmaybe%determinedeitherfromthefunctionZbus=zbuild(zdata)%orthefunctionZbus=zbuildpi(linedata,gendata,yload).%Theprogrampromptstheusertoenterthefaultedbusnumber%andthefaultimpedanceZf.Theprefaultbusvoltagesare%definedbythereservedVectorV.ThearrayVmaybedefinedor%itisreturnedfromthepowerflowprogramslfgauss,lfnewton,%decoupleorperturb.IfVdoesnotexisttheprefaultbusvoltages%areautomaticallysettoperunit.Theprogramobtainsthe%totalfaultcurrent,busvoltagesandlinecurrentsduringthefault.functionlgfault(zdata0,Zbus0,zdata1,Zbus1,zdata2,Zbus2,V)ifexist('zdata2')~=1zdata2=zdata1;else,endifexist('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);X1=zdata1(:,4);R2=zdata1(:,3);X2=zdata1(:,4);fork=1:nbr0ifR0(k)==inf|X0(k)==infR0(k)=;X0(k)=;else,endendZB1=R1+j*X1;ZB0=R0+j*X0;ZB2=R2+j*X2;ifexist('V')==1iflength(V)==nbusV0=V;else,endelse,V0=ones(nbus,1)+j*zeros(nbus,1);endfprintf('\nLine-to-groundfaultanalysis\n')ff=999;whileff>0nf=input('EnterFaultedBusNo.->');rtn=isempty(nf);ifrtn==1;nf=-1;endwhilenf<=0|nf>nbusfprintf('FaultedbusNo.mustbebetween1&%g\n',nbus)nf=input('EnterFaultedBusNo.->');rtn=isempty(nf);ifrtn==1;nf=-1;endendrtz=1;whilertz==1fprintf('\nEnterFaultImpedanceZf=R+j*Xin')Zf=input('complexform(forboltedfaultenter0).Zf=');rtz=isempty(Zf);endfprintf('\n')fprintf('Singlelineto-groundfaultatbusNo.%g\n',nf)a=cos(2*pi/3)+j*sin(2*pi/3);sctm=[111;1a^2a;1aa^2];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('Totalfaultcurrent=%perunit\n\n',Ifabcm(1))fprintf('BusVoltagesduringthefaultinperunit\n\n')fprintf('Bus-------VoltageMagnitude-------\n')fprintf('No.PhaseaPhasebPhasec\n')forn=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)fprintf('%',abs(Va(n))),fprintf('%',abs(Vb(n)))fprintf('%\n',abs(Vc(n)))endfprintf('\n')fprintf('LinecurrentsforfaultatbusNo.%g\n\n',nf)fprintf('FromTo-----LineCurrentMagnitude----\n')fprintf('BusBusPhaseaPhasebPhasec\n')forn=1:nbusforI=1:nbrifnl(I)==n|nr(I)==nifnl(I)==nk=nr(I);elseifnr(I)==nk=nl(I);endifk~=0Ink1(n,k)=(Vf1(n)-Vf1(k))/ZB1(I);Ink2(n,k)=(Vf2(n)-Vf2(k))/ZB2(I);else,endelse,endendforI=1:nbr0ifnl0(I)==n|nr0(I)==nifnl0(I)==nk=nr0(I);elseifnr0(I)==nk=nl0(I);endifk~=0Ink0(n,k)=(Vf0(n)-Vf0(k))/ZB0(I);else,endelse,endendforI=1:nbrifnl(I)==n|nr(I)==nifnl(I)==nk=nr(I);elseifnr(I)==nk=nl(I);endifk~=0Inkabc=sctm*[Ink0(n,k);Ink1(n,k);Ink2(n,k)];Inkabcm=abs(Inkabc);th=angle(Inkabc);ifreal(Inkabc(1))>0fprintf('%7g',n),fprintf('%10g',k),fprintf('%',abs(Inkabc(1))),fprintf('%',abs(Inkabc(2)))fprintf('%\n',abs(Inkabc(3)))elseifreal(Inkabc(1))==0&imag(Inkabc(1))<0fprintf('%7g',n),fprintf('%10g',k),fprintf('%',abs(Inkabc(1))),fprintf('%',abs(Inkabc(2)))fprintf('%\n',abs(Inkabc(3)))else,endelse,endelse,endendifn==nffprintf('%7g',n),fprintf('F'),fprintf('%',Ifabcm(1)),fprintf('%',Ifabcm(2))fprintf('%\n',Ifabcm(3))else,endendresp=0;whilestrcmp(resp,'n')~=1&strcmp(resp,'N')~=1&strcmp(resp,'y')~=1&strcmp(resp,'Y')~=1resp=input('AnotherfaultlocationEnter''y''or''n''withinsinglequote->');ifstrcmp(resp,'n')~=1&strcmp(resp,'N')~=1&strcmp(resp,'y')~=1&strcmp(resp,'Y')~=1fprintf('\nIncorrectreply,tryagain\n\n'),endendifresp=='y'|resp=='Y'nf=999;elseff=0;endend（2）程序代碼：%Theprogramllfaultisdesignedfortheline-to-line%faultanalysisofapowersystemnetwork.Theprogramrequires%thepositive-andnegative-sequencebusimpedancematrices,%Zbus1,andbusimpedancesmatricesmaybedefined%bytheuser,obtainedbytheinversionofYbusoritmaybe%determinedeitherfromthefunctionZbus=zbuild(zdata)%orthefunctionZbus=zbuildpi(linedata,gendata,yload).%Theprogrampromptstheusertoenterthefaultedbusnumber%andthefaultimpedanceZf.Theprefaultbusvoltagesare%definedbythereservedVectorV.ThearrayVmaybedefinedor%itisreturnedfromthepowerflowprogramslfgauss,lfnewton,%decoupleorperturb.IfVdoesnotexisttheprefaultbusvoltages%areautomaticallysettoperunit.Theprogramobtainsthe%totalfaultcurrent,busvoltagesandlinecurrentsduringthefault.functionllfault(zdata1,Zbus1,zdata2,Zbus2,V)ifexist('zdata2')~=1zdata2=zdata1;else,endifexist('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));ifexist('V')==1iflength(V)==nbusV0=V;else,endelse,V0=ones(nbus,1)+j*zeros(nbus,1);endfprintf('\nLine-to-linefaultanalysis\n')ff=999;whileff>0nf=input('EnterFaultedBusNo.->');rtn=isempty(nf);ifrtn==1;nf=-1;endwhilenf<=0|nf>nbusfprintf('FaultedbusNo.mustbebetween1&%g\n',nbus)nf=input('EnterFaultedBusNo.->');rtn=isempty(nf);ifrtn==1;nf=-1;endendrtz=1;whilertz==1fprintf('\nEnterFaultImpedanceZf=R+j*Xin')Zf=input('complexform(forboltedfaultenter0).Zf=');rtz=isempty(Zf);endfprintf('\n')fprintf('Line-to-linefaultatbusNo.%g\n',nf)a=cos(2*pi/3)+j*sin(2*pi/3);sctm=[111;1a^2a;1aa^2];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('Totalfaultcurrent=%perunit\n\n',Ifabcm(2))fprintf('BusVoltagesduringthefaultinperunit\n\n')fprintf('Bus-------VoltageMagnitude-------\n')fprintf('No.PhaseaPhasebPhasec\n')forn=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('%',abs(Va(n))),fprintf('%',abs(Vb(n)))fprintf('%\n',abs(Vc(n)))endfprintf('\n')fprintf('LinecurrentsforfaultatbusNo.%g\n\n',nf)fprintf('FromTo-----LineCurrentMagnitude----\n')fprintf('BusBusPhaseaPhasebPhasec\n')forn=1:nbusforI=1:nbrifnl(I)==n|nr(I)==nifnl(I)==nk=nr(I);elseifnr(I)==nk=nl(I);endifk~=0Ink0(n,k)=0;Ink1(n,k)=(Vf1(n)-Vf1(k))/ZB1(I);Ink2(n,k)=(Vf2(n)-Vf2(k))/ZB2(I);Inkabc=sctm*[Ink0(n,k);Ink1(n,k);Ink2(n,k)];Inkabcm=abs(Inkabc);th=angle(Inkabc);ifreal(Inkabc(2))<0fprintf('%7g',n),fprintf('%10g',k),fprintf('%',abs(Inkabc(1))),fprintf('%',abs(Inkabc(2)))fprintf('%\n',abs(Inkabc(3)))elseifreal(Inkabc(2))==0&imag(Inkabc(2))>0fprintf('%7g',n),fprintf('%10g',k),fprintf('%',abs(Inkabc(1))),fprintf('%',abs(Inkabc(2)))fprintf('%\n',abs(Inkabc(3)))else,endelse,endelse,endendifn==nffprintf('%7g',n),fprintf('F'),fprintf('%',Ifabcm(1)),fprintf('%',Ifabcm(2))fprintf('%\n',Ifabcm(3))else,endendresp=0;whilestrcmp(resp,'n')~=1&strcmp(resp,'N')~=1&strcmp(resp,'y')~=1&strcmp(resp,'Y')~=1resp=input('AnotherfaultlocationEnter''y''or''n''withinsinglequote->');ifstrcmp(resp,'n')~=1&strcmp(resp,'N')~=1&strcmp(resp,'y')~=1&strcmp(resp,'Y')~=1fprintf('\nIncorrectreply,tryagain\n\n'),endendifresp=='y'|resp=='Y'nf=999;elseff=0;endend（3）程序代碼：%Theprogramdlgfaultisdesignedforthedoubleline-to-ground%faultanalysisofapowersystemnetwork.Theprogramrequires%thepositive-,negative-orzero-sequencebusimpedancematrices,%Zbus1Zbus2,andZbus0.Thebusimpedancesmatricesmaybedefined%bytheuser,obtainedbytheinversionofYbusoritmaybe%determinedeitherfromthefunctionZbus=zbuild(zdata)%orthefunctionZbus=zbuildpi(linedata,gendata,yload).%Theprogrampromptstheusertoenterthefaultedbusnumber%andthefaultimpedanceZf.Theprefaultbusvoltagesare%definedbythereservedVectorV.ThearrayVmaybedefinedor%itisreturnedfromthepowerflowprogramslfgauss,lfnewton,%decoupleorperturb.IfVdoesnotexisttheprefaultbusvoltages%areautomaticallysettoperunit.Theprogramobtainsthe%totalfaultcurrent,busvoltagesandlinecurrentsduringthefault.functiondlgfault(zdata0,Zbus0,zdata1,Zbus1,zdata2,Zbus2,V)ifexist('zdata2')~=1zdata2=zdata1;else,endifexist('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);X1=zdata1(:,4);R2=zdata2(:,3);X2=zdata2(:,4);fork=1:nbr0ifR0(k)==inf|X0(k)==infR0(k)=;X0(k)=9;else,endendZB1=R1+j*X1;ZB0=R0+j*X0;ZB2=R2+j*X2;ifexist('V')==1iflength(V)==nbusV0=V;else,endelse,V0=ones(nbus,1)+j*zeros(nbus,1);endfprintf('\nDoubleline-to-groundfaultanalysis\n')ff=999;whileff>0nf=input('EnterFaultedBusNo.->');rtn=isempty(nf);ifrtn==1;nf=-1;endwhilenf<=0|nf>nbusfprintf('FaultedbusNo.mustbebetween1&%g\n',nbus)nf=input('EnterFaultedBusNo.->');rtn=isempty(nf);ifrtn==1;nf=-1;endendrtz=1;whilertz==1fprintf('\nEnterFaultImpedanceZf=R+j*Xin')Zf=input('complexform(forboltedfaultenter0).Zf=');rtz=isempty(Zf);endfprintf('\n')fprintf('Doubleline-to-groundfaultatbusNo.%g\n',nf)a=cos(2*pi/3)+j*sin(2*pi/3);sctm=[111;1a^2a;1aa^2];Z11=Zbus2(nf,nf)*(Zbus0(nf,nf)+3*Zf)/(Zbus2(nf,nf)+Zbus0(nf,nf)+3*Zf);Ia1=V0(nf)/(Zbus1(nf,nf)+Z11);Ia2=-(V0(nf)-Zbus1(nf,nf)*Ia1)/Zbus2(nf,nf);Ia0=-(V0(nf)-Zbus1(nf,nf)*Ia1)/(Zbus0(nf,nf)+3*Zf);I012=[Ia0;Ia1;Ia2];Ifabc=sctm*I012;Ifabcm=abs(Ifabc);Ift=Ifabc(2)+Ifabc(3);Iftm=abs(Ift);fprintf('Totalfaultcurrent=%perunit\n\n',Iftm)fprintf('BusVoltagesduringthefaultinperunit\n\n')fprintf('Bus-------VoltageMagnitude-------\n')fprintf('No.PhaseaPhasebPhasec\n')forn=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)fprintf('%',abs(Va(n))),fprintf('%',abs(Vb(n)))fprintf('%\n',abs(Vc(n)))endfprintf('\n')fprintf('LinecurrentsforfaultatbusNo.%g\n\n',nf)fprintf('FromTo-----LineCurrentMagnitude----\n')fprintf('BusBusPhaseaPhasebPhasec\n')forn=1:nbusforI=1:nbrifnl(I)==n|nr(I)==nifnl(I)==nk=nr(I);elseifnr(I)==nk=nl(I);endifk~=0Ink1(n,k)=(Vf1(n)-Vf1(k))/ZB1(I);Ink2(n,k)=(Vf2(n)-Vf2(k))/ZB2(I);else,endelse,endendforI=1:nbr0ifnl0(I)==n|nr0(I)==nifnl0(I)==nk=nr0(I);elseifnr0(I)==nk=nl0(I);endifk~=0Ink0(n,k)=(Vf0(n)-Vf0(k))/ZB0(I);else,endelse,endendforI=1:nbrifnl(I)==n|nr(I)==nifnl(I)==nk=nr(I);elseifnr(I)==nk=nl(I);endifk~=0Inkabc=sctm*[Ink0(n,k);Ink1(n,k);Ink2(n,k)];Inkabcm=abs(Inkabc);th=angle(Inkabc);ifreal(Inkabc(2))<0fprintf('%7g',n),fprintf('%10g',k),fprintf('%',abs(Inkabc(1))),fprintf('%',abs(Inkabc(2)))fprintf('%\n',abs(Inkabc(3)))elseifreal(Inkabc(2))==0&imag(Inkabc(2))>0fprintf('%7g',n),fprintf('%10g',k),fprintf('%',abs(Inkabc(1))),fprintf('%',abs(Inkabc(2)))fprintf('%\n',abs(Inkabc(3)))else,endelse,endelse,endendifn==nffprintf('%7g',n),fprintf('F'),fprintf('%',Ifabcm(1)),fprintf('%',Ifabcm(2))fprintf('%\n',Ifabcm(3))else,endendresp=0;whilestrcmp(resp,'n')~=1&strcmp(resp,'N')~=1&strcmp(resp,'y')~=1&strcmp(resp,'Y')~=1resp=input('AnotherfaultlocationEnter''y''or''n''withinsinglequote->');ifstrcmp(resp,'n')~=1&strcmp(resp,'N')~=1&strcmp(resp,'y')~=1&strcmp(resp,'Y')~=1fprintf('\nIncorrectreply,tryagain\n\n'),endendifresp=='y'|resp=='Y'nf=999;elseff=0;endend三、短路故障計算實例一個簡單的11節(jié)點電力系統(tǒng)網(wǎng)絡如圖1所示，發(fā)電機的正序、零序電抗以及中性點接地電抗見表1。變壓器連接方式如圖1所示,節(jié)點11和節(jié)點7之間的變壓器為△-Y連接，接地電抗為。線路和變壓器正序、零序電抗標幺值見表2，基準功率為100MVA。忽略電阻、并聯(lián)電抗和負荷，并假設所有的負序電抗等于正序電抗。試用函數(shù)zbuild求正序和零序網(wǎng)絡的節(jié)點阻抗矩陣,假定故障前所有的節(jié)點電壓都為，用lgfault，llfault,dlgfault和symfault計算下列故障情況下的故障電流、節(jié)點電壓和線路電流。(a)母線2處發(fā)生金屬性單相接地故障;(b)母線4處發(fā)生金屬性兩相短路故障(c)母線6處發(fā)生金屬性兩相接地短路故障;(d)母線8處發(fā)生三相對稱故障;圖111節(jié)點電力系統(tǒng)連接圖圖211節(jié)點電力系統(tǒng)Simulink仿真圖表1發(fā)電機次暫態(tài)電抗發(fā)電機次暫態(tài)電抗節(jié)點編號X1X0Xn11011表2線路和變壓器參數(shù)線路和變壓器參數(shù)節(jié)點編號節(jié)點編號X1X0122325263436464941057687871189在MATLAB的命令窗口輸入正序、零序阻抗數(shù)據(jù)以及求解指令：zdata1=[010100111223252634