电力系统三相短路计算地MATLAB代码
电力系统三相短路计算的MATLAB代码
![电力系统三相短路计算的MATLAB代码](https://img.taocdn.com/s3/m/20b5ed2c4693daef5ef73dfb.png)
. . . .电力系统三相短路计算a main.mclear tim%打开文件[dfile,pathname]=uigetfile('*.m','Select Data File');if pathname == 0error(' you must select a valid data file')elselfile =length(dfile);eval(dfile(1:lfile-2));end%定义输出文件output_file=fopen('output.dat','w');%开始计时tic;%求解节点导纳矩阵,其中Ymatrix1是考虑了变比,且支路未近似的导纳矩阵;Ymatrix2是近似变比为1,但是支路未近似计算的节点导纳矩阵;Ymatrix3是近似变比为1,采取近似支路参数1的导纳矩阵;Ymatrix4是近似变比为1,采取近似支路参数2的导纳矩阵。
Y = Ymatrix2(bus,line);%对故障点进行导纳修正fixY = FixY(Y,bus,fault);%求注入电流Iinj = Inode(bus,calcSettings);U = fixY\Iinj;%得到故障支路与其他支路电流Bcurrent = Ibranch( line,U,fault,Y );%如果发生支路三相短路,那么对应该支路的电流修正为-999999-j999999Ib = ReviseBcurrent( fault,Bcurrent );%结束计时tim=toc;fprintf('\n程序运行结果');fprintf('\n计算完成,共用时%4.4fs,相关结果已保存在output.dat\n',tim);%输出结果fprintf_result(output_file, Ib);fprintf_result1(Ib);b FixY.mfunction fixY = FixY( Y,bus,fault )%对形成的导纳矩阵进行故障点的修正[nb,mb]=size(bus);[nf,mf]= size(fault);fixY = Y;%对发电机节点导纳修正for k=1:nbbusType=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%对节点短路和支路短路的导纳矩阵进行修正for k=1:nfnodeI=fault(k,1);nodeJ=fault(k,2);dis=fault(k,3);if (nodeI==0)fixY(nodeJ,nodeJ) = 999999+1i*999999;continue;endif (nodeJ==0)fixY(nodeI,nodeI) = 999999+1i*999999;continue;endif (dis==0)&&(nodeI*nodeJ~=0)fixY(nodeI,nodeI) = 999999+1i*999999;continue;endif (dis==1)&&(nodeI*nodeJ~=0)fixY(nodeJ,nodeJ) = 999999+1i*999999;continue;endif (dis~=1)&&(dis~=0)&&(nodeI*nodeJ~=0)fixY(nodeI,nodeI) = 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;endendendc fprintf_result.mfunction [ output_args ] = fprintf_result( output_file, Ib )%将得到的短路电流输入到输出文件中[n,m]=size(Ib);fprintf( output_file, ' No. No. vector of I value of I\n');for k=1:nI=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);endif(I02<0)I02=abs(I02);fprintf( output_file, '%3d %3d %10.6f-j%10.6f %10.6f',I,J,I01,I02,I1);endfprintf( output_file, '\n');endendd fprintf_result1.mfunction [ output_args ] = fprintf_result1( Ib )%UNTITLED ÇëÔÚ´Ë´¦ÊäÈ뺯Êý¸ÅÒª[n,m]=size(Ib);fprintf(' No. No. vector of I value of I\n');for k=1:nI=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);endif(I02<0)I02=abs(I02);fprintf('%3d %3d %10.6f-j%10.6f %10.6f',I,J,I01,I02,I1);endfprintf('\n');endende Ibranch.mfunction Bcurrent = Ibranch( line,U,fault,Y )%计算短路电流%记录短路故障参数,如短路节点,如为支路短路,记录距离节点的距离%此段计算采用的支路参数未近似,如果计算近似的时候需要修改[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)faultNode = 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:nli=line(k,1);j=line(k,2);Ui=U(i);if j~=0Uj=U(j);elseUj=0;endif line(k,2)==0Ym=line(k,5)+1i*line(k,6);Iij=Ui*Ym;Iij1=abs(Iij);endif line(k,2)~=0Zt=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);endBcurrent(k,1)=i;Bcurrent(k,2)=j;Bcurrent(k,3)=Iij;Bcurrent(k,4)=Iij1;end%如果为节点短路,修正短路点的电流大小if(faultNode~=0)Bcurrent(nl+1,1) = faultNode;Bcurrent(nl+1,2) = faultNode;Ifault = 0;branchCurrent=0;for k=1:nlI=line(k,1);J=line(k,2);if(I*J==0)continue;endbranchCurrent = (U(I)-U(J))/(line(k,3)+1i*line(k,4));if (I==faultNode)Ifault = Ifault - branchCurrent ;elseif (J==faultNode)Ifault = Ifault + branchCurrent ;endendendBcurrent(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)*Y(faultI,faultJ)/(1-dis);Bcurrent(nl+1,4) = abs(Bcurrent(nl+1,3));endendf Inode.mfunction Iinj = Inode( bus,calcSettings )%计算节点注入电流[nb,mb]=size(bus);Iinj = zeros(nb,1);for k=1:nbbusType=bus(k,7);if(calcSettings(1)==1)v = 1;elsev = bus(k,2);end%对发电机节点电流进行修正if (busType==1)Iinj(bus(k,1),1) = Iinj(bus(k,1),1) + v/1i/bus(k,8);endendendg ReviseBcurrent.mfunction Ib = ReviseBcurrent( fault,Bcurrent )%如果发生支路短路,对原来的计算电流进行修正,使该支路短路电流输出为-999999-j999999 clear faultIfaultJdis[nt,mt]=size(Bcurrent);Ib=zeros(nt,mt);faultI=fault(1,1);faultJ=fault(1,2);dis=fault(1,3);for k=1:nt-1i=Bcurrent(k,1);j=Bcurrent(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;endif (faultI*faultJ~=0)&&(dis~=1)&&(dis~=0)&&(i==faultJ)&&(j==faultI)Ib(k,1)=i;Ib(k,2)=j;Ib(k,3)=-999999-1i*999999;Ib(k,4)=-999999;endIb(nt,:)=Bcurrent(nt,:);endh Ymatrix1.mfunction Y = Ymatrix1( bus,line )%考虑变压器,并且支路参数不近似的节点导纳矩阵[nb,mb]=size(bus);[nl,ml]=size(line);Y=zeros(nb,nb);for k=1:nlI=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);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);endif (K==0)&&(J==0)Y(I,I)=Y(I,I)+Ym;endif K>0Y(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);endif K<0Y(I,I)=Y(I,I)+Yt+Ym;Y(J,J)=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,支路参数不等效[nb,mb]=size(bus);[nl,ml]=size(line);Y=zeros(nb,nb);for k=1:nlI=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~=0Y(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);endif J==0Y(I,I)=Y(I,I)+Ym;endendendj Ymatrix3.mfunction Y = Ymatrix3( bus,line )%考虑变压器变比为1,采用支路参数近似1[nb,mb]=size(bus);[nl,ml]=size(line);Y=zeros(nb,nb);for k=1:nlI=line(k,1);J=line(k,2);Zt=line(k,3)+1i*line(k,4);Yt=imag(1/Zt);Ym=imag(line(k,5)+1i*line(k,6));if J~=0Y(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);endif J==0Y(I,I)=Y(I,I)+Ym;endendendk Ymatrix4.mfunction Y = Ymatrix4( bus,line )%变压器变比近似为1,采用支路等效参数2 [nb,mb]=size(bus);[nl,ml]=size(line);Y=zeros(nb,nb);for k=1:nlI=line(k,1);J=line(k,2);Zt=1i*line(k,4);Yt=1/Zt;Ym=1i*line(k,6);if J~=0Y(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);endif J==0Y(I,I)=Y(I,I)+Ym;endendend。
MATLAB实验:同步发电机三相短路响应过程仿真
![MATLAB实验:同步发电机三相短路响应过程仿真](https://img.taocdn.com/s3/m/c740111a3a3567ec102de2bd960590c69ec3d800.png)
MATLAB实验:同步发电机三相短路响
应过程仿真
简介
本实验通过使用MATLAB软件对同步发电机的三相短路响应过程进行仿真。
同步发电机是发电厂中常见的发电设备,了解其短路响应过程对于保证电网的稳定性很重要。
实验过程
1. 确定同步发电机的参数,包括额定功率、额定电压、额定频率、定子、转子等参数。
2. 通过在MATLAB中建立电力系统模型,将同步发电机和其他电网元件(例如传输线、发电机组)连接起来。
3. 设计相应的控制策略,实现短路故障的切除、保护与恢复。
4. 运行仿真程序,模拟同步发电机在发生三相短路故障时的响应过程。
5. 分析仿真结果,包括电流、电压、转速等参数的变化情况,以及系统的稳定性和安全性。
实验目标
通过进行该实验,我们可以了解同步发电机在发生三相短路故障时的响应过程,包括电流和电压的变化情况。
这有助于我们更好地了解电力系统的稳定性和安全性,并为实际电力系统中的故障分析和保护设计提供参考。
结论
通过对同步发电机三相短路响应过程的仿真实验,我们可以获得电流和电压的变化情况,并进一步分析电力系统的稳定性和安全性。
这对于电力系统的运行和维护至关重要,有助于提高电网的可靠性和安全性。
【注意】本文档中的内容仅供参考,详细的实验步骤和具体参数需根据实际情况进行确定和调整。
电力系统三相短路计算的MATLAB代码.doc
![电力系统三相短路计算的MATLAB代码.doc](https://img.taocdn.com/s3/m/6867c297192e45361066f5f3.png)
电力系统三相短路计算a main.mclear tim%打开文件[dfile,pathname]=uigetfile('*.m' , 'Select Data File');if pathname == 0error( ' you must select a valid data file')elselfile =length(dfile);eval(dfile(1:lfile-2));end%定义输出文件output_file=fopen('output.dat', 'w' );%开始计时tic;%求解节点导纳矩阵,其中Ymatrix1 是考虑了变比,且支路未近似的导纳矩阵;Ymatrix2 是近似变比为 1,但是支路未近似计算的节点导纳矩阵;Ymatrix3 是近似变比为 1,采取近似支路参数 1的导纳矩阵; Ymatrix4 是近似变比为 1,采取近似支路参数2的导纳矩阵。
Y = Ymatrix2(bus,line);%对故障点进行导纳修正fixY = FixY(Y,bus,fault);%求注入电流Iinj = Inode(bus,calcSettings);U = fixY\Iinj;%得到故障支路与其他支路电流Bcurrent = Ibranch( line,U,fault,Y );%如果发生支路三相短路,那么对应该支路的电流修正为-999999-j999999Ib = ReviseBcurrent( fault,Bcurrent );%结束计时tim=toc;fprintf( '\n 程序运行结果 ' );fprintf( '\n 计算完成,共用时 %4.4fs, 相关结果已保存在 output.dat\n',tim);%输出结果fprintf_result(output_file, Ib);fprintf_result1(Ib);b FixY.mfunction fixY = FixY( Y,bus,fault )%对形成的导纳矩阵进行故障点的修正[nb,mb]=size(bus);[nf,mf]= size(fault);fixY = Y;%对发电机节点导纳修正for k=1:nbbusType=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%对节点短路和支路短路的导纳矩阵进行修正for k=1:nfnodeI=fault(k,1);nodeJ=fault(k,2);dis=fault(k,3);if (nodeI==0)fixY(nodeJ,nodeJ) = 999999+1i*999999;continue ;endif(nodeJ==0)fixY(nodeI,nodeI) = 999999+1i*999999;continue ;if(dis==0)&&(nodeI*nodeJ~=0)fixY(nodeI,nodeI) = 999999+1i*999999;continue ;endif(dis==1)&&(nodeI*nodeJ~=0)fixY(nodeJ,nodeJ) = 999999+1i*999999;continue ;endif(dis~=1)&&(dis~=0)&&(nodeI*nodeJ~=0)fixY(nodeI,nodeI) = 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;endendc fprintf_result.mfunction [ output_args ] = fprintf_result( output_file, Ib )%将得到的短路电流输入到输出文件中[n,m]=size(Ib);fprintf(output_file,' No. No. vector of I value of I\n' );for k=1:nI=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);endif (I02<0)I02=abs(I02);fprintf( output_file,'%3d %3d %10.6f-j%10.6f %10.6f',I,J,I01,I02,I1);endfprintf( output_file,'\n' );endendd fprintf_result1.mfunction [ output_args ] = fprintf_result1( Ib )%UNTITLED? ? ? ú′? ′|ê? è? oˉêy ? ? òa[n,m]=size(Ib);fprintf( ' No. No. vector of I value of I\n' );for k=1:nI=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);endif (I02<0)I02=abs(I02);fprintf( '%3d %3d %10.6f-j%10.6f %10.6f' ,I,J,I01,I02,I1);endfprintf('\n' );endende Ibranch.mfunction Bcurrent = Ibranch( line,U,fault,Y )%计算短路电流%记录短路故障参数,如短路节点,如为支路短路,记录距离节点的距离%此段计算采用的支路参数未近似,如果计算近似的时候需要修改[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)faultNode = 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:nli=line(k,1);j=line(k,2);Ui=U(i);if j~=0Uj=U(j);elseUj=0;endif line(k,2)==0Ym=line(k,5)+1i*line(k,6);Iij=Ui*Ym;Iij1=abs(Iij);endif line(k,2)~=0Zt=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);endBcurrent(k,1)=i;Bcurrent(k,2)=j;Bcurrent(k,3)=Iij;Bcurrent(k,4)=Iij1;end%如果为节点短路,修正短路点的电流大小if (faultNode~=0)Bcurrent(nl+1,1) = faultNode;Bcurrent(nl+1,2) = faultNode;Ifault = 0;branchCurrent=0;for k=1:nlI=line(k,1);J=line(k,2);if (I*J==0)continue ;endbranchCurrent = (U(I)-U(J))/(line(k,3)+1i*line(k,4));if(I==faultNode)Ifault = Ifault - branchCurrent ;else if (J==faultNode)Ifault = Ifault + branchCurrent ;endendendBcurrent(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)*Y(faultI,faultJ)/(1-dis);Bcurrent(nl+1,4) = abs(Bcurrent(nl+1,3));endendf Inode.mfunction Iinj = Inode( bus,calcSettings )%计算节点注入电流[nb,mb]=size(bus);Iinj = zeros(nb,1);for k=1:nbbusType=bus(k,7);if (calcSettings(1)==1)v = 1;elsev = bus(k,2);end%对发电机节点电流进行修正if (busType==1)Iinj(bus(k,1),1) = Iinj(bus(k,1),1) + v/1i/bus(k,8);endendendg ReviseBcurrent.mfunction Ib = ReviseBcurrent( fault,Bcurrent )%如果发生支路短路,对原来的计算电流进行修正,使该支路短路电流输出为 -999999-j999999clear faultI faultJ dis[nt,mt]=size(Bcurrent);Ib=zeros(nt,mt);faultI=fault(1,1);faultJ=fault(1,2);dis=fault(1,3);for k=1:nt-1i=Bcurrent(k,1);j=Bcurrent(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;endif(faultI*faultJ~=0)&&(dis~=1)&&(dis~=0)&&(i==faultJ)&&(j==faultI) Ib(k,1)=i;Ib(k,2)=j;Ib(k,3)=-999999-1i*999999;Ib(k,4)=-999999;endIb(nt,:)=Bcurrent(nt,:);endh Ymatrix1.mfunction Y = Ymatrix1( bus,line )%考虑变压器,并且支路参数不近似的节点导纳矩阵[nb,mb]=size(bus);[nl,ml]=size(line);Y=zeros(nb,nb);for k=1:nlI=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);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);endif (K==0)&&(J==0) Y(I,I)=Y(I,I)+Ym;endif K>0Y(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);endif K<0Y(I,I)=Y(I,I)+Yt+Ym;Y(J,J)=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,支路参数不等效[nb,mb]=size(bus);[nl,ml]=size(line);Y=zeros(nb,nb);for k=1:nlI=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~=0Y(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);endif J==0Y(I,I)=Y(I,I)+Ym;endendendj Ymatrix3.mfunction Y = Ymatrix3( bus,line )%考虑变压器变比为 1,采用支路参数近似 1 [nb,mb]=size(bus);[nl,ml]=size(line);Y=zeros(nb,nb);for k=1:nlI=line(k,1);J=line(k,2);Zt=line(k,3)+1i*line(k,4);Yt=imag(1/Zt);Ym=imag(line(k,5)+1i*line(k,6));if J~=0Y(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);endif J==0Y(I,I)=Y(I,I)+Ym;endendendk Ymatrix4.mfunction Y = Ymatrix4( bus,line )%变压器变比近似为 1,采用支路等效参数 2 [nb,mb]=size(bus);[nl,ml]=size(line);Y=zeros(nb,nb);for k=1:nlI=line(k,1);J=line(k,2);Zt=1i*line(k,4);Yt=1/Zt;Ym=1i*line(k,6);if J~=0Y(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);endif J==0Y(I,I)=Y(I,I)+Ym;endend end。
电力系统三相短路的编程计算
![电力系统三相短路的编程计算](https://img.taocdn.com/s3/m/58c0415cf5335a8103d2202c.png)
电力系统分析课程专题报告学生姓名:班级:学号:指导教师:所在单位:提交日期:评分电力系统三相短路的编程计算某某某(学院, )摘要:在电力系统中,三相短路故障造成的危害是最大的,发生的几率也最高,故短路计算对电力系统的稳定运行具有十分重要的意义。
作为电力系统三大计算之一,计算三相短路故障发生时的短路电流、各节点电压、各支路电流是短路计算的基本内容。
在电力系统短路电流的工程计算中,由于快速继电保护的应用,最重要的是计算短路电流基频交流分量的初始值,即次暂态电流I ''。
在给定电源电势时,实际上就是稳态交流电路的求解。
本文基于教材例3-2应用MATLAB 编程计算三相短路故障的电流电压情况,并于例3-3进行进一步的验证和完善。
关键词:三相短路;MATLAB中图分类号:TM 713 文献标识码:A0 引言在电力系统的四种短路类型中,三相短路是其中最严重的,其短路电流可达数万安以至十几万安,随之产生的热效应和电动力效应将使电气设备遭受严重破坏。
因此,计算短路电流的主要应用目的是电力系统设计中的电气设备选择,短路计算已成为电力系统运行分析、设计计算的重要环节。
实际电力系统短路电流交流分量初始值的计算,小型系统可以手算,而对于结构复杂的大型系统,短路电流计算量较大,用计算机进行辅助计算成为大势所趋。
1 解析法求三相短路电流1.1 参数说明(1) 为了元件参数标幺值计算方便,取基准容量B S 为A MV ⋅60,可设任意值,但必须唯一值参与计算。
(2) 取基准电压B U 为平均额定电压AV U ,基于例3-2中系统的额定电压等级有kv 10、kv 110,平均额定电压分别为kv 115、kv 5.10,平均额定电压与线路额定电压相差5%,为简化计算,故取平均额定电压。
(3) I ''为次暂态短路电流有效值,短路电流周期分量的初值等于时间0=t 时的有效值。
满足产生最大短路电流的三个条件下的最大次暂态短路电流作为计算依据。
基于matlab的电力系统短路电流计算
![基于matlab的电力系统短路电流计算](https://img.taocdn.com/s3/m/272d42c49ec3d5bbfd0a742e.png)
湖北民族学院信息工程学院题目: 基于matlab的电力系统短路电流计算专业:电气工程及其自动化班级: 0308407学号: *********学生姓名:指导教师:2011年6 月1 日信息工程学院课程设计任务书年月日信息工程学院课程设计成绩评定表摘要随着电力工业的发展,电力系统的规模越来越大,在这种情况下,许多大型的电力科研实验很难进行,尤其是电力系统中对设备和人员等危害最大的事故故障,尤其是短路故障,而在分析解决事故故障时要不断的实验,在现实设备中很难实现,一是实际的条件难以满足;二是从系统的安全角度来讲也是不允许进行实验的。
考虑这两种情况,寻求一种最接近于电力系统实际运行状况的数字仿真工具十分重要,而MATLAB软件中的SIMULINK是用来对动态系统进行建模、仿真和分析的集成开发环境,是结合了框图界面和交互仿真能力的非线性动态系统仿真工具,为解决具体的工程问题提供了更为快速、准确和简洁的途径。
关键词:短路电流计算,MATLAB,仿真AbstractAlong with the development of the electric power industry, the scale of the power system is more and more big, in this case, many large power research is difficult to, especially in the power system, equipment and personnel to the harm such as the biggest accident, especially fault fault location, and on the analysis of the accident to solve the fault of the experiment, in the reality constantly in equipment, it is difficult to accomplish a is practical conditions to meet; The security of the system from the perspective is not allowed in the experiment. Consider the two kinds of circumstances, for one of the most close to power system actual the operation condition of digital simulation tool is very important, and MATLAB software SIMULINK is used for the dynamic system modeling, simulation and analysis of the integrated development environment, is combined with the block diagram interface and interactive simulation of nonlinear dynamic system ability of simulation tools, to solve the specific engineering problem, provides a more rapid, accurate and simple way.Keywords: short-circuit current calculation, MATLAB, the simulatio目录摘要 (4)1 概述 (6)1.1短路产生的原因 (6)1.2短路的危害 (6)1.3短路故障分析的内容和目的 (6)1.4防范短路电流的措施 (6)2短路计算 (8)2.1简单不对称故障的分析 (8)2.2短路电流的计算过程 (8)2.2.1选择基准值 (10)2.2.2计算系统各元件阻抗标幺值 (10)2.2.3求短路电流的周期分量及冲击电流 (10)3用MA TLAB计算短路电流 (12)3.1MA TLAB简介 (12)3.1.1MATLAB应用 (12)3.2Simulink简介 (12)3.2.1Simulink功能 (13)3.2.2Simulink特点 (13)3.3用MATLAB计算短路电流的实现 (14)3.3.1短路计算内容概述 (14)3.3.2电力系统短路电流计算仿真运行 (15)3.3.3仿真参数设置 (15)3.3.4仿真结果 (16)4 总结 (18)参考文献 (19)1 概述1.1短路产生的原因(1).元件损坏例如绝缘材料的自然老化,设计,安装维护不良所带来的设备缺陷发展成短路等,(2).气象条件恶化例如雷击造成的闪络放电或避雷器动作,架空线路由于大风或导线覆冰引起电杆倒塌(3).违规操作,例如运行人员带负荷拉刀闸,线路或设备检修后未拆除接地线就加上电压等;(4). 其他,例如挖沟损伤电缆,鸟兽跨接在裸露的载流部分等。
电力系统短路故障分析的MATLAB辅助程序设计短路计算程序
![电力系统短路故障分析的MATLAB辅助程序设计短路计算程序](https://img.taocdn.com/s3/m/671ca63530b765ce0508763231126edb6f1a7608.png)
电力系统短路故障分析的MATLAB辅助程序设计短路计算程序电力系统短路故障分析是电力系统设计和运行过程中非常重要的一环。
短路故障会导致电力系统各个部分的电压、电流和功率的突然变化,对设备的保护和稳定运行产生不利影响甚至引起事故。
因此,进行短路计算和故障分析非常必要。
MATLAB是一种功能强大的数值计算和数据可视化工具,对于电力系统短路计算和故障分析也可以发挥重要的作用。
下面将介绍如何使用MATLAB设计一个简单的电力系统短路计算程序。
首先,我们需要建立一个电力系统的模型。
电力系统可以用图模型表示,其中节点表示发电机、变压器、负荷等设备,边表示导线、变压器等电力连接。
我们可以使用MATLAB中的图模型工具箱创建电力系统模型,并且设置各个节点和边的属性,例如电压、电流、阻抗等。
然后,我们需要编写短路计算程序。
短路计算可以分为对称故障和不对称故障两种情况。
对称故障是指短路故障发生在电力系统的正常运行条件下,例如三相短路。
不对称故障是指短路故障发生在电力系统的不正常运行条件下,例如单相接地短路。
对于对称故障,我们可以使用节点电流法进行计算。
首先,应用基尔霍夫电流定律,根据电压和阻抗计算电流。
然后,根据节点电流方程和电流方程计算电流分布。
最后,根据电流分布计算短路电流和故障点的电压。
对于不对称故障,我们可以使用仿真方法进行计算。
首先,需要设置故障位置和故障类型,例如A相到地短路。
然后,根据故障位置和类型修改节点和边的参数,例如将故障位置的阻抗设置为零。
最后,使用数值方法求解电力系统的动态响应,得到短路电流和故障点的电压。
在MATLAB中,可以使用矩阵运算和数值求解函数实现短路计算。
例如,可以使用矩阵乘法和矩阵求逆函数计算节点电流和电流分布。
可以使用ODE求解器求解动态响应方程。
可以使用MATLAB的绘图函数绘制电力系统的电流分布和故障点的电压。
总结起来,电力系统短路故障分析的MATLAB辅助程序设计涉及建立电力系统模型、编写短路计算程序并使用MATLAB的数值计算和数据可视化工具进行计算和分析。
matlab基于节点阻抗矩阵的三相短路计算
![matlab基于节点阻抗矩阵的三相短路计算](https://img.taocdn.com/s3/m/02364c2d11a6f524ccbff121dd36a32d7275c740.png)
matlab基于节点阻抗矩阵的三相短路计算MATLAB基于节点阻抗矩阵的三相短路计算三相短路是电力系统中最常见的故障类型之一,也是最严重的一种故障,其产生的电流会对设备造成故障、损坏电力设备,甚至会导致火灾等事故。
因此,对电力系统进行三相短路计算及分析非常必要,MATLAB是一款通用的工具软件,可用于电力系统的短路计算中,本文就基于节点阻抗矩阵介绍MATLAB的三相短路计算。
一、节点阻抗矩阵节点阻抗矩阵是一种直观、简单、易于理解的方法,其基本思想是将电力系统中每个节点的短路电流计算单独列成一个向量,向量中每个元素都代表着该节点与其他节点之间的电流响应系数。
节点阻抗矩阵根据电力系统的拓扑结构所形成,元素值来自于两个节点之间的阻抗和电导之和。
二、节点阻抗矩阵的计算1、定义节点位置和电力设备参数在MATLAB中,首先需要定义电力系统中所有节点的位置参数(x,y)、所有支路的编号、阻抗参数值等,定义方法如下:node_position=[1 1; 1 1.5; 1.5 1; 1.5 1.5]; %四个节点的位置line=[1 2 0.02 0.04; 1 3 0.01 0.03; 2 4 0.03 0.06; 3 40.02 0.04];%四条线路的起始节点、结束节点、电阻、电抗2、构建节点阻抗矩阵根据节点位置和电力设备参数,可以通过以下语句构建节点阻抗矩阵:[~,Y_node]=Admittance_Matrix(line,4); %采用自定义函数计算导纳矩阵Z_node=inv(Y_node);%计算节点阻抗矩阵其中,Admittance_Matrix函数是一个自定义函数,用于求取系统的导纳矩阵,在导纳矩阵中,对称元素等于发电机到负载电工的导纳,非对称元素等于接线点的线路阻抗之和,函数的具体实现方法可以查阅MATLAB帮助文档。
三、三相短路计算有了节点阻抗矩阵,就可以进行三相短路计算,MATLAB中可以通过以下步骤进行计算:1、定义故障考虑的节点编号和负荷首先需要确定故障考虑的节点,即需要计算的节点,有多少个节点,就需要计算多少次,在节点矩阵中的位置,判定方法如下:fault_node=2;%节点2作为故障节点bus=[1 0; 2 10+4*j; 3 10+j; 4 15];%分别代表总发电机、故障节点、非故障节点1和非故障节点2的编号和电阻2、计算故障前的稳态电压在进行三相短路计算之前,需要先计算故障前的稳态电压值,详细计算方式可以参考MATLAB帮助文档或其他相关资料,计算发电机电动势和短路电流,具体计算方法如下:[num_bus,~]=size(bus);V_node=zeros(num_bus,num_bus);for i=1:num_busfor j=1:num_busV_node(i,j)=(bus(i,2)-bus(j,2))/Z_node(i,j);endend%计算电动势和短路电流Es=bus(1,2)-V_node(1,2)*Z_node(1,2);Fault_Pre_Curr=(Es-bus(fault_node,2))/Z_node(1,fault_node);3、计算短路电流和故障后电流得到故障前稳态电压后,可以根据一定的公式计算短路电流和故障后电流,具体计算方法如下:Z_fault=0.03+0.02*j;%故障阻抗I_fault=bus(fault_node,2)/(Z_node(fault_node,fault_node) +Z_fault);I_fault_phase=I_fault/(3^0.5);4、计算故障后电压最后,可以根据故障后电流计算出故障时间的电压,公式如下:V_fault_node=bus(fault_node,2)-I_fault*Z_node(fault_node,fault_node);总结本文简要介绍了MATLAB基于节点阻抗矩阵的三相短路计算方法,通过以上的步骤,可以较为准确地计算出电力系统中三相短路的电流和电压,为电力系统的安全稳定运行提供了重要保障。
三相短路和单相接地短路MATLAB
![三相短路和单相接地短路MATLAB](https://img.taocdn.com/s3/m/11d17785f705cc175527097f.png)
电力系统三相短路和单相接地短路
实验目的:
1.理解掌握短路的类型及在短路故障后的影响。
2.运用 MATLAB 的电力系统工具箱对三相短路和单相接地短路进行建模,
并分别观察分析其电压、电流的波形,并得出结论。
实验内容:
在二相电力系统中,大多数故障都是由于短路故障引起的,在发生短路故障的情况下,电力系统从一种状态剧烈变化到另一种状态,产生复杂的暂态现象。
在三相系统中,可能发生的短路有:三相短路、两相短路、一相短路接地和单相接地短路。
建立如图 3.1 理想情况下小型电力系统的模型。
图 3.1
参数设置:
电源模块图 3.2:
图 3.2三相电压电流测量模块图 3.3
图 3.3输电线路图 3.4:
图 3.4输电线路模块图 3.5:
图 3.5
万用表模块图 3.6:
(电压)
(电流)
图 3.6万用表 1 图 3.7
图 3.7三相序分量分析图 3.8:
图 3.8算法模块图 3.9:
图 3.9
单项接地:
故障模块 3.10:
图 3.10
对于测量模块和三相模块相同。
基于MATLAB的三相短路仿真与分析(2020年10月整理).pdf
![基于MATLAB的三相短路仿真与分析(2020年10月整理).pdf](https://img.taocdn.com/s3/m/f8d171b77cd184254a35357a.png)
研究生课程论文
(2014—2015 学年第 1 学期)
课程名称: 电力系统运行与控制
课程类型:
专业课
授课教师:
学 时: 36 学 分: 2.0
论文得分 批阅人签字 批阅意见:
论文题目 基于 MATLAB 的三相短路仿真与分析 姓名:
学号:
年级:
专业:
电气工程
学院:
电气学院
注意事项:
1、 以上各项由研究生本人认真填写; 2、 研究生课程论文应符合一般学术规范,具
1
学海无涯
目录
绪论 ..............................................................................................................1 第一章 供电系统仿真模型建立 .................................................................2
在三相系统中,可能发生的短路有:三相短路、两相短路、两相短路接地和单 相接地短路。三相短路也称为对称短路,系统各相与正常运行时一样仍处于对称状 态。其他类型的短路都不是对称短路。
电力系统的运行经验表明,在各类类型的短路中,单相短路占大多数,两相短 路较少,三相短路的机会最少。三相短路虽然很少发生,但情况较严重,应给以足 够的重视。况且,从短路计算方法来看,一切不对称短路的计算,在采用对称分量 法后,都归结为对称短路的计算。因此,对三相短路的研究是有其重要意义的。
kim −1 2
=
Im 2
1+ 2(kim −1)2
短路电流的最大有效值主要用于检验开关电器等设备切断短路电流的能力。
matlab实验 电力系统短路分析
![matlab实验 电力系统短路分析](https://img.taocdn.com/s3/m/e920f525ed630b1c59eeb5d1.png)
实验二 短路电流计算程序的实现一、三相短路电流计算程序计算短路电流周期分量,如I ''(I ')时,实际上就是求解交流电路的稳态电流,其数学模型也就是网络的线性代数方程,一般选用节点电压方程。
方程的系数矩阵是对称的。
在短路电流计算中变化的量往往是方程的常数项,需要多次求解线性方程组。
1.等值网络图2-1给出了不计负荷情况下计算短路电流I ''的等值网络。
在图2-1(a )中G 代表发电机端电压节点,发电机等值电势和电抗分别为E '' 和dx '',D 表示负荷节点,f 点为直接短路点。
应用叠加原理如图2-1所示。
正常运行方式为空载运行,网络中各点电压均为1;在故障分量网络中。
只需作故障分量的计算。
由图2-1的故障分量网络可见,这个网络与潮流计算的网络的差别在于发电机节点上多接了对地电抗dx ''。
当然如果短路计算中可以忽略线路电阻和电纳,而且不计变压器的实际变比,则短路计算网络较潮流计算网络简化,而且网络本身是纯感性的。
1E ''x1E '' x 1E '' x 1-=1=图2-1 在不计负荷情况下计算短路电流I ″的等值电路2. 用节电阻抗矩阵计算短路电流如果已经形成了故障分量网络的节点阻抗矩阵,则矩阵中的对角元素就是网络从f 点看进去的等值阻抗,又称为f 点的自阻抗。
fi Z 为f 点与i 点的互阻抗,均用大写Z 表示。
由节点方程中的第f 个方程:n fn f ff f f I Z I Z I Z U ++++=11。
ff Z 为其它节电电流为零时,节点f 的电压和电流之比,即网络对f 点的等值阻抗。
根据故障分量网络,直接应用戴维南定理可求得直接短路电流(由故障点流出)为fff ffz Z U I +=0(2-1)式中,f z 为接地阻抗;0f U 为f 点短路前的电压。
matlab电力系统快速解耦法潮流计算及短路计算程序
![matlab电力系统快速解耦法潮流计算及短路计算程序](https://img.taocdn.com/s3/m/3380bf77f56527d3240c844769eae009581ba23e.png)
电力系统快速解耦法潮流分析及短路计算一.程序设计的基本思想:(1)由于电力系统潮流分析中要利用到矩阵运算,复数运算,故采用matlab编程。
采用文件输入,将系统的各个参数以文件的形式输入,便于程序的通用化。
(2)本程序共有两个输入文件,分别为线路参数的文件,和已知的节点状态文件(PQ)(3)为了使程序不仅仅局限于计算9节点网络,在形成节点导纳的函数Yn()中,利用循环,找出线路首节点中的最大编号数,自动确定节点导纳矩阵的维数。
故对于任意n节点网络,均可以计算出节点导纳矩阵(4)在(3)的前提下,为了使程序支持系统增加节点,增加负荷等造成的PQ参数改变,或者PQ表的加长。
对程序做了如下优化。
首先,程序执行的基础是PQ表中平衡节点在第一行,接下来是PV节点,最后是PQ节点,如果系统添加节点,或者删除节点,均在PQ表的末端操作,会造成PQ表的顺序不是平衡节点、PV节点、PQ节点的顺序。
故引入了seqencing()函数,其作用就是不论输入的PQ表是什么顺序,在程序读入后均按平衡-PV-PQ的顺序排列。
其次,顺序打乱的PQ表必须与支路参数表对应,故在Yn()函数中加入了两段循环体,使之对应(见相应函数体注释)(5)在满足了上述4个条件后,程序便可以通用化了。
当然,由于水平有限,且程序未能由大量数据测试,故缺陷在所难免,这里仅是做了通用化的尝试。
在本文最末附加了该程序通用化的实例。
二、潮流计算框图三.定义相应的函数1.形成节点导纳矩阵的函数Yn()function Y=Yn(x,y)%定义一名为Yn的函数,其功能是自动识别输入表中节点的个数,形成相应的节点导纳矩阵[fid,message]=fopen(x,'r') ; %从x文件中读入支路参数if fid==-1; %判断文件是否正确打开error(message);end;[HeadPoint,HeadNumber, EndPoint,EndNumber,R,X,B,k]=textread(x,'%s %d %s %d %f %f %f %f'); %将读入的参数处理为以列为向量的数组fclose(fid);%关闭文件L=length(HeadNumber); %确定输入表的行数[Pointstyle,PointNumber,Ps,Qs,Uk,Ok]=seqencing(y);%调用seqencing函数,引入y文件中的PQ参数A=PointNumber;for i=1:L; %通过以下两循环体,实现PQ参数与支路参数的编号对应for j=1:L;if HeadNumber(i)==j;HeadNumber(i)=A(j);break;end;end;end;for i=1:L;for j=1:L;if EndNumber(i)==j;EndNumber(i)=A(j);break;end;end;end;Y=zeros(L,L); %根据txt文件中数据表的长度建立空的节点导纳矩阵for i=1:Lm=HeadNumber(i);n=EndNumber(i);if k(i)==0; %判断是否何种元件,为输电线元件if n~=0;Y(m,m)=Y(m,m)+1j*B(i)+1/(R(i)+1j*X(i));Y(n,n)=Y(n,n)+1j*B(i)+1/(R(i)+1j*X(i));Y(m,n)=Y(m,n)-1/(R(i)+1j*X(i));Y(n,m)=Y(n,m)-1/(R(i)+1j*X(i));elseY(m,m)=Y(m,m)+R(i)+1j*X(i);end;else %为变压器元件if n~=0;Y(m,m)=Y(m,m)+1/(R(i)+1j*X(i));Y(m,n)=Y(m,n)-1/(k(i)*(R(i)+1j*X(i)));Y(n,n)=Y(n,n)+1/(k(i)*k(i)*(R(i)+1j*X(i)));Y(n,m)=Y(n,m)-1/(k(i)*(R(i)+1j*X(i)));elseY(m,m)=Y(m,m)+R(i)+1j*X(i);end;end;end;maxm=HeadNumber(1);%通过下面两个循环体,确定输入表中节点编号的最大值,及为节点导纳矩阵的维数for i=1:L;if maxm<=HeadNumber(i);maxm=HeadNumber(i);end;end;maxn=EndNumber(1);for i=1:L;if maxn<=EndNumber(i);maxn=EndNumber(i);end;end;Y=Y(1:max(maxm,maxn),1:max(maxm,maxn));%形成导纳矩阵2.对不满足要求的PQ参数表进行排序的函数seqencing()function [Pointstyle,PointNumber,Ps,Qs,Uk,Ok]=seqencing(y)%定义名为seqencing的函数,其功能是在系统添加节点,或输入的PQ参数的顺序不满足要求时,对PQ参数表进行重新排序,保证平衡节点放在第一行,接下来是PV节点,最后是PQ节点[fid,message]=fopen(y,'r'); %从y文件中读入PQ参数if fid==-1; %判断文件是否正确打开error(message);end;[Pointstyle,PointNumber,Ps,Qs,Uk,Ok]=textread(y,'%f %f %f %f %f %f');fclose(fid);L=length(PointNumber);%通过以下两个循环体,完成对PQ输入表的重新排序,其思想是,在PQ参数之前加入一列Pointstyle用于标识节点类型,平衡节点为0,PV节点为1,PQ节点为2,以Pointstyle列为基准进行排序for i=1:L;for j=1:L-i;if Pointstyle(j)>Pointstyle(j+1);t=Pointstyle(j+1);Pointstyle(j+1)=Pointstyle(j);Pointstyle(j)=t;t=PointNumber(j+1);PointNumber(j+1)=PointNumber(j);PointNumber(j)=t;t=Ps(j+1);Ps(j+1)=Ps(j);Ps(j)=t;t=Qs(j+1);Qs(j+1)=Qs(j);Qs(j)=t;t=Uk(j+1);Uk(j+1)=Uk(j);Uk(j)=t;end;end;end;3、形成解耦算法B’矩阵的函数 formB1()function B1=formB1(x,y)%定义名为B1的函数形成解耦算法中的B’矩阵,得到的B’矩阵用B1表示[fid,message]=fopen(x,'r') ; %从x文件中读入支路参数if fid==-1; %判断文件是否正确打开error(message);end;[HeadPoint,HeadNumber,EndPoint,EndNumber,R,X,B,k]=textread(x,'%s %d %s %d %f %f %f %f'); %将读入的参数处理为以列为向量的数组fclose(fid);L=length(HeadNumber);[Pointstyle,PointNumber,Ps,Qs,Uk,Ok]=seqencing(y);%调用seqencing函数,引入y文件中的PQ参数A=PointNumber;%通过以下两循环体,实现PQ参数与支路参数的编号对应for i=1:L;for j=1:L;if HeadNumber(i)==j;HeadNumber(i)=A(j);break;end;end;end;for i=1:L;for j=1:L;if EndNumber(i)==j;EndNumber(i)=A(j);break;end;end;end;B1=zeros(L,L);for i=1:L %以行为单位,通过循环,用支路参数对B1进行修改,形成B’矩阵 m=HeadNumber(i);n=EndNumber(i);B1(m,m)=B1(m,m)-1/X(i);B1(n,n)=B1(n,n)-1/X(i);B1(m,n)=B1(m,n)+1/X(i);B1(n,m)=B1(n,m)+1/X(i);endmaxm=HeadNumber(1);for i=1:L;if maxm<=HeadNumber(i);maxm=HeadNumber(i);end;end;maxn=EndNumber(1);for i=1:L;if maxn<=EndNumber(i);maxn=EndNumber(i);end;end;B1=B1(2:max(maxm,maxn),2:max(maxm,maxn)); %形成B’矩阵4、形成解耦算法B’’矩阵的函数 formB11()function B11=formB11(x,y)%定义名为B11的函数形成解耦算法中B’'矩阵,用B11表示从x文件中读入支路参数确定Y,从y文件中读入PQ参数确定B11的维数,即除去平衡节点和pv节点,此处要求PQ参数录入时,将平衡节点和PQ节点放在前排,这一要求在Yn函数中通过seqencing函数已经满足Y=Yn(x,y);B=imag(Y);[Pointstyle,PointNumber,Ps,Qs,Uk,Ok]=seqencing(y);i=1;j=1;while Pointstyle(i)<=1;i=i+1;j=j+1;end;B11=B(j:end,j:end); %形成B’’矩阵5、计算正常情况下系统节点电压的函数 powerflow()function [U0,O0]=powerflow(x,y)%定义名为powerflow的函数,利用快速解耦算法来计算正常情况下系统内各个节点的电压和相角[Pointstyle,PointNumber,Ps,Qs,Uk,Ok]=seqencing(y);% 调用seqencing函数对PQ参数表进行排序Y=Yn(x,y); %形成节点导纳矩阵,Yn为n维B1=formB1(x,y); %形成解耦算法中的B矩阵,B1为n-1维B11=formB11(x,y); %形成解耦算法中的B'矩阵,B'为m维G=real(Y); %取Y的实部B=imag(Y); %取Y的虚部U0=Uk;O0=Ok;L=length(PointNumber);P=zeros(L,1);Q=zeros(L,1);dP=zeros(L,1);dQ=zeros(L,1);number=1;i=1;k=1;while Pointstyle(i)<=1; %通过k值确定系统中PQ节点的个数i=i+1;k=k+1;end;while number<100 %定义迭代次数上限为100次for i=2:L;sum1=0;for j=1:L;sum1=sum1+U0(j)*(G(i,j)*cos(O0(i)-O0(j))+B(i,j)*sin(O0(i)-O0(j))); %潮流方程,n-1维end;dP(i)=Ps(i)-U0(i)*sum1;endfor i=k:L;sum2=0;for j=1:L;sum2=sum2+U0(j)*(G(i,j)*sin(O0(i)-O0(j))-B(i,j)*cos(O0(i)-O0(j))); %潮流方程,m 维end;dQ(i)=Qs(i)-U0(i)*sum2;enddP1=dP(2:L)./U0(2:L);dQ1=dQ(k:L)./U0(k:L);a=max(norm(dP1,inf));b=max(norm(dQ1,inf));if max(a,b)<0.00001 %判断是否收敛break;disp(‘迭代’)disp(k);disp(‘次后收敛’);else %如不收敛,dO=-inv(B1)*dP1; %dO为n-1维dU=-inv(B11)*dQ1; %dU为m维zero1=zeros(k-1,1);zero2=[0];DU=[zero1;dU];DO=[zero2;dO];U0=U0+DU;O0=O0+DO;number=number+1;end;if number==100;disp('迭代100次后不收敛,迭代结束');end;end;四.对相应的系统进行潮流分析和短路计算定义完上述函数之后,可直接调用函数形成导纳矩阵,计算正常情况下的节点电压,进行短路计算计算短路电流,短路后各个节点电压以及支路潮流分布。
电力系统三相短路的编程计算
![电力系统三相短路的编程计算](https://img.taocdn.com/s3/m/1bc40426793e0912a21614791711cc7930b7787e.png)
电力系统三相短路的编程计算简介:电力系统三相短路是指电力系统发生故障时,三相电流短路连接,导致系统电流急剧增加,可能引发设备烧毁和电力系统瘫痪的危险。
因此,进行电力系统短路计算是非常重要的。
本文将介绍如何利用编程语言进行电力系统三相短路计算。
1.三相短路计算原理三相短路计算是指在给定电网拓扑、负荷数据以及故障点位置的情况下,通过计算得到短路电流的数值和相位。
计算过程主要包括以下几个步骤:1.确定故障位置:根据电网的拓扑结构确定故障点的位置。
2.分解复杂电网为简化电路:将复杂的电网拓扑结构分解为简化的电路模型。
3.计算电流分布:通过求解电流方程,得到各支路的电流分布。
4.计算短路电流:根据电流分布和故障点位置,计算短路电流大小和相位。
2.编程实现步骤2.1数据输入首先,需要从用户获取输入数据,包括电网拓扑结构、负荷数据和故障点位置。
可以通过编程语言提供的输入函数或者读取数据文件的方式获取数据。
2.2电网拓扑解析将输入的电网拓扑结构进行解析,构建电路模型。
可以采用图论算法进行解析,比如深度优先遍历或广度优先遍历,将电网拓扑结构转化为图的邻接矩阵或邻接表。
2.3电流分布计算利用电路分析的方法,根据电路模型和负荷数据,计算每个节点的电压值。
可以采用基尔霍夫电压法或者毕奥萨法进行计算。
2.4短路电流计算根据故障点位置和电流分布,计算各个支路的电流值和相位。
可以采用基尔霍夫电流法来计算短路电流。
2.5结果输出将计算得出的短路电流值和相位进行输出,可以利用编程语言提供的输出函数或者将结果保存到文件中。
3.编程计算实例以下是一个使用Python语言实现电力系统三相短路计算的示例代码:```pythonimport numpy as np#输入电网拓扑结构、负荷数据和故障点位置等数据topology = np.array([[0, 1, 1, 0],[1,0,1,1],[1,1,0,1],[0,1,1,0]])loads = np.array([100, 200, 150, 100])fault_location = 2#电压计算voltages = np.linalg.solve(topology, loads)#短路电流计算short_circuit_current = loads[fault_location] /np.sum(voltages)#输出结果print("短路电流大小:", short_circuit_current)```以上代码使用numpy库进行矩阵计算,首先输入电网拓扑结构、负荷数据和故障点位置等数据,然后利用numpy库提供的线性方程求解函数linalg.solve(计算电压值。
电力系统短路故障分析的MATLAB辅助程序设计短路计算程序
![电力系统短路故障分析的MATLAB辅助程序设计短路计算程序](https://img.taocdn.com/s3/m/f35f47a484868762caaed556.png)
电力系统短路故障分析的MATLAB辅助程序设计电力系统短路故障可分为三相对称短路故障(three-phase balanced faults)和不对称短路故障(unbalanced faults )。
不对称短路故障又分为单相接地短路故障(single line-to-ground fault)、两相短路故障 (line-to-line fault)以及两相接地短路故障(double line-to-ground fault)。
根据故障分析结果可以对继电保护装置、自动装置进行整定计算,我们可以建立算法来形成节点阻抗矩阵,利用节点阻抗矩阵来计算短路故障情况下的节点电压和线路电流。
一、三相对称短路故障进行三相短路计算需要两个程序:zbuild /zbuildpi和symfault程序,zbuild、zbuildpi用来在MATLAB中形成节点阻抗矩阵,symfault用来计算三相对称故障。
Zbus=zbuild(zdata)这里的参数zdata是一个(e×4)阶矩阵,e是拓扑图的总支路数目。
第一列和第二列为元素两端的节点编号,第三列和第四列分别是线路的电阻、电抗的标幺值。
连接在0节点和发电机节点之间的发电机阻抗可能是次暂态电抗、暂态电抗或同步电抗,而且这个矩阵中还包含并联电抗器和负荷阻抗。
Zbus=zbuildpi(linedata,gendata,yload)这个函数与潮流计算程序是相容的,第一个参数linedat a与潮流计算程序中的文件是一致的。
第一列和第二列为节点编号;第三列到第五列分别是线路的电阻、电抗以及1/2线路电纳值,这三项都为在统一基准容量下的标幺值;最后一列是变压器分接头位置,对线路来说,必须输入1;线路无输入顺序。
发电机参数不包含在Linedata参数中,而是包含在第二个参数gendata中,gendata是一个g×4阶矩阵,g是发电机总数。
第一列和第二列为0节点、发电机节点编号,第三列和第四列为发电机的暂态电阻和暂态电抗。
电力系统短路故障的Matlab算法
![电力系统短路故障的Matlab算法](https://img.taocdn.com/s3/m/3e846eff0242a8956bece4b5.png)
2011-2012学年度下学期电力系统分析课程设计电力系统短路故障的计算机算法程序设计专业:电气工程及其自动化班级:学号:学生姓名:指导教师:2012年 2 月 28 日信息工程学院课程设计任务书目录课程设计说明 (4)1 任务提出与方案论证 (5)1.1 任务提出 (5)1.2 方案论证 (6)2 设计思路 (7)2.1 短路电流的计算方法 (7)2.2 利用节点阻抗矩阵计算短路电流 (8)3 详细设计 (10)3.1分析与计算 (10)3.2程序主框图及主要数据变量说明 (13)4 总结 (18)参考文献 (19)课程设计说明本文根据电力系统三相对称短路的特点,建立了合理的三相短路的数学模型,在此基础上,形成电力系统短路电流实用汁算方法即节点阻抗矩阵的支路追加法。
编制了对任意一个电力系统在任意点发生短路故障时三相短路电流及其分布的通用计算机程序,该办法适用予各种复杂结构的电力系统,从一个侧面展示了计算机应用于电力系统的广阔前景。
电力系统的短路故障是严重的,而又是发生几率最多的故障,一般说来,最严重的短路是三相短路。
当发生短路时,其短路电流可达数万安以至十几万安,它们所产生的热效应和电动力效应将使电气设备遭受严重破环。
为此,当发生短路时,继电保护装置必须迅速切除故障线路,以避免故障部分继续遭受危害,并使非故障部分从不正常运行情况下解脱出来,这要求电气设备必须有足够的机械强度和热稳定度,开关电气设备必须具备足够的开断能力,即必须经得起可能最大短路的侵扰而不致损坏。
因此,电力系统短路电流计算是电力系统运行分析,设计计算的重要环节,许多电业设计单位和个人倾注极大精力从事这一工作的研究。
在电力设计中选择电气设备必须计算短路电流, 短路电流的计算是电气专业设计不可缺少的环节,是电力设计中最重要的计算之一。
传统的设计方法中, 短路电流的计算是以手工计算形式进行的 ,先通过手工方法化简电力网络, 求出各电源点对短路点的转移阻抗, 从而求出计算用电抗XJS , 再查找运算曲线, 以求得短路电流的周期分量。
发电机出口三相短路计算matlab
![发电机出口三相短路计算matlab](https://img.taocdn.com/s3/m/0465550fff4733687e21af45b307e87101f6f8ea.png)
发电机出口三相短路计算matlab发电机出口三相短路是指在电机出口端发生发电机绕组之间的相间短路故障。
这种故障可能会导致电机发生过电流和电磁力矩的异常增加,严重情况甚至会导致设备损坏。
因此,对发电机出口三相短路进行计算和分析非常重要。
在进行发电机出口三相短路计算时,可以借助Matlab工具进行相应的数值计算和仿真分析。
Matlab工具提供了一系列的电力系统仿真和计算函数,可以帮助我们进行发电机出口三相短路计算。
首先,我们可以使用Matlab中的复数运算功能来计算短路电流。
假设发电机出口处的电动势为E,短路电阻为Z,短路电流为I。
根据欧姆定律,我们可以得到I = E / Z。
在Matlab中,可以使用复数数学函数来计算电压和阻抗的复数表示,然后进行除法计算得到电流。
其次,我们需要考虑发电机出口三相短路的影响范围和电流传播路径。
可以使用Matlab中的电网等值方法来建立发电机的等值电路模型。
例如,可以使用Matlab中的电网等值工具箱来建立配电网的等值模型,并将发电机的等值电路与之相连。
然后,可以使用示波器模块来显示电流的传播路径和短路电流的分布情况。
此外,我们还可以借助Matlab的图形绘制功能来进行电流和电压波形的展示和分析。
通过绘制发电机出口端的电流和电压波形,可以清楚地观察到短路故障对电机运行的影响。
可以使用Matlab中的信号处理工具箱来进行频谱分析和谐波分析,以计算短路电流的谐波含量和幅值。
此外,我们还可以使用Matlab中的仿真工具来模拟发电机出口三相短路故障,并对故障后的电流和电压进行仿真分析。
例如,可以使用Matlab中的Simulink工具箱来建立电力系统的仿真模型,将发电机等值电路和负载模型相连,并对短路故障进行仿真。
可以通过仿真结果来判断发电机短路故障对系统稳定性和设备保护的影响。
综上所述,使用Matlab进行发电机出口三相短路计算非常方便。
通过Matlab提供的复数运算功能、电网等值方法、图形绘制功能和仿真工具,可以进行短路电流的计算、电流传播路径的分析、波形展示和仿真分析。
Matlab-Three-PhaseSource(三相电源)
![Matlab-Three-PhaseSource(三相电源)](https://img.taocdn.com/s3/m/481273c1250c844769eae009581b6bd97f19bcdb.png)
Matlab-Three-PhaseSource(三相电源)内部R-L 感性阻抗的三相电源的实现库电源模块库电源描述三相电源模块执⾏带有内部R L 阻抗的平衡的三相电压来源。
三电压来源在中与中性点连接成Y 或者其他⽅式。
你可以通过规定电源短路容量或者是直接规定电源内部的电阻和电感R 和L 或者是间接地设置X 与R ⽐率。
对话框和参数参数表线电压有效值内部线电压单位:VA 相初相⾓指定内部电压的相⾓为A 相相⾓,⽤度数表⽰。
三相电压的相序为正序。
这样,B 相和C 相的电压分别落后A 相120°和240°。
频率三相电源的频率单位:Hz内部连接⽅式三相电源的内部连接。
改变电源连接,电源模块图标就跟着改变。
选择其中⼀个连接⽅式:通过短路电流容量指定阻抗选择它表明你想通过短路电流容量和X/R 的⽐值来指定阻抗值。
基准电压下三相短路电流容量基准电压下三相电感短路功率单位:VA ,⽤来计算内部电感值。
只有选中通过短路电流容量指定阻抗选项时,才能够使⽤该参数。
在给定基准相电压V base (单位:V )和电感三相短路功率Psc (单位:VA )时,内部电感L (单位:H )可由下式得到:212base sc V L P fπ=? 基准电压相基准电压,单位:V ,⽤来指定三相短路电流容量。
基准电压值通常就是名义电压源的⼤⼩。
只有选中通过短路电流容量指定阻抗选项时,才能够使⽤该参数。
X/R ⽐值在标称电源频率或者X/R ⽐值内部电源阻抗品质因数下的X/R ⽐值。
只有选中通过短路电流容量指定阻抗选项时,才能够使⽤该参数。
在给定内部电阻R (单位:Ω)和电感三相短路功率X (单位:Ω)时,X/R ⽐值可由下式得到:2(/)/X fL R X R X Rπ==只有不选中通过短路电流容量指定阻抗选项时,才能够使⽤该参数。
电源电阻单位:Ω电感只有不选中通过短路电流容量指定阻抗选项时,才能够使⽤该参数。
电感单位:H潮流计算表表中的这些参数是 Powergui模块潮流计算的⼯具,这些潮流计算参数只能给模型进⾏初始化,对块模型和仿真参数没有影响。
(完整版)短路计算matlab程序
![(完整版)短路计算matlab程序](https://img.taocdn.com/s3/m/2cd5d9e4f46527d3240ce0ce.png)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Power System Analysis' Project: Matlab Program %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%clear allclose allclca=exp(j*2*pi/3);A=[1 1 1 ;1 a^2 a ;1 a a^2];Sref=900e6;Uref=300e3;Iref=Sref/(sqrt(3)*Uref);% data from the different documents:Positive;Negative;Zero;% run power flow:opt = mpoption('ENFORCE_Q_LIMS', 1);results=runpf(case9gs,opt);Vp=[results.bus(:,8).*(cos(pi/180*results.bus(:,9))+j*sin(pi/180*result s.bus(:,9)))];clc;% Menu:Bucle=1;while (Bucle==1)% Menu:disp('_________________________________________________________________ _______________________')disp('%%%%%% Program to calculate the symetrical and unsymetrical short circuit faults %%%%%%');disp(' 1) Symetrical fault (three phase short circuit fault)');disp(' 2) Unsymetrical fault (A single phase to ground fault)');disp(' 0) For exit the program');obtion=input('Choose the option you want to be calculated by typing its number:','s');disp('_________________________________________________________________ _______________________')% Calculation of the faultswitch (obtion)%%case'1'disp('You Have choose: SYMETRICAL FAULT');disp(' ');disp('Please follow the menu');% Ask for the fault's bus number and the Rf value:disp(' ');disp('*******************************************************************************')disp('* DATA *')disp('*******************************************************************************')bus=input ('Please, number the bus where the three phase short circuitfault occur. Bus number:');%bus=bus+1; % Since node zero is also analyzed, the bus number Xwill have a X+1 indexRf=input('Please write the value of the arc resistance Rf per phase.Rf=');% Calculation:If=zeros(8,1);If(bus)=Vp(bus)/(Zpos(bus,bus)+Rf);Vs=Vp-Zpos*If;% Showing the Results:disp(' ');disp('*******************************************************************************')disp('* RESULTS *')disp('*******************************************************************************')FaultCurrent=sprintf('The short circuit currents is: %gA. withangle %g?, Iref*abs(If(bus)),180/pi*angle(If(bus)));disp(FaultCurrent);figure();compass(If(bus),'blue');hold oncompass(If(bus)*a,'red');compass(If(bus)*a^2,'green');text(-abs(If(bus)),-abs(If(bus)),'Short Circuit current');disp(' ')disp ('The Voltages at all the buses in the system are:')disp ('_____________________________________________')disp (' Bus VAngle[Degree]')disp ('_____________________________________________')VoltageTokke= sprintf(' Tokke %gV %g?,Uref*abs(Vs(1)),180/pi*angle(Vs(1)));disp(VoltageTokke)VoltageVinje= sprintf(' Vinje %gV %g?,Uref*abs(Vs(2)),180/pi*angle(Vs(2)));disp(VoltageVinje)VoltageSonga= sprintf(' Songa %gV %g?,Uref*abs(Vs(3)),180/pi*angle(Vs(3)));disp(VoltageSonga)VoltageVemork= sprintf(' Vemork %gV %g?,Uref*abs(Vs(4)),180/pi*angle(Vs(4)));disp(VoltageVemork)VoltageRjukan= sprintf(' Rjunkan %gV %g?,Uref*abs(Vs(5)),180/pi*angle(Vs(5)));disp(VoltageRjukan)VoltageFlesaker=sprintf(' Flesaker %gV %g?,Uref*abs(Vs(6)),180/pi*angle(Vs(6)));disp(VoltageFlesaker) VoltageTveiten= sprintf(' Tveiten %gV %g?,Uref*abs(Vs(7)),180/pi*angle(Vs(7)));disp(VoltageTveiten) VoltageRod= sprintf(' Rod %gV %g?,Uref*abs(Vs(8)),180/pi*angle(Vs(8)));disp(VoltageRod)disp ('_____________________________________________')figure()subplot(2,4,1)compass(Vs(1),'b');hold oncompass(Vs(1)*a,'r');compass(Vs(1)*a^2,'g');text(-abs(Vs(1)),-abs(Vs(1)),'Tokke`s Voltage vectors');subplot(2,4,2)compass(Vs(2));hold oncompass(Vs(2)*a,'r');compass(Vs(2)*a^2,'g');text(-abs(Vs(2)),-abs(Vs(2)),'Vinje`s Voltage vectors');subplot(2,4,3)compass(Vs(3));hold oncompass(Vs(3)*a,'r');compass(Vs(3)*a^2,'g');text(-abs(Vs(3)),-abs(Vs(3)),'Songa`s Voltage vectors');subplot(2,4,4)compass(Vs(4));hold oncompass(Vs(4)*a,'r');compass(Vs(4)*a^2,'g');text(-abs(Vs(4)),-abs(Vs(4)),'Vemork`s Voltage vectors');subplot(2,4,5)compass(Vs(5));hold oncompass(Vs(5)*a,'r');compass(Vs(5)*a^2,'g');text(-abs(Vs(5)),-abs(Vs(5)),'Rjukan`s Voltage vectors');subplot(2,4,6)compass(Vs(6));hold oncompass(Vs(6)*a,'r');compass(Vs(6)*a^2,'g');text(-abs(Vs(6)),-abs(Vs(6)),'Flesaker`s Voltage vectors');subplot(2,4,7)compass(Vs(7));hold oncompass(Vs(7)*a,'r');compass(Vs(7)*a^2,'g');text(-abs(Vs(7)),-abs(Vs(7)),'Tveiten`s Voltage vectors');subplot(2,4,8)compass(Vs(8));hold oncompass(Vs(8)*a,'r');compass(Vs(8)*a^2,'g');text(-abs(Vs(8)),-abs(Vs(8)),'Rod`s Voltage vectors');disp(' ')disp('*******************************************************************************')disp('* CURRENTS *')disp('*******************************************************************************')Q=input ('Do you want to know any of the currents through two connectedpoints? Y/N:','s');k=0;while Q=='Y'k=k+1;B1(k)=input('Write the first bus number:');B2(k)=input ('Write the second bus number:');if Ypos(B1(k),B2(k))==0disp(' ')disp('Warning: *************There is no connectionbetween*************')disp(' ')k=k-1;elseIbranch(k)=(Vs(B1(k))-Vs(B2(k)))*(-Ypos(B1(k),B2(k)));endQ=input ('Do you want to know any more of the currents? Y/N:','s');enddisp(' ')disp('*******************************************************************************')disp('* CURRENTS RESULTS *')disp('*******************************************************************************')disp('')if k==0disp('There are no currents results');endfor K=1:kCurrent=sprintf ('From bus no %g to bus no %g the current flowingis: %gA %g?,B1(K),B2(K),Iref*abs(Ibranch(K)),180/pi*angle(Ibranch( K)));disp (Current)figure();compass(Ibranch(K),'blue');hold oncompass(Ibranch(K)*a,'red');compass(Ibranch(K)*a^2,'green');text(-abs(Ibranch(K)),-abs(Ibranch(K)),'Current betweenselected nodes');end%%case'2'disp('You Have choose: UNSYMETRICAL FAULT');disp('The program only calcutes single line to ground faults, it is consider to happen in phase "a"')disp(' ')disp('Please follow the menu');% Ask for the fault's bus number and the Rf value:disp(' ');disp('*******************************************************************************')disp('* DATA *')disp('*******************************************************************************')bus=input ('Please, number the bus where the single phase to groundfault occur. Bus number:');%bus=bus+1; % Since node zero is also analyzed, the bus number Xwill have a X+1 indexRf=input('Please write the value of the arc resistance Rf per phase.Rf=');% Calculation of the currents:Ifphaseground=3*Vp(bus)/(Zpos(bus,bus)+Zneg(bus,bus)+Zzero(bus,bus)+3*Rf);%Fase aIfa=zeros(8,1);Ifa(bus)=Ifphaseground;Ifa1=Ifa/3;Ifa2=Ifa/3;Ifa0=Ifa/3;for k=1:8Ifault=A*[Ifa0(k);Ifa1(k);Ifa2(k)];IF(k,:)=Ifault';end% Voltages;%Fase a:V1=Vp-Zpos*Ifa1;V2=-Zneg*Ifa2;V0=-Zzero*Ifa0;for k=1:8U=A*[V0(k);V1(k);V2(k)]V(k,:)=transpose(U)end% Showing the Results:disp(' ');disp('*******************************************************************************')disp('* RESULTS *')disp('*******************************************************************************')disp(' ')FaultCurrent=sprintf('The short circuit currents is: %g A. withangle %g ?, Iref*abs(Ifa(bus)),180/pi*angle(Ifa(bus)));disp(FaultCurrent);figure();compass(IF(bus,1),'red');hold oncompass(IF(bus,2),'blue');compass(IF(bus,3),'green');text(-abs(IF(bus,1)),-abs(IF(bus,1)),'Short Circuit current');disp(' ')disp ('The Voltages at all the buses in the system are:')disp('__________________________________________________________________________________________________')disp(' Bus Va VbVc');disp('__________________________________________________________________________________________________')VoltageTokke= sprintf('Tokke %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(1,1)),180/pi*angle(V(1,1)),Uref/sqrt(3)*abs(V(1,2) ),180/pi*angle(V(1,2)),Uref/sqrt(3)*abs(V(1,3)),180/pi*angle(V(1,3)));d isp(VoltageTokke)VoltageVinje= sprintf('Vinje %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(1,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(2,2) ),180/pi*angle(V(2,2)),Uref/sqrt(3)*abs(V(2,3)),180/pi*angle(V(2,3)));d isp(VoltageVinje)VoltageSonga= sprintf('Songa %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(3,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(3,2) ),180/pi*angle(V(3,2)),Uref/sqrt(3)*abs(V(3,3)),180/pi*angle(V(3,3)));d isp(VoltageSonga)VoltageVemork= sprintf('Vemork %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(4,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(4,2) ),180/pi*angle(V(4,2)),Uref/sqrt(3)*abs(V(4,3)),180/pi*angle(V(4,3)));d isp(VoltageVemork)VoltageRjukan= sprintf('Rjukan %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(5,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(5,2) ),180/pi*angle(V(5,2)),Uref/sqrt(3)*abs(V(5,3)),180/pi*angle(V(5,3)));d isp(VoltageRjukan)VoltageFlesaker=sprintf('Flesaker %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(6,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(6,2) ),180/pi*angle(V(6,2)),Uref/sqrt(3)*abs(V(6,3)),180/pi*angle(V(6,3)));d isp(VoltageFlesaker)VoltageTveiten= sprintf('Tveiten %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(7,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(7,2) ),180/pi*angle(V(7,2)),Uref/sqrt(3)*abs(V(7,3)),180/pi*angle(V(7,3)));d isp(VoltageTveiten)VoltageRod= sprintf('Rod %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(8,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(8,2) ),180/pi*angle(V(8,2)),Uref/sqrt(3)*abs(V(8,3)),180/pi*angle(V(8,3)));d isp(VoltageRod)disp('_________________________________________________________________ _________________________________')figure()subplot(2,4,1)compass(V(1,2),'b');hold oncompass(V(1,1),'r');compass(V(1,3),'g');text(-abs(V(1,2)),-abs(V(1,2)),'Tokke`s Voltage vectors');subplot(2,4,2)compass(V(2,2));hold oncompass(V(2,1),'r');compass(V(2,3),'g');text(-abs(V(2,2)),-abs(V(2,2)),'Vinje`s Voltage vectors');subplot(2,4,3)compass(V(3,2));hold oncompass(V(3,1),'r');compass(V(4,3),'g');text(-abs(V(3,2)),-abs(V(3,2)),'Songa`s Voltage vectors');subplot(2,4,4)compass(V(4,2));hold oncompass(V(4,1),'r');compass(V(4,3),'g');text(-abs(V(4,2)),-abs(V(4,2)),'Vemork`s Voltage vectors');subplot(2,4,5)compass(V(5,2));hold oncompass(V(5,1),'r');compass(V(5,3),'g');text(-abs(V(5,2)),-abs(V(5,2)),'Rjukan`s Voltage vectors');subplot(2,4,6)compass(V(6,2));hold oncompass(V(6,1),'r');compass(V(6,3),'g');text(-abs(V(6,2)),-abs(V(6,2)),'Flesaker`s Voltage vectors');subplot(2,4,7)compass(V(7,2));hold oncompass(V(7,1),'r');compass(V(7,3),'g');text(-abs(V(7,2)),-abs(V(7,3)),'Tveiten`s Voltage vectors');subplot(2,4,8)compass(V(8,2));hold oncompass(V(8,1),'r');compass(V(8,3),'g');text(-abs(V(8,2)),-abs(V(8,2)),'Rod`s Voltage vectors');disp(' ')disp('*******************************************************************************')disp('* CURRENTS *')disp('*******************************************************************************')Q=input ('Do you want to know any of the currents through two connectedpoints? Y/N:','s');k=0;while Q=='Y'k=k+1;B1(k)=input('Write the first bus number:');B2(k)=input ('Write the second bus number:');if Ypos(B1(k),B2(k))==0disp(' ')disp('Warning: *************There is no connectionbetween*************')disp(' ')k=k-1;elseIbranch0=(V0(B1(k))-V0(B2(k)))*(-Yzero((B1(k)),(B2(k))))Ibranch1=(V1(B1(k))-V1(B2(k)))*(-Ypos((B1(k)),(B2(k))))Ibranch2=(V2(B1(k))-V2(B2(k)))*(-Yneg((B1(k)),(B2(k))))IbranchK=A*[Ibranch0;Ibranch1;Ibranch2]Ibranch(k,:)=transpose(IbranchK)endQ=input ('Do you want to know any more of the currents? Y/N:','s');enddisp(' ')disp('*******************************************************************************')disp('* CURRENTS RESULTS *')disp('*******************************************************************************')disp('')if k==0disp('There are no currents results');enddisp('__________________________________________________________________________________________________')disp(' Buses Ia Ib Ic');disp('__________________________________________________________________________________________________')for K=1:kCurrent= sprintf(' %g-%g %gA and %g? %gAand %g? %gA and %g?',B1(K),B2(K),Iref*abs(Ibranch(K,1)),180/pi*angle(Ibranch(K,1)),Iref*abs(Ibranch(K,2)),180/pi*angle(Ibranch(K,2)),Iref*abs(Ibranch(K,3)),180/pi*angle(Ibranch(K,3)));disp(Current)figure();compass(Ibranch(K,1),'red');hold oncompass(Ibranch(K,2),'blue');compass(Ibranch(K,3),'green');text(-abs(Ibranch(K,1)),-abs(Ibranch(K,1)),'Current betweenselected nodes');enddisp('__________________________________________________________________________________________________')case'0'Bucle=0;disp(' ')disp(' ')disp('********************************************************************************')disp('* This is the result of our project *');disp('* We are: jy Wang*');disp('* Kalle*')disp('* Aravinda*');disp('* Pablo*');disp('********************************************************************************')otherwisedisp('Warning: *********************Wrong choose (try again thepromgram)*********************');end;input('') % until you don't press "enter" it doesn't continuoend;。
(完整版)短路计算matlab程序
![(完整版)短路计算matlab程序](https://img.taocdn.com/s3/m/2cd5d9e4f46527d3240ce0ce.png)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Power System Analysis' Project: Matlab Program %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%clear allclose allclca=exp(j*2*pi/3);A=[1 1 1 ;1 a^2 a ;1 a a^2];Sref=900e6;Uref=300e3;Iref=Sref/(sqrt(3)*Uref);% data from the different documents:Positive;Negative;Zero;% run power flow:opt = mpoption('ENFORCE_Q_LIMS', 1);results=runpf(case9gs,opt);Vp=[results.bus(:,8).*(cos(pi/180*results.bus(:,9))+j*sin(pi/180*result s.bus(:,9)))];clc;% Menu:Bucle=1;while (Bucle==1)% Menu:disp('_________________________________________________________________ _______________________')disp('%%%%%% Program to calculate the symetrical and unsymetrical short circuit faults %%%%%%');disp(' 1) Symetrical fault (three phase short circuit fault)');disp(' 2) Unsymetrical fault (A single phase to ground fault)');disp(' 0) For exit the program');obtion=input('Choose the option you want to be calculated by typing its number:','s');disp('_________________________________________________________________ _______________________')% Calculation of the faultswitch (obtion)%%case'1'disp('You Have choose: SYMETRICAL FAULT');disp(' ');disp('Please follow the menu');% Ask for the fault's bus number and the Rf value:disp(' ');disp('*******************************************************************************')disp('* DATA *')disp('*******************************************************************************')bus=input ('Please, number the bus where the three phase short circuitfault occur. Bus number:');%bus=bus+1; % Since node zero is also analyzed, the bus number Xwill have a X+1 indexRf=input('Please write the value of the arc resistance Rf per phase.Rf=');% Calculation:If=zeros(8,1);If(bus)=Vp(bus)/(Zpos(bus,bus)+Rf);Vs=Vp-Zpos*If;% Showing the Results:disp(' ');disp('*******************************************************************************')disp('* RESULTS *')disp('*******************************************************************************')FaultCurrent=sprintf('The short circuit currents is: %gA. withangle %g?, Iref*abs(If(bus)),180/pi*angle(If(bus)));disp(FaultCurrent);figure();compass(If(bus),'blue');hold oncompass(If(bus)*a,'red');compass(If(bus)*a^2,'green');text(-abs(If(bus)),-abs(If(bus)),'Short Circuit current');disp(' ')disp ('The Voltages at all the buses in the system are:')disp ('_____________________________________________')disp (' Bus VAngle[Degree]')disp ('_____________________________________________')VoltageTokke= sprintf(' Tokke %gV %g?,Uref*abs(Vs(1)),180/pi*angle(Vs(1)));disp(VoltageTokke)VoltageVinje= sprintf(' Vinje %gV %g?,Uref*abs(Vs(2)),180/pi*angle(Vs(2)));disp(VoltageVinje)VoltageSonga= sprintf(' Songa %gV %g?,Uref*abs(Vs(3)),180/pi*angle(Vs(3)));disp(VoltageSonga)VoltageVemork= sprintf(' Vemork %gV %g?,Uref*abs(Vs(4)),180/pi*angle(Vs(4)));disp(VoltageVemork)VoltageRjukan= sprintf(' Rjunkan %gV %g?,Uref*abs(Vs(5)),180/pi*angle(Vs(5)));disp(VoltageRjukan)VoltageFlesaker=sprintf(' Flesaker %gV %g?,Uref*abs(Vs(6)),180/pi*angle(Vs(6)));disp(VoltageFlesaker) VoltageTveiten= sprintf(' Tveiten %gV %g?,Uref*abs(Vs(7)),180/pi*angle(Vs(7)));disp(VoltageTveiten) VoltageRod= sprintf(' Rod %gV %g?,Uref*abs(Vs(8)),180/pi*angle(Vs(8)));disp(VoltageRod)disp ('_____________________________________________')figure()subplot(2,4,1)compass(Vs(1),'b');hold oncompass(Vs(1)*a,'r');compass(Vs(1)*a^2,'g');text(-abs(Vs(1)),-abs(Vs(1)),'Tokke`s Voltage vectors');subplot(2,4,2)compass(Vs(2));hold oncompass(Vs(2)*a,'r');compass(Vs(2)*a^2,'g');text(-abs(Vs(2)),-abs(Vs(2)),'Vinje`s Voltage vectors');subplot(2,4,3)compass(Vs(3));hold oncompass(Vs(3)*a,'r');compass(Vs(3)*a^2,'g');text(-abs(Vs(3)),-abs(Vs(3)),'Songa`s Voltage vectors');subplot(2,4,4)compass(Vs(4));hold oncompass(Vs(4)*a,'r');compass(Vs(4)*a^2,'g');text(-abs(Vs(4)),-abs(Vs(4)),'Vemork`s Voltage vectors');subplot(2,4,5)compass(Vs(5));hold oncompass(Vs(5)*a,'r');compass(Vs(5)*a^2,'g');text(-abs(Vs(5)),-abs(Vs(5)),'Rjukan`s Voltage vectors');subplot(2,4,6)compass(Vs(6));hold oncompass(Vs(6)*a,'r');compass(Vs(6)*a^2,'g');text(-abs(Vs(6)),-abs(Vs(6)),'Flesaker`s Voltage vectors');subplot(2,4,7)compass(Vs(7));hold oncompass(Vs(7)*a,'r');compass(Vs(7)*a^2,'g');text(-abs(Vs(7)),-abs(Vs(7)),'Tveiten`s Voltage vectors');subplot(2,4,8)compass(Vs(8));hold oncompass(Vs(8)*a,'r');compass(Vs(8)*a^2,'g');text(-abs(Vs(8)),-abs(Vs(8)),'Rod`s Voltage vectors');disp(' ')disp('*******************************************************************************')disp('* CURRENTS *')disp('*******************************************************************************')Q=input ('Do you want to know any of the currents through two connectedpoints? Y/N:','s');k=0;while Q=='Y'k=k+1;B1(k)=input('Write the first bus number:');B2(k)=input ('Write the second bus number:');if Ypos(B1(k),B2(k))==0disp(' ')disp('Warning: *************There is no connectionbetween*************')disp(' ')k=k-1;elseIbranch(k)=(Vs(B1(k))-Vs(B2(k)))*(-Ypos(B1(k),B2(k)));endQ=input ('Do you want to know any more of the currents? Y/N:','s');enddisp(' ')disp('*******************************************************************************')disp('* CURRENTS RESULTS *')disp('*******************************************************************************')disp('')if k==0disp('There are no currents results');endfor K=1:kCurrent=sprintf ('From bus no %g to bus no %g the current flowingis: %gA %g?,B1(K),B2(K),Iref*abs(Ibranch(K)),180/pi*angle(Ibranch( K)));disp (Current)figure();compass(Ibranch(K),'blue');hold oncompass(Ibranch(K)*a,'red');compass(Ibranch(K)*a^2,'green');text(-abs(Ibranch(K)),-abs(Ibranch(K)),'Current betweenselected nodes');end%%case'2'disp('You Have choose: UNSYMETRICAL FAULT');disp('The program only calcutes single line to ground faults, it is consider to happen in phase "a"')disp(' ')disp('Please follow the menu');% Ask for the fault's bus number and the Rf value:disp(' ');disp('*******************************************************************************')disp('* DATA *')disp('*******************************************************************************')bus=input ('Please, number the bus where the single phase to groundfault occur. Bus number:');%bus=bus+1; % Since node zero is also analyzed, the bus number Xwill have a X+1 indexRf=input('Please write the value of the arc resistance Rf per phase.Rf=');% Calculation of the currents:Ifphaseground=3*Vp(bus)/(Zpos(bus,bus)+Zneg(bus,bus)+Zzero(bus,bus)+3*Rf);%Fase aIfa=zeros(8,1);Ifa(bus)=Ifphaseground;Ifa1=Ifa/3;Ifa2=Ifa/3;Ifa0=Ifa/3;for k=1:8Ifault=A*[Ifa0(k);Ifa1(k);Ifa2(k)];IF(k,:)=Ifault';end% Voltages;%Fase a:V1=Vp-Zpos*Ifa1;V2=-Zneg*Ifa2;V0=-Zzero*Ifa0;for k=1:8U=A*[V0(k);V1(k);V2(k)]V(k,:)=transpose(U)end% Showing the Results:disp(' ');disp('*******************************************************************************')disp('* RESULTS *')disp('*******************************************************************************')disp(' ')FaultCurrent=sprintf('The short circuit currents is: %g A. withangle %g ?, Iref*abs(Ifa(bus)),180/pi*angle(Ifa(bus)));disp(FaultCurrent);figure();compass(IF(bus,1),'red');hold oncompass(IF(bus,2),'blue');compass(IF(bus,3),'green');text(-abs(IF(bus,1)),-abs(IF(bus,1)),'Short Circuit current');disp(' ')disp ('The Voltages at all the buses in the system are:')disp('__________________________________________________________________________________________________')disp(' Bus Va VbVc');disp('__________________________________________________________________________________________________')VoltageTokke= sprintf('Tokke %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(1,1)),180/pi*angle(V(1,1)),Uref/sqrt(3)*abs(V(1,2) ),180/pi*angle(V(1,2)),Uref/sqrt(3)*abs(V(1,3)),180/pi*angle(V(1,3)));d isp(VoltageTokke)VoltageVinje= sprintf('Vinje %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(1,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(2,2) ),180/pi*angle(V(2,2)),Uref/sqrt(3)*abs(V(2,3)),180/pi*angle(V(2,3)));d isp(VoltageVinje)VoltageSonga= sprintf('Songa %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(3,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(3,2) ),180/pi*angle(V(3,2)),Uref/sqrt(3)*abs(V(3,3)),180/pi*angle(V(3,3)));d isp(VoltageSonga)VoltageVemork= sprintf('Vemork %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(4,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(4,2) ),180/pi*angle(V(4,2)),Uref/sqrt(3)*abs(V(4,3)),180/pi*angle(V(4,3)));d isp(VoltageVemork)VoltageRjukan= sprintf('Rjukan %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(5,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(5,2) ),180/pi*angle(V(5,2)),Uref/sqrt(3)*abs(V(5,3)),180/pi*angle(V(5,3)));d isp(VoltageRjukan)VoltageFlesaker=sprintf('Flesaker %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(6,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(6,2) ),180/pi*angle(V(6,2)),Uref/sqrt(3)*abs(V(6,3)),180/pi*angle(V(6,3)));d isp(VoltageFlesaker)VoltageTveiten= sprintf('Tveiten %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(7,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(7,2) ),180/pi*angle(V(7,2)),Uref/sqrt(3)*abs(V(7,3)),180/pi*angle(V(7,3)));d isp(VoltageTveiten)VoltageRod= sprintf('Rod %gV and %g? %gVand %g? %gV and %g?',Uref/sqrt(3)*abs(V(8,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(8,2) ),180/pi*angle(V(8,2)),Uref/sqrt(3)*abs(V(8,3)),180/pi*angle(V(8,3)));d isp(VoltageRod)disp('_________________________________________________________________ _________________________________')figure()subplot(2,4,1)compass(V(1,2),'b');hold oncompass(V(1,1),'r');compass(V(1,3),'g');text(-abs(V(1,2)),-abs(V(1,2)),'Tokke`s Voltage vectors');subplot(2,4,2)compass(V(2,2));hold oncompass(V(2,1),'r');compass(V(2,3),'g');text(-abs(V(2,2)),-abs(V(2,2)),'Vinje`s Voltage vectors');subplot(2,4,3)compass(V(3,2));hold oncompass(V(3,1),'r');compass(V(4,3),'g');text(-abs(V(3,2)),-abs(V(3,2)),'Songa`s Voltage vectors');subplot(2,4,4)compass(V(4,2));hold oncompass(V(4,1),'r');compass(V(4,3),'g');text(-abs(V(4,2)),-abs(V(4,2)),'Vemork`s Voltage vectors');subplot(2,4,5)compass(V(5,2));hold oncompass(V(5,1),'r');compass(V(5,3),'g');text(-abs(V(5,2)),-abs(V(5,2)),'Rjukan`s Voltage vectors');subplot(2,4,6)compass(V(6,2));hold oncompass(V(6,1),'r');compass(V(6,3),'g');text(-abs(V(6,2)),-abs(V(6,2)),'Flesaker`s Voltage vectors');subplot(2,4,7)compass(V(7,2));hold oncompass(V(7,1),'r');compass(V(7,3),'g');text(-abs(V(7,2)),-abs(V(7,3)),'Tveiten`s Voltage vectors');subplot(2,4,8)compass(V(8,2));hold oncompass(V(8,1),'r');compass(V(8,3),'g');text(-abs(V(8,2)),-abs(V(8,2)),'Rod`s Voltage vectors');disp(' ')disp('*******************************************************************************')disp('* CURRENTS *')disp('*******************************************************************************')Q=input ('Do you want to know any of the currents through two connectedpoints? Y/N:','s');k=0;while Q=='Y'k=k+1;B1(k)=input('Write the first bus number:');B2(k)=input ('Write the second bus number:');if Ypos(B1(k),B2(k))==0disp(' ')disp('Warning: *************There is no connectionbetween*************')disp(' ')k=k-1;elseIbranch0=(V0(B1(k))-V0(B2(k)))*(-Yzero((B1(k)),(B2(k))))Ibranch1=(V1(B1(k))-V1(B2(k)))*(-Ypos((B1(k)),(B2(k))))Ibranch2=(V2(B1(k))-V2(B2(k)))*(-Yneg((B1(k)),(B2(k))))IbranchK=A*[Ibranch0;Ibranch1;Ibranch2]Ibranch(k,:)=transpose(IbranchK)endQ=input ('Do you want to know any more of the currents? Y/N:','s');enddisp(' ')disp('*******************************************************************************')disp('* CURRENTS RESULTS *')disp('*******************************************************************************')disp('')if k==0disp('There are no currents results');enddisp('__________________________________________________________________________________________________')disp(' Buses Ia Ib Ic');disp('__________________________________________________________________________________________________')for K=1:kCurrent= sprintf(' %g-%g %gA and %g? %gAand %g? %gA and %g?',B1(K),B2(K),Iref*abs(Ibranch(K,1)),180/pi*angle(Ibranch(K,1)),Iref*abs(Ibranch(K,2)),180/pi*angle(Ibranch(K,2)),Iref*abs(Ibranch(K,3)),180/pi*angle(Ibranch(K,3)));disp(Current)figure();compass(Ibranch(K,1),'red');hold oncompass(Ibranch(K,2),'blue');compass(Ibranch(K,3),'green');text(-abs(Ibranch(K,1)),-abs(Ibranch(K,1)),'Current betweenselected nodes');enddisp('__________________________________________________________________________________________________')case'0'Bucle=0;disp(' ')disp(' ')disp('********************************************************************************')disp('* This is the result of our project *');disp('* We are: jy Wang*');disp('* Kalle*')disp('* Aravinda*');disp('* Pablo*');disp('********************************************************************************')otherwisedisp('Warning: *********************Wrong choose (try again thepromgram)*********************');end;input('') % until you don't press "enter" it doesn't continuoend;。
电力系统短路电流计算matlab程序
![电力系统短路电流计算matlab程序](https://img.taocdn.com/s3/m/81749eb1ed3a87c24028915f804d2b160b4e86ff.png)
电力系统短路电流计算matlab程序%电力系统极坐标下的牛顿-拉夫逊法潮流计算disp('电力系统极坐标下的牛顿-拉夫逊法潮流计算:');clearn=input('请输入结点数:n=');n1=input('请输入PV结点数:n1=');n2=input('请输入PQ结点数:n2=');isb=input('请输入平衡结点:isb=');pr=input('请输入精确度:pr=');K=input('请输入变比矩阵看:K=');C=input('请输入支路阻抗矩阵:C=');y=input('请输入支路导纳矩阵:y=');U=input('请输入结点电压矩阵:U=');S=input('请输入各结点的功率:S=');Z=zeros(1,n);N=zeros(n1+n2,n2);L=zeros(n2,n2);QT1=zeros( 1,n1+n2); for m=1:nfor R=1:nC(m,m)=C(m,m)+y(m,R);if K(m,R)~=0C(m,m)=C(m,m)+1/((K(m,R)*C(m,R))/(K(m,R)-1));C(R,R)=C(R,R)+1/((K(m,R)^2*C(m,R))/(1-K(m,R)));C(m,R)=C(m,R)*K(m,R);C(R,m)=C(m,R);endendendfor m=1:nfor R=1:nif m~=RZ(m)=Z(m)+1/C(m,R); endendendfor m=1:nfor R=1:nif m==RY(m,m)=C(m,m)+Z(m); elseY(m,R)=-1/C(m,R);endendenddisp('结点导纳矩阵:');disp(Y);disp('迭代中的雅克比矩阵:'); G=real(Y);B=imag(Y);O=angle(U);U1=abs(U);k=0;PR=1;P=real(S);Q=imag(S);while PR>prfor m=1:n2UD(m)=U1(m);endfor m=1:n1+n2for R=1:nPT(R)=U1(m)*U1(R)*(G(m,R)*cos(O(m)-O(R))+B(m,R)*sin(O(m)-O(R)));endPT1(m)=sum(PT);PP(m)=P(m)-PT1(m);PP1(k+1,m)=PP(m);endfor m=1:n2for R=1:nQT(R)=U1(m)*U1(R)*(G(m,R)*sin(O(m)-O(R))-B(m,R)*cos(O(m)-O(R)));endQT1(m)=sum(QT);QQ(m)=Q(m)-QT1(m);QQ1(k+1,m)=QQ(m);endPR1=max(abs(PP));PR2=max(abs(QQ));PR=max(PR1,PR2);for m=1:n1+n2for R=1:n1+n2if m==RH(m,m)=U1(m)^2*B(m,m)+QT1(m);elseH(m,R)=-U1(m)*U1(R)*(G(m,R)*sin(O(m)-O(R))-B(m,R)*cos(O(m)-O(R)));endendendfor m=1:n1+n2for R=1:n2if m==RN(m,m)=-U1(m)^2*G(m,m)-PT1(m);elseN(m,R)=-U1(m)*U1(R)*(G(m,R)*cos(O(m)-O(R))+B(m,R)*sin(O(m)-O(R)));endendendfor m=1:n2for R=1:n1+n2if m==RJ(m,m)=U1(m)^2*G(m,m)-PT1(m);elseJ(m,R)=U1(m)*U1(R)*(G(m,R)*cos(O(m)-O(R))+B(m,R)*sin(O(m)-O(R)));endendendfor m=1:n2for R=1:n2if m==RL(m,m)=U1(m)^2*B(m,m)-QT1(m);elseL(m,R)=-U1(m)*U1(R)*(G(m,R)*sin(O(m)-O(R))-B(m,R)*cos(O(m)-O(R)));endendendJJ=[H N;J L];PQ=[PP';QQ'];DA=-inv(JJ)*PQ;DA1=DA';for m=1:n1+n2OO(m)=DA1(m);endfor m=n:n1+n2+n2UU1(m-n1-n2)=DA1(m); endUD2=diag(UD);UU=UU1*UD2;for m=1:n1+n2O(m)=O(m)+OO(m); endfor m=1:n2U1(m)=U1(m)+UU(m); endfor m=1:n1+n2o(k+1,m)=180/pi*O(m); endfor m=1:n2u(k+1,m)=U1(m);endk=k+1;endfor m=1:nb(m)=U1(m)*cos(O(m)); c(m)=U1(m)*sin(O(m)); endfor R=1:nPH1(R)=U(isb)*conj(Y(isb,R))*conj(U(R));endPH=sum(PH1);for m=1:nfor R=1:nif m~=RC1(m,R)=1/C(m,R);elseC1(m,m)=C(m,m);endendendfor m=1:nfor R=1:nif (C(m,R)~=inf)&(m~=R)SS(m,R)=U1(m)^2*conj(C1(m,m))+U(m)*(conj(U(m))-conj(U(R)))*conj(C1(m,R));endendenddisp('迭代中的△P:');disp(PP1);disp('迭代中的△Q:');disp(QQ1);disp('迭代中相角:');disp(o);disp('迭代中电压的模:');disp(u);disp('平衡结点的功率:');disp(PH);disp('全部线路功率分布:');disp(SS);。
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
实用文档电力系统三相短路计算a main.mclear tim%打开文件[dfile,pathname]=uigetfile('*.m','Select Data File');if pathname == 0error(' you must select a valid data file')elselfile =length(dfile);eval(dfile(1:lfile-2));end%定义输出文件output_file=fopen('output.dat','w');%开始计时tic;%求解节点导纳矩阵,其中Ymatrix1是考虑了变比,且支路未近似的导纳矩阵;Ymatrix2是近似变比为1,但是支路未近似计算的节点导纳矩阵;Ymatrix3是近似变比为1,采取近似支路参数1的导纳矩阵;Ymatrix4是近似变比为1,采取近似支路参数2的导纳矩阵。
Y = Ymatrix2(bus,line);%对故障点进行导纳修正fixY = FixY(Y,bus,fault);%求注入电流Iinj = Inode(bus,calcSettings);U = fixY\Iinj;%得到故障支路与其他支路电流Bcurrent = Ibranch( line,U,fault,Y );%如果发生支路三相短路,那么对应该支路的电流修正为-999999-j999999Ib = ReviseBcurrent( fault,Bcurrent );%结束计时tim=toc;fprintf('\n程序运行结果');fprintf('\n计算完成,共用时%4.4fs,相关结果已保存在output.dat\n',tim);%输出结果fprintf_result(output_file, Ib);fprintf_result1(Ib);b FixY.mfunction fixY = FixY( Y,bus,fault )%对形成的导纳矩阵进行故障点的修正[nb,mb]=size(bus);[nf,mf]= size(fault);fixY = Y;%对发电机节点导纳修正for k=1:nbbusType=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%对节点短路和支路短路的导纳矩阵进行修正for k=1:nfnodeI=fault(k,1);nodeJ=fault(k,2);dis=fault(k,3);if (nodeI==0)fixY(nodeJ,nodeJ) = 999999+1i*999999;continue;endif (nodeJ==0)fixY(nodeI,nodeI) = 999999+1i*999999;continue;endif (dis==0)&&(nodeI*nodeJ~=0)fixY(nodeI,nodeI) = 999999+1i*999999;continue;endif (dis==1)&&(nodeI*nodeJ~=0)fixY(nodeJ,nodeJ) = 999999+1i*999999;continue;endif (dis~=1)&&(dis~=0)&&(nodeI*nodeJ~=0)fixY(nodeI,nodeI) = 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;endendendc fprintf_result.mfunction [ output_args ] = fprintf_result( output_file, Ib )%将得到的短路电流输入到输出文件中[n,m]=size(Ib);fprintf( output_file, ' No. No. vector of I value of I\n');for k=1:nI=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);endif(I02<0)I02=abs(I02);fprintf( output_file, '%3d %3d %10.6f-j%10.6f %10.6f',I,J,I01,I02,I1);endfprintf( output_file, '\n');endendd fprintf_result1.mfunction [ output_args ] = fprintf_result1( Ib )%UNTITLED ÇëÔÚ´Ë´¦ÊäÈ뺯Êý¸ÅÒª[n,m]=size(Ib);fprintf(' No. No. vector of I value of I\n');for k=1:nI=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);endif(I02<0)I02=abs(I02);fprintf('%3d %3d %10.6f-j%10.6f %10.6f',I,J,I01,I02,I1);endfprintf('\n');endende Ibranch.mfunction Bcurrent = Ibranch( line,U,fault,Y )%计算短路电流%记录短路故障参数,如短路节点,如为支路短路,记录距离节点的距离%此段计算采用的支路参数未近似,如果计算近似的时候需要修改[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)faultNode = 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:nli=line(k,1);j=line(k,2);Ui=U(i);if j~=0Uj=U(j);elseUj=0;endif line(k,2)==0Ym=line(k,5)+1i*line(k,6);Iij=Ui*Ym;Iij1=abs(Iij);endif line(k,2)~=0Zt=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);endBcurrent(k,1)=i;Bcurrent(k,2)=j;Bcurrent(k,3)=Iij;Bcurrent(k,4)=Iij1;end%如果为节点短路,修正短路点的电流大小if(faultNode~=0)Bcurrent(nl+1,1) = faultNode;Bcurrent(nl+1,2) = faultNode;Ifault = 0;branchCurrent=0;for k=1:nlI=line(k,1);J=line(k,2);if(I*J==0)continue;endbranchCurrent = (U(I)-U(J))/(line(k,3)+1i*line(k,4));if (I==faultNode)Ifault = Ifault - branchCurrent ;else if (J==faultNode)Ifault = Ifault + branchCurrent ;endendendBcurrent(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)*Y(faultI,faultJ)/(1-dis);Bcurrent(nl+1,4) = abs(Bcurrent(nl+1,3));endendf Inode.mfunction Iinj = Inode( bus,calcSettings )%计算节点注入电流[nb,mb]=size(bus);Iinj = zeros(nb,1);for k=1:nbbusType=bus(k,7);if(calcSettings(1)==1)v = 1;elsev = bus(k,2);end%对发电机节点电流进行修正if (busType==1)Iinj(bus(k,1),1) = Iinj(bus(k,1),1) + v/1i/bus(k,8);endendendg ReviseBcurrent.mfunction Ib = ReviseBcurrent( fault,Bcurrent )%如果发生支路短路,对原来的计算电流进行修正,使该支路短路电流输出为-999999-j999999 clear faultI faultJ dis[nt,mt]=size(Bcurrent);Ib=zeros(nt,mt);faultI=fault(1,1);faultJ=fault(1,2);dis=fault(1,3);for k=1:nt-1i=Bcurrent(k,1);j=Bcurrent(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;endif (faultI*faultJ~=0)&&(dis~=1)&&(dis~=0)&&(i==faultJ)&&(j==faultI) Ib(k,1)=i;Ib(k,2)=j;Ib(k,3)=-999999-1i*999999;Ib(k,4)=-999999;endIb(nt,:)=Bcurrent(nt,:);endh Ymatrix1.mfunction Y = Ymatrix1( bus,line )%考虑变压器,并且支路参数不近似的节点导纳矩阵[nb,mb]=size(bus);[nl,ml]=size(line);Y=zeros(nb,nb);for k=1:nlI=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);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);endif (K==0)&&(J==0)Y(I,I)=Y(I,I)+Ym;endif K>0Y(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);endif K<0Y(I,I)=Y(I,I)+Yt+Ym;Y(J,J)=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 )用matlab实现电力系统潮流计算晏鸣宇%考虑变压器变比近似为1,支路参数不等效[nb,mb]=size(bus);[nl,ml]=size(line);Y=zeros(nb,nb);for k=1:nlI=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~=0Y(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);endif J==0Y(I,I)=Y(I,I)+Ym;endendendj Ymatrix3.mfunction Y = Ymatrix3( bus,line )%考虑变压器变比为1,采用支路参数近似1[nb,mb]=size(bus);[nl,ml]=size(line);Y=zeros(nb,nb);for k=1:nlI=line(k,1);J=line(k,2);Zt=line(k,3)+1i*line(k,4);Yt=imag(1/Zt);Ym=imag(line(k,5)+1i*line(k,6));if J~=0Y(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);endif J==0Y(I,I)=Y(I,I)+Ym;endend用matlab实现电力系统潮流计算晏鸣宇endk Ymatrix4.mfunction Y = Ymatrix4( bus,line )%变压器变比近似为1,采用支路等效参数2[nb,mb]=size(bus);[nl,ml]=size(line);Y=zeros(nb,nb);for k=1:nlI=line(k,1);J=line(k,2);Zt=1i*line(k,4);Yt=1/Zt;Ym=1i*line(k,6);if J~=0Y(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);endif J==0Y(I,I)=Y(I,I)+Ym;endendend。