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

用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 ÇëÔÚ´Ë´¦ÊäÈ뺯Êý¸ÅÒª
[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 )
％考虑变压器变比近似为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
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
9