電力系統(tǒng)三相短路計算的MATLAB代碼

1、電力系統(tǒng)三相短路計算a main.mclear tim%打開文件dfile,pathname=uigetfile('*.m','Select Data File');if pathname = 0 error(' you must select a valid data file')else lfile =length(dfile); eval(dfile(1:lfile-2);end%定義輸出文件output_file=fopen('output.dat','w');%開始計時tic; %求解節(jié)點導(dǎo)納矩陣，其中Y

2、matrix1是考慮了變比，且支路未近似的導(dǎo)納矩陣；Ymatrix2是近似變比為1，但是支路未近似計算的節(jié)點導(dǎo)納矩陣；Ymatrix3是近似變比為1，采取近似支路參數(shù)1的導(dǎo)納矩陣；Ymatrix4是近似變比為1，采取近似支路參數(shù)2的導(dǎo)納矩陣。Y = Ymatrix2(bus,line);%對故障點進行導(dǎo)納修正fixY = FixY(Y,bus,fault);%求注入電流Iinj = Inode(bus,calcSettings);U = fixYIinj;%得到故障支路與其他支路電流Bcurrent = Ibranch( line,U,fault,Y );%如果發(fā)生支路三相短路，那么對應(yīng)該支路

3、的電流修正為-999999-j999999Ib = ReviseBcurrent( fault,Bcurrent );%結(jié)束計時tim=toc;fprintf('n程序運行結(jié)果');fprintf('n計算完成，共用時%4.4fs,相關(guān)結(jié)果已保存在output.datn',tim);%輸出結(jié)果fprintf_result(output_file, Ib);fprintf_result1(Ib);b FixY.mfunction fixY = FixY( Y,bus,fault )%對形成的導(dǎo)納矩陣進行故障點的修正nb,mb=size(bus);nf,mf= si

4、ze(fault);fixY = Y;%對發(fā)電機節(jié)點導(dǎo)納修正for k=1:nb busType=bus(k,7); if (busType=1) fixY(bus(k,1),bus(k,1) = fixY(bus(k,1),bus(k,1) + 1/1i/bus(k,8); endend%對節(jié)點短路和支路短路的導(dǎo)納矩陣進行修正for k=1:nf nodeI=fault(k,1); nodeJ=fault(k,2); dis=fault(k,3); if (nodeI=0) fixY(nodeJ,nodeJ) = 999999+1i*999999; continue; end if (nod

5、eJ=0) fixY(nodeI,nodeI) = 999999+1i*999999; continue; end if (dis=0)&&(nodeI*nodeJ=0) fixY(nodeI,nodeI) = 999999+1i*999999; continue; end if (dis=1)&&(nodeI*nodeJ=0) fixY(nodeJ,nodeJ) = 999999+1i*999999; continue; end if (dis=1)&&(dis=0)&&(nodeI*nodeJ=0) fixY(nodeI,no

6、deI) = fixY(nodeI,nodeI) - fixY(nodeI,nodeJ)/dis; fixY(nodeJ,nodeJ) = fixY(nodeJ,nodeJ) - fixY(nodeI,nodeJ)/(1-dis); fixY(nodeI,nodeJ)=0; fixY(nodeJ,nodeI)=0; endendend c fprintf_result.mfunction output_args = fprintf_result( output_file, Ib )%將得到的短路電流輸入到輸出文件中n,m=size(Ib); fprintf( output_file, '

7、; No. No. vector of I value of In'); for k=1:n I=Ib(k,1); J=Ib(k,2); I01=real(Ib(k,3); I02=imag(Ib(k,3); I1=Ib(k,4); if(I02>=0) fprintf( output_file, '%3d %3d %10.6f+j%10.6f %10.6f',I,J,I01,I02,I1); end if(I02<0) I02=abs(I02); fprintf( output_file, '%3d %3d %10.6f-j%10.6f %10.6

8、f',I,J,I01,I02,I1); end fprintf( output_file, 'n');endendd fprintf_result1.mfunction output_args = fprintf_result1( Ib )%UNTITLED ÇëÔÚ´Ë´¦ÊäÈ뺯Êý¸ÅÒª n,m=size(Ib); fprintf(' No. No. v

9、ector of I value of In');for k=1:n I=Ib(k,1); J=Ib(k,2); I01=real(Ib(k,3); I02=imag(Ib(k,3); I1=Ib(k,4); if(I02>=0) fprintf('%3d %3d %10.6f+j%10.6f %10.6f',I,J,I01,I02,I1); end if(I02<0) I02=abs(I02); fprintf('%3d %3d %10.6f-j%10.6f %10.6f',I,J,I01,I02,I1); end fprintf('

10、;n');endende Ibranch.mfunction Bcurrent = Ibranch( line,U,fault,Y )%計算短路電流%記錄短路故障參數(shù)，如短路節(jié)點，如為支路短路，記錄距離節(jié)點的距離%此段計算采用的支路參數(shù)未近似，如果計算近似的時候需要修改nl,ml=size(line);Bcurrent=zeros(nl+1,4);faultI=fault(1,1);faultJ=fault(1,2);dis=fault(1,3);faultNode = 0;if(faultI=0) faultNode = faultJ;endif(faultJ=0) faultNod

11、e = faultI;endif(dis=1)&&(faultI*faultJ=0) faultNode = faultJ;endif(dis=0)&&(faultI*faultJ=0) faultNode = faultI;endif(faultNode=0) Bcurrent(nl+1,1) = faultNode; Bcurrent(nl+1,2) = faultNode; Iij = 0; Iij1=0;end%計算非故障支路的短路電流for k=1:nl i=line(k,1); j=line(k,2); Ui=U(i); if j=0 Uj=U(j)

12、; else Uj=0; end if line(k,2)=0 Ym=line(k,5)+1i*line(k,6); Iij=Ui*Ym; Iij1=abs(Iij); end if line(k,2)=0 Zt=line(k,3)+1i*line(k,4); Yt=1/Zt; Ym=line(k,5)+1i*line(k,6); Iij=(Ui-Uj)*Yt+Ui*Ym; Iij1=abs(Iij); end Bcurrent(k,1)=i; Bcurrent(k,2)=j; Bcurrent(k,3)=Iij; Bcurrent(k,4)=Iij1;end%如果為節(jié)點短路，修正短路點的電流

13、大小if(faultNode=0) Bcurrent(nl+1,1) = faultNode; Bcurrent(nl+1,2) = faultNode; Ifault = 0; branchCurrent=0; for k=1:nl I=line(k,1); J=line(k,2); if(I*J=0) continue; end branchCurrent = (U(I)-U(J)/(line(k,3)+1i*line(k,4); if (I=faultNode) Ifault = Ifault - branchCurrent ; else if (J=faultNode) Ifault

14、= Ifault + branchCurrent ; end end end Bcurrent(nl+1,3) = Ifault; Bcurrent(nl+1,4) = abs(Bcurrent(nl+1,3);end %如果為支路短路，修正短路支路的短路電流大小if(dis=0)&&(dis=1)&&(faultI*faultJ=0) Bcurrent(nl+1,1) = faultI; Bcurrent(nl+1,2) = faultJ; Bcurrent(nl+1,3) = U(faultI)*Y(faultI,faultJ)/dis + U(faultJ

15、)*Y(faultI,faultJ)/(1-dis); Bcurrent(nl+1,4) = abs(Bcurrent(nl+1,3);endendf Inode.mfunction Iinj = Inode( bus,calcSettings )%計算節(jié)點注入電流nb,mb=size(bus);Iinj = zeros(nb,1);for k=1:nb busType=bus(k,7); if(calcSettings(1)=1) v = 1; else v = bus(k,2); end%對發(fā)電機節(jié)點電流進行修正 if (busType=1) Iinj(bus(k,1),1) = Iinj

16、(bus(k,1),1) + v/1i/bus(k,8); endendendg ReviseBcurrent.mfunction Ib = ReviseBcurrent( fault,Bcurrent )%如果發(fā)生支路短路，對原來的計算電流進行修正，使該支路短路電流輸出為-999999-j999999clear faultI faultJ disnt,mt=size(Bcurrent);Ib=zeros(nt,mt);faultI=fault(1,1);faultJ=fault(1,2);dis=fault(1,3); for k=1:nt-1 i=Bcurrent(k,1); j=Bcur

17、rent(k,2); Ib(k,:)=Bcurrent(k,:); if (faultI*faultJ=0)&&(dis=1)&&(dis=0)&&(i=faultI)&&(j=faultJ) Ib(k,1)=i; Ib(k,2)=j; Ib(k,3)=-999999-1i*999999; Ib(k,4)=-999999; end if (faultI*faultJ=0)&&(dis=1)&&(dis=0)&&(i=faultJ)&&(j=faultI) Ib(k,1

18、)=i; Ib(k,2)=j; Ib(k,3)=-999999-1i*999999; Ib(k,4)=-999999; end Ib(nt,:)=Bcurrent(nt,:);endh Ymatrix1.mfunction Y = Ymatrix1( bus,line )%考慮變壓器，并且支路參數(shù)不近似的節(jié)點導(dǎo)納矩陣nb,mb=size(bus);nl,ml=size(line);Y=zeros(nb,nb);for k=1:nl I=line(k,1); J=line(k,2); Zt=line(k,3)+1i*line(k,4); Yt=1/Zt; Ym=line(k,5)+1i*line

19、(k,6); K=line(k,7); if (K=0)&&(J=0) Y(I,I)=Y(I,I)+Yt+Ym; Y(J,J)=Y(J,J)+Yt+Ym; Y(I,J)=Y(I,J)-Yt; Y(J,I)=Y(I,J); end if (K=0)&&(J=0) Y(I,I)=Y(I,I)+Ym; end if K>0 Y(I,I)=Y(I,I)+Yt+Ym; Y(J,J)=Y(J,J)+Yt/(K*K); Y(I,J)=Y(I,J)-Yt/K; Y(J,I)=Y(I,J); end if K<0 Y(I,I)=Y(I,I)+Yt+Ym; Y(J,J)

20、=Y(J,J)+K*K*Yt; Y(I,J)=Y(I,J)+K*Yt; Y(J,I)=Y(I,J); endendendi Ymatrix2.mfunction Y = Ymatrix2( bus,line )%考慮變壓器變比近似為1，支路參數(shù)不等效nb,mb=size(bus);nl,ml=size(line);Y=zeros(nb,nb);for k=1:nl I=line(k,1); J=line(k,2); Zt=line(k,3)+1i*line(k,4); Yt=1/Zt; Ym=line(k,5)+1i*line(k,6); if J=0 Y(I,I)=Y(I,I)+Yt+Ym; Y(J,J)=Y(J,J)+Yt+Ym; Y(I,J)=Y(I,J)-Yt; Y(J,I)=Y(I,J); end if J=0 Y(I,I)=Y(I