llgfault.m

一个用matlab编写的用于故障诊断的源程序

function done=llgfault(nb,ng,nl,nt)
%formation of positive sequence Ybus
Yp=Ypositive(nb,ng,nl,nt);
%formation of zero sequence Ybus
Yz=Yzero(nb,ng,nl,nt);
%formation of positive and negative sequence Zbus
Zp=inv(Yp);
Zn=Zp;
%formation of zero sequence Zbus
check=diag(Yz);
k=find(check==0)
if length(k)==0
    Zz=inv(Yz);
end
if length(k)~=0
    Zz=pinv(Yz);
end
%preparation of pre-fault voltage vectors
xxx=ones(nb);
vp0=xxx(:,1);
yyy=zeros(nb);
vn0=yyy(:,1);
vz0=vn0;
%calculation of fault current
fid=fopen('faultdata.txt','r');
b=textread('faultdata.txt');
fclose(fid);
bus=b(1);
imp=b(3);
%calculation of Thevenin impedances
Zth=((Zz(bus,bus)+3*imp)*j*Zn(bus,bus)*j/((Zz(bus,bus)+3*imp+Zn(bus,bus))*j))+Zp(bus,bus)*j;
ifp=vp0(bus)/Zth;
ifn=-ifp*(Zz(bus,bus)+3*imp)*j/(Zz(bus,bus)+3*imp+Zn(bus,bus));
ifz=-ifp*(Zn(bus,bus))*j/(Zz(bus,bus)+3*imp+Zn(bus,bus));
%post fault voltage calculations
for i=1:nb
    vpf(i)=vp0(i)-Zp(i,bus)*ifp*j;
    vnf(i)=vn0(i)-Zn(i,bus)*ifn*j;
    vzf(i)=vz0(i)-Zz(i,bus)*ifz*j;
end
vpf_mag=abs(vpf);
vpf_phase=angle(vpf);
vnf_mag=abs(vnf);
vnf_phase=angle(vnf);
vzf_mag=abs(vzf);
vzf_phase=angle(vzf);

%calculation of post fault flows in generators, transformers and line
%reading of elements and classifying them as lines, transformers or
%generators
%the convention for the status of elements is as follows
% 1=lines
% 2=generators
% 3=transformers without phase shift
% 0=transformers with phase shift
fid=fopen('linedata.txt','r');
bb=textread('linedata.txt');
fclose(fid);
x=bb(:,1);
y=bb(:,2);
zp=bb(:,3);
zz=bb(:,4);
ele_num=1:length(x);
start_bus=x;
end_bus=y;
for i=1:length(x)
    status(i)=1;
end
for i=1:nl
    poslineflowp(i)=(vpf_mag(x(i))*(cos(vpf_phase(x(i)))+sin(vpf_phase(x(i)))*j)-vpf_mag(y(i))*(cos(vpf_phase(y(i)))+sin(vpf_phase(y(i)))*j))/(zp(i)*j);
    poslineflown(i)=(vnf_mag(x(i))*(cos(vnf_phase(x(i)))+sin(vnf_phase(x(i)))*j)-vnf_mag(y(i))*(cos(vnf_phase(y(i)))+sin(vnf_phase(y(i)))*j))/(zp(i)*j);
    poslineflowz(i)=(vzf_mag(x(i))*(cos(vzf_phase(x(i)))+sin(vzf_phase(x(i)))*j)-vzf_mag(y(i))*(cos(vzf_phase(y(i)))+sin(vzf_phase(y(i)))*j))/(zz(i)*j);
    neglineflowp(i)=(vpf_mag(y(i))*(cos(vpf_phase(y(i)))+sin(vpf_phase(y(i)))*j)-vpf_mag(x(i))*(cos(vpf_phase(x(i)))+sin(vpf_phase(x(i)))*j))/(zp(i)*j);
    neglineflown(i)=(vnf_mag(y(i))*(cos(vnf_phase(y(i)))+sin(vnf_phase(y(i)))*j)-vnf_mag(x(i))*(cos(vnf_phase(x(i)))+sin(vnf_phase(x(i)))*j))/(zp(i)*j);
    neglineflowz(i)=(vzf_mag(y(i))*(cos(vzf_phase(y(i)))+sin(vzf_phase(y(i)))*j)-vzf_mag(x(i))*(cos(vzf_phase(x(i)))+sin(vzf_phase(x(i)))*j))/(zz(i)*j);
 end
 poslineflowp_mag=abs(poslineflowp);
 poslineflowp_phase=angle(poslineflowp);
 poslineflown_mag=abs(poslineflown);
 poslineflown_phase=angle(poslineflown);
 poslineflowz_mag=abs(poslineflowz);
 poslineflowz_phase=angle(poslineflowz);
 neglineflowp_mag=abs(neglineflowp);
 neglineflowp_phase=angle(neglineflowp);
 neglineflown_mag=abs(neglineflown);
 neglineflown_phase=angle(neglineflown);
 neglineflowz_mag=abs(neglineflowz);
 neglineflowz_phase=angle(neglineflowz)

fid=fopen('transformerdata.txt','r');
bb=textread('transformerdata.txt');
fclose(fid);
x=bb(:,1);
y=bb(:,3);
z=bb(:,5);
zz=bb(:,6);
hv=bb(:,2);
lv=bb(:,4);
hvbuses=x;
ele_num=1:length([ele_num y']);
start_bus=[start_bus' x'];
end_bus=[end_bus' y'];
for i=1:length(x)
    if (hv(i)==1 & lv(i)==2) | (hv(i)==2 & lv(i)==1) | (hv(i)==0 & lv(i)==2) | (hv(i)==2 & lv(i)==0)
        status(i+nl)=0;
    else
        status(i+nl)=3;
    end
end
s=length(status);
for i=1:nt
    postransflowp(i)=(vpf_mag(x(i))*(cos(vpf_phase(x(i)))+sin(vpf_phase(x(i)))*j)-vpf_mag(y(i))*(cos(vpf_phase(y(i)))+sin(vpf_phase(y(i)))*j))/(z(i)*j);
    postransflown(i)=(vnf_mag(x(i))*(cos(vnf_phase(x(i)))+sin(vnf_phase(x(i)))*j)-vnf_mag(y(i))*(cos(vnf_phase(y(i)))+sin(vnf_phase(y(i)))*j))/(z(i)*j);
    postransflowz(i)=(vzf_mag(x(i))*(cos(vzf_phase(x(i)))+sin(vzf_phase(x(i)))*j)-vzf_mag(y(i))*(cos(vzf_phase(y(i)))+sin(vzf_phase(y(i)))*j))*(-j)*Yz(x(i),y(i));
    negtransflowp(i)=(vpf_mag(y(i))*(cos(vpf_phase(y(i)))+sin(vpf_phase(y(i)))*j)-vpf_mag(x(i))*(cos(vpf_phase(x(i)))+sin(vpf_phase(x(i)))*j))/(z(i)*j);
    negtransflown(i)=(vnf_mag(y(i))*(cos(vnf_phase(y(i)))+sin(vnf_phase(y(i)))*j)-vnf_mag(x(i))*(cos(vnf_phase(x(i)))+sin(vnf_phase(x(i)))*j))/(z(i)*j);
    negtransflowz(i)=(vzf_mag(y(i))*(cos(vzf_phase(y(i)))+sin(vzf_phase(y(i)))*j)-vzf_mag(x(i))*(cos(vzf_phase(x(i)))+sin(vzf_phase(x(i)))*j))*(-j)*Yz(x(i),y(i));
end
postransflowp_mag=abs(postransflowp);
postransflowp_phase=angle(postransflowp);
postransflown_mag=abs(postransflown);
postransflown_phase=angle(postransflown);
postransflowz_mag=abs(postransflowz);
postransflowz_phase=angle(postransflowz);
negtransflowp_mag=abs(negtransflowp);
negtransflowp_phase=angle(negtransflowp);
negtransflown_mag=abs(negtransflown);
negtransflown_phase=angle(negtransflown);
negtransflowz_mag=abs(negtransflowz);
negtransflowz_phase=angle(negtransflowz);

fid=fopen('gendata.txt','r');
bb=textread('gendata.txt');
fclose(fid);
x=bb(:,1);
z=bb(:,2);
zz=bb(:,3);
gnd=bb(:,4);
ele_num=1:length([ele_num y']);
%length(ele_num)
start_bus=[start_bus x'];
for i=1:length(x)
    temp(i)=nb+1; %nb+1 is taken as the code for the reference bus i.e ground
end
end_bus=[end_bus temp];
for i=1:ng
    status(s+i)=2;
end
for i=1:ng
    posgenflowp(i)=(vpf_mag(x(i))*(cos(vpf_phase(x(i)))+sin(vpf_phase(x(i)))*j)-1)/(z(i)*j);
    posgenflown(i)=(vnf_mag(x(i))*(cos(vnf_phase(x(i)))+sin(vnf_phase(x(i)))*j))/(z(i)*j);
    posgenflowz(i)=(vzf_mag(x(i))*(cos(vzf_phase(x(i)))+sin(vzf_phase(x(i)))*j))/((zz(i)+3*gnd(i))*j);
    neggenflowp(i)=-(vpf_mag(x(i))*(cos(vpf_phase(x(i)))+sin(vpf_phase(x(i)))*j)-1)/(z(i)*j);
    neggenflown(i)=-(vnf_mag(x(i))*(cos(vnf_phase(x(i)))+sin(vnf_phase(x(i)))*j))/(z(i)*j);
    neggenflowz(i)=-(vzf_mag(x(i))*(cos(vzf_phase(x(i)))+sin(vzf_phase(x(i)))*j))/((zz(i)+3*gnd(i))*j);
end
posgenflowp_mag=abs(posgenflowp);
posgenflowp_phase=angle(posgenflowp);
posgenflown_mag=abs(posgenflown);
posgenflown_phase=angle(posgenflown);
posgenflowz_mag=abs(posgenflowz);
posgenflowz_phase=angle(posgenflowz);
neggenflowz_mag=abs(neggenflowz);
neggenflowz_phase=angle(neggenflowz);
neggenflowp_mag=abs(neggenflowp);
neggenflowp_phase=angle(neggenflowp);
neggenflown_mag=abs(neggenflown);
neggenflown_phase=angle(neggenflown);
    
posflowp=[poslineflowp postransflowp posgenflowp];
posflown=[poslineflown postransflown posgenflown];
posflowz=[poslineflowz postransflowz posgenflowz];

posflowp_phase=[poslineflowp_phase postransflowp_phase posgenflowp_phase];
posflown_phase=[poslineflown_phase postransflown_phase posgenflown_phase];
posflowz_phase=[poslineflowz_phase postransflowz_phase posgenflowz_phase];

negflowp=[neglineflowp negtransflowp neggenflowp];
negflown=[neglineflown negtransflown neggenflown];
negflowz=[neglineflowz negtransflowz neggenflowz];

negflowp_phase=[neglineflowp_phase negtransflowp_phase neggenflowp_phase];
negflown_phase=[neglineflown_phase negtransflown_phase neggenflown_phase];
negflowz_phase=[neglineflowz_phase negtransflowz_phase neggenflowz_phase];

%correction for star-delta connected transformers
for i=1:(nb+1)
    M(i)=1000; %phase adjusting parameter of each bus
end
M(1)=0;
k=1; %running variable for element number
kk=find(M==1000);
while length(kk)~=0
    flag=k;
    ib=start_bus(k);
    jb=end_bus(k);
    
      
    if M(ib)==1000
        if M(jb)==1000
            k=k+1;
        end
        if M(jb)~=1000
            if status(k)~=0
                M(ib)=M(jb);
                k=k+1;
            end
            temp=find(hvbuses==ib);
            if k<=(length(ele_num)) & k==flag
            if status(k)==0
                if (length(temp)==0)
                    M(ib)=M(jb)-1;
                    k=k+1;
                else
                    M(ib)=M(jb)+1;
                    k=k+1;
                end