llfault.m

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

function done=llfault(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);
ifp=vp0(bus)/((Zp(bus,bus)+Zn(bus,bus)+imp)*j);
ifn=-ifp;
ifz=0;
%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
            end
        end
        end
        kk=find(M==1000);
        if length(kk)==0
            %M
            break;
        end
        if length(kk)~=0
            if k<=length(ele_num)
                %k=k+1;
                %M
                continue;
            end
            if k>length(ele_num)
                k=1;
                %M
                continue;
            end
        end
    end
    
    if M(ib)~=1000
        if M(jb)~=1000
            k=k+1;
        end
        if M(jb)==1000
            if status(k)~=0
                M(jb)=M(ib);
                k=k+1;
            end
            temp=find(hvbuses==jb);
            if k<=(length(ele_num)) & k==flag
            if status(k)==0
                if (length(temp)==0)
                    M(jb)=M(ib)-1;
                    k=k+1;
                else
                    M(jb)=M(ib)+1;
                    k=k+1;
                end
            end
        end
        end
        kk=find(M==1000);
        if length(kk)==0
            %M
            break;
        end
        if length(kk)~=0
            if k<=length(ele_num)
                %k=k+1;
                %M
                continue;
            end
            if k>length(ele_num)
                k=1;
                %M
                continue;
            end
        end
    end
end

fid=fopen('faultdata.txt','r');
b=textread('faultdata.txt');
fclose(fid);
FB=b(1);
flag=M(FB);
for i=1:(nb+1)
    M(i)=M(i)-flag;
end
M;
%correction of voltage and flows
value=input('Enter the phase shift that you desire in degrees. The choice that you have are 30 and 90 degrees.');
for i=1:nb
    if value==30
        vpf_phase(i)=vpf_phase(i)+M(i)*(value*pi/180);
        vnf_phase(i)=vnf_phase(i)-M(i)*(value*pi/180);
    end
    if value==90
        vpf_phase(i)=vpf_phase(i)-M(i)*(value*pi/180);
        vnf_phase(i)=vnf_phase(i)+M(i)*(value*pi/180);
    end
end
for i=1:length(posflowp_phase)
    if value==30
        posflowp_phase(i)=posflowp_phase(i)+M(start_bus(i))*(value*pi/180);
        posflown_phase(i)=posflown_phase(i)-M(start_bus(i))*(value*pi/180);
    end
    if value==90
        posflowp_phase(i)=posflowp_phase(i)-M(start_bus(i))*(value*pi/180);
        posflown_phase(i)=posflown_phase(i)+M(start_bus(i))*(value*pi/180);
    end
    if value==30
        negflowp_phase(i)=negflowp_phase(i)+M(end_bus(i))*(value*pi/180);
        negflown_phase(i)=negflown_phase(i)-M(end_bus(i))*(value*pi/180);
    end
    if value==90
        negflowp_phase(i)=negflowp_phase(i)-M(end_bus(i))*(value*pi/180);
        negflown_phase(i)=negflown_phase(i)+M(end_bus(i))*(value*pi/180);
    end
end

%conversion to phase quantities
a=cos(2*pi/3)+sin(2*pi/3)*j;
asq=cos(4*pi/3)+sin(4*pi/3)*j;
for i=1:nb
    va(i)=vzf(i)+vpf(i)+vnf(i);
    vb(i)=vzf(i)+asq*vpf(i)+a*vnf(i);
    vc(i)=vzf(i)+a*vpf(i)+asq*vnf(i);
end
for i=1:length(ele_num)
    posflowa(i)=posflowz(i)+posflowp(i)+posflown(i);
    posflowb(i)=posflowz(i)+asq*posflowp(i)+a*posflown(i);
    posflowc(i)=posflowz(i)+a*posflowp(i)+asq*posflown(i);
end
for i=1:length(ele_num)
    negflowa(i)=negflowz(i)+negflowp(i)+negflown(i);
    negflowb(i)=negflowz(i)+asq*negflowp(i)+a*negflown(i);
    negflowc(i)=negflowz(i)+a*negflowp(i)+asq*negflown(i);
end

%data entry into output file
fid=fopen('output.txt','w');
fprintf(fid,'The positive sequence Ybus\n\n');
for i=1:nb
    for k=1:nb
        fprintf(fid,'%-10.3f\t',Yp(i,k));
    end
    fprintf(fid,'\n\n');
end
fprintf(fid,'The positive sequence Zbus\n\n');
for i=1:nb
    for k=1:nb
        fprintf(fid,'%-10.3f\t',Zp(i,k));
    end
    fprintf(fid,'\n\n');
end

fprintf(fid,'The zero sequence Ybus\n\n');
for i=1:nb
    for k=1:nb
        fprintf(fid,'%-10.3f\t',Yz(i,k));
    end
    fprintf(fid,'\n\n');
end
fprintf(fid,'The zero sequence Zbus\n\n');
for i=1:nb
    for k=1:nb
        fprintf(fid,'%-10.3f\t',Zz(i,k));
    end
    fprintf(fid,'\n\n');
end

fprintf(fid,'Fault current\n\n');
fprintf(fid,'Magnitude   Angle(radians)\n');
fprintf(fid,'I1\n');
fprintf(fid,'%-10.3f %-10.3f\n\n',abs(ifp),angle(ifp));
fprintf(fid,'I2\n');
fprintf(fid,'%-10.3f %-10.3f\n\n',abs(ifn),angle(ifn));
fprintf(fid,'I0\n');
fprintf(fid,'%-10.3f %-10.3f\n\n',abs(ifz),angle(ifz));
fprintf(fid,'Printing of sequence components\n');
fprintf(fid,'Post Fault Voltages\n\n');
fprintf(fid,'Bus-number  Pos-seq-Magnitude  Pos-seq-Angle(radians)  Neg-seq-Magnitude  Neg-seq-Angle(radians)  Zero-seq-Magnitude  Zero-seq-Angle(radians)\n');
for i=1:length(vpf)
    fprintf(fid,'%-10.3f  %-10.3f  %-10.3f  %-10.3f  %-10.3f  %-10.3f  %-10.3f\n',i,abs(vpf(i)),vpf_phase(i),abs(vnf(i)),vnf_phase(i),abs(vzf(i)),angle(vzf(i)));
end
fprintf(fid,'\n\n');
fprintf(fid,'Post fault flows\n\n');
fprintf(fid,'The extra bus represents the reference bus\n');
fprintf(fid,'Start-bus End-bus pos-Magnitude pos-Phase neg-mag neg-phase zero-mag zero-phase\n');
fprintf(fid,'This represents the I to J flows\n\n');
for i=1:length(posflowp)
    fprintf(fid,'%-10.3f %-10.3f %-10.3f %-10.3f %-10.3f %-10.3f %-10.3f %-10.3f\n',start_bus(i),end_bus(i),abs(posflowp(i)),angle(posflowp(i)),abs(posflown(i)),angle(posflown(i)),abs(posflowz(i)),angle(posflowz(i)));
end
fprintf(fid,'\n\n');
fprintf(fid,'This represents the J to I flows\n\n');
fprintf(fid,'Start-bus End-bus pos-Magnitude pos-Phase neg-mag neg-phase zero-mag zero-phase\n');
kk=find(end_bus==(nb+1));
for i=1:length(negflowp)
    test=find(kk==i);
    if length(test)==0
        fprintf(fid,'%-10.3f %-10.3f %-10.3f %-10.3f %-10.3f %-10.3f %-10.3f %-10.3f\n',end_bus(i),start_bus(i),abs(negflowp(i)),angle(negflowp(i)),abs(negflown(i)),angle(negflown(i)),abs(negflowz(i)),angle(negflowz(i)));
    end
end
fprintf(fid,'\n\n');

fprintf(fid,'Printing of phase quantities\n');
fprintf(fid,'Bus-number  A-Magnitude  A-Angle(radians)  B-Magnitude  B-Angle(radians)  C-Magnitude  C-Angle(radians)\n');
for i=1:length(vpf)
    fprintf(fid,'%-10.3f  %-10.3f  %-10.3f  %-10.3f  %-10.3f  %-10.3f  %-10.3f\n',i,abs(va(i)),angle(va(i)),abs(vb(i)),angle(vb(i)),abs(vc(i)),angle(vc(i)));
end
fprintf(fid,'\n\n');
fprintf(fid,'Post fault flows\n\n');
fprintf(fid,'The extra bus represents the reference bus\n');
fprintf(fid,'Start-bus End-bus A-Magnitude A-Phase B-mag B-phase C-mag C-phase\n');
fprintf(fid,'This represents the I to J flows\n\n');
for i=1:length(posflowp)
    fprintf(fid,'%-10.3f %-10.3f %-10.3f %-10.3f %-10.3f %-10.3f %-10.3f %-10.3f\n',start_bus(i),end_bus(i),abs(posflowa(i)),angle(posflowa(i)),abs(posflowb(i)),angle(posflowb(i)),abs(posflowc(i)),angle(posflowc(i)));
end
fprintf(fid,'\n\n');
fprintf(fid,'This represents the J to I flows\n\n');
fprintf(fid,'Start-bus End-bus A-Magnitude A-Phase B-mag B-phase C-mag C-phase\n');
kk=find(end_bus==(nb+1));
for i=1:length(negflowp)
    test=find(kk==i);
    if length(test)==0
        fprintf(fid,'%-10.3f %-10.3f %-10.3f %-10.3f %-10.3f %-10.3f %-10.3f %-10.3f\n',end_bus(i),start_bus(i),abs(negflowa(i)),angle(negflowa(i)),abs(negflowb(i)),angle(negflowb(i)),abs(negflowc(i)),angle(negflowc(i)));
    end
end
fprintf(fid,'\n\n');
fclose(fid);

done=1;