llgfault.m

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

            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;