电力系统三相短路计算的MATLAB代码

电力系统三相短路计算

a main.m

clear 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;

%求解节点导纳矩阵，其中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-j999999

Ib = 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.m

function fixY = FixY( Y,bus,fault )

%对形成的导纳矩阵进行故障点的修正

[nb,mb]=size(bus);

[nf,mf]= size(fault);

fixY = Y;

%对发电机节点导纳修正

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);

end

end

%对节点短路和支路短路的导纳矩阵进行修正

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 (nodeJ==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,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;

end

end

end

c fprintf_result.m

function [ 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: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.6f',I,J,I01,I02,I1);

end

fprintf( output_file, '\n');

end

end

d fprintf_result1.m

function [ 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: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('\n');

end

end

e Ibranch.m

function 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;

end

if(faultJ==0)

faultNode = faultI;

end

if(dis==1)&&(faultI*faultJ~=0)

faultNode = faultJ;

end

if(dis==0)&&(faultI*faultJ~=0)

faultNode = faultI;

end

if(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);

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

%如果为节点短路，修正短路点的电流大小

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 = 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)*Y(faultI,faultJ)/(1-dis);

Bcurrent(nl+1,4) = abs(Bcurrent(nl+1,3));

end

end

f Inode.m

function Iinj = Inode( bus,calcSettings )

%计算节点注入电流

[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

%对发电机节点电流进行修正

if (busType==1)

Iinj(bus(k,1),1) = Iinj(bus(k,1),1) + v/1i/bus(k,8);

end

end

end

g ReviseBcurrent.m

function 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-1

i=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;

end

if (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;

end

Ib(nt,:)=Bcurrent(nt,:);

end

h Ymatrix1.m

function Y = Ymatrix1( bus,line )

%考虑变压器，并且支路参数不近似的节点导纳矩阵[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);

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)=Y(J,J)+K*K*Yt;

Y(I,J)=Y(I,J)+K*Yt;

Y(J,I)=Y(I,J);

end

end

end

i Ymatrix2.m

function Y = Ymatrix2( bus,line )

用matlab实现电力系统潮流计算晏鸣宇

%考虑变压器变比近似为1，支路参数不等效

[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,I)+Ym;

end

end

end

j Ymatrix3.m

function Y = Ymatrix3( bus,line )

%考虑变压器变比为1，采用支路参数近似1

[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=imag(1/Zt);

Ym=imag(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,I)+Ym;

end

end

用matlab实现电力系统潮流计算晏鸣宇end

k Ymatrix4.m

function Y = Ymatrix4( bus,line )

%变压器变比近似为1，采用支路等效参数2

[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=1i*line(k,4);

Yt=1/Zt;

Ym=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,I)+Ym;

end

end

end