电力系统三相短路计算的MATLAB代码.docx
《电力系统三相短路计算的MATLAB代码.docx》由会员分享，可在线阅读，更多相关《电力系统三相短路计算的MATLAB代码.docx（15页珍藏版）》请在冰点文库上搜索。
电力系统三相短路计算的MATLAB代码
电力系统三相短路计算
amain.m
cleartim
%打开文件
[dfile,pathname]=uigetfile（'*.m','SelectDataFile'）;
ifpathname==0
error（'youmustselectavaliddatafile'）
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）;
bFixY.m
functionfixY=FixY（Y,bus,fault）
%对形成的导纳矩阵进行故障点的修正
[nb,mb]=size（bus）;
[nf,mf]=size（fault）;
fixY=Y;
%对发电机节点导纳修正
fork=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
%对节点短路和支路短路的导纳矩阵进行修正
fork=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
cfprintf_result.m
function[output_args]=fprintf_result（output_file,Ib）
%将得到的短路电流输入到输出文件中
[n,m]=size（Ib）;
fprintf（output_file,'No.No.vectorofIvalueofI\n'）;
fork=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
dfprintf_result1.m
function[output_args]=fprintf_result1（Ib）
%UNTITLEDÇëÔÚ´Ë´¦ÊäÈ뺯Êý¸ÅÒª
[n,m]=size（Ib）;
fprintf（'No.No.vectorofIvalueofI\n'）;
fork=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
eIbranch.m
functionBcurrent=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
%计算非故障支路的短路电流
fork=1:
nl
i=line（k,1）;
j=line（k,2）;
Ui=U（i）;
ifj~=0
Uj=U（j）;
else
Uj=0;
end
ifline（k,2）==0
Ym=line（k,5）+1i*line（k,6）;
Iij=Ui*Ym;
Iij1=abs（Iij）;
end
ifline（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;
fork=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;
elseif（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
fInode.m
functionIinj=Inode（bus,calcSettings）
%计算节点注入电流
[nb,mb]=size（bus）;
Iinj=zeros（nb,1）;
fork=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
gReviseBcurrent.m
functionIb=ReviseBcurrent（fault,Bcurrent）
%如果发生支路短路，对原来的计算电流进行修正，使该支路短路电流输出为-999999-j999999
clearfaultIfaultJdis
[nt,mt]=size（Bcurrent）;
Ib=zeros（nt,mt）;
faultI=fault（1,1）;
faultJ=fault（1,2）;
dis=fault（1,3）;
fork=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
hYmatrix1.m
functionY=Ymatrix1（bus,line）
%考虑变压器，并且支路参数不近似的节点导纳矩阵
[nb,mb]=size（bus）;
[nl,ml]=size（line）;
Y=zeros（nb,nb）;
fork=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
ifK>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
ifK<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
iYmatrix2.m
functionY=Ymatrix2（bus,line）
%考虑变压器变比近似为1，支路参数不等效
[nb,mb]=size（bus）;
[nl,ml]=size（line）;
Y=zeros（nb,nb）;
fork=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）;
ifJ~=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
ifJ==0
Y（I,I）=Y（I,I）+Ym;
end
end
end
jYmatrix3.m
functionY=Ymatrix3（bus,line）
%考虑变压器变比为1，采用支路参数近似1
[nb,mb]=size（bus）;
[nl,ml]=size（line）;
Y=zeros（nb,nb）;
fork=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））;
ifJ~=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
ifJ==0
Y（I,I）=Y（I,I）+Ym;
end
end
end
kYmatrix4.m
functionY=Ymatrix4（bus,line）
%变压器变比近似为1，采用支路等效参数2
[nb,mb]=size（bus）;
[nl,ml]=size（line）;
Y=zeros（nb,nb）;
fork=1:
nl
I=line（k,1）;
J=line（k,2）;
Zt=1i*line（k,4）;
Yt=1/Zt;
Ym=1i*line（k,6）;
ifJ~=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
ifJ==0
Y（I,I）=Y（I,I）+Ym;
end
end
end
升级会员 