threephfault.m

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

function done=threephfault(nb,ng,nl,nt)
%formation of positive sequence Ybus
Yp=Ypositive(nb,ng,nl,nt);
%formation of Zbus
Zbus=inv(Yp);
%initialisation of all voltage to 1 pu at an angle zero asuming no loading
v0_matrix=ones(nb);
v0=v0_matrix(:,1); %creates a coulumn vector of voltages with 1 pu voltage
%calculation of fault current
fid=fopen('faultdata.txt','r');
b=textread('faultdata.txt');
fclose(fid);
bus=b(1);
imp=b(3);
fault_current=v0(bus)/((Zbus(bus,bus)+imp)*j);
%post fault voltage calculations
for i=1:length(v0)
    vf(i)=v0(i)-Zbus(i,bus)*fault_current*j;
end
vf_mag=abs(vf);
vf_phase=angle(vf);

%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);
ele_num=1:length(x);
start_bus=x;
end_bus=y;
for i=1:length(x)
    status(i)=1;
end
for i=1:nl
    poslineflow(i)=(vf_mag(x(i))*(cos(vf_phase(x(i)))+sin(vf_phase(x(i)))*j)-vf_mag(y(i))*(cos(vf_phase(y(i)))+sin(vf_phase(y(i)))*j))/(zp(i)*j);
    neglineflow(i)=(vf_mag(y(i))*(cos(vf_phase(y(i)))+sin(vf_phase(y(i)))*j)-vf_mag(x(i))*(cos(vf_phase(x(i)))+sin(vf_phase(x(i)))*j))/(zp(i)*j);
 end
 poslineflow_mag=abs(poslineflow);
 poslineflow_phase=angle(poslineflow);
 neglineflow_mag=abs(neglineflow);
 neglineflow_phase=angle(neglineflow);

fid=fopen('transformerdata.txt','r');
bb=textread('transformerdata.txt');
fclose(fid);
x=bb(:,1);
y=bb(:,3);
z=bb(:,5);
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
    postransflow(i)=(vf_mag(x(i))*(cos(vf_phase(x(i)))+sin(vf_phase(x(i)))*j)-vf_mag(y(i))*(cos(vf_phase(y(i)))+sin(vf_phase(y(i)))*j))/(z(i)*j);
    negtransflow(i)=(vf_mag(y(i))*(cos(vf_phase(y(i)))+sin(vf_phase(y(i)))*j)-vf_mag(x(i))*(cos(vf_phase(x(i)))+sin(vf_phase(x(i)))*j))/(z(i)*j);
end
postransflow_mag=abs(postransflow);
postransflow_phase=angle(postransflow);
negtransflow_mag=abs(negtransflow);
negtransflow_phase=angle(negtransflow);

fid=fopen('gendata.txt','r');
bb=textread('gendata.txt');
fclose(fid);
x=bb(:,1);
z=bb(:,2);
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
    posgenflow(i)=(vf_mag(x(i))*(cos(vf_phase(x(i)))+sin(vf_phase(x(i)))*j)-1)/(z(i)*j);
    neggenflow(i)=-(vf_mag(x(i))*(cos(vf_phase(x(i)))+sin(vf_phase(x(i)))*j)-1)/(z(i)*j);
end
%posgenflow
%neggenflow
posgenflow_mag=abs(posgenflow);
posgenflow_phase=angle(posgenflow);
neggenflow_mag=abs(neggenflow);
neggenflow_phase=angle(neggenflow);
    
posflow_mag=[poslineflow_mag postransflow_mag posgenflow_mag];
posflow_phase=[poslineflow_phase postransflow_phase posgenflow_phase];
negflow_mag=[neglineflow_mag negtransflow_mag neggenflow_mag];
negflow_phase=[neglineflow_phase negtransflow_phase neggenflow_phase];
status
%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;
    %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
    vf_phase(i)=vf_phase(i)+M(i)*(value*pi/180);
end
for i=1:length(posflow_phase)
    posflow_phase(i)=posflow_phase(i)+M(start_bus(i))*(value*pi/180);
    negflow_phase(i)=negflow_phase(i)+M(end_bus(i))*(value*pi/180);
end
%negflow_mag
%negflow_phase

%data entry into output file
fid=fopen('output.txt','w');
fprintf(fid,'The 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 Zbus\n\n');
for i=1:nb
    for k=1:nb
        fprintf(fid,'%-10.3f\t',Zbus(i,k));
    end
    fprintf(fid,'\n\n');
end
fprintf(fid,'Fault current\n\n');
fprintf(fid,'Magnitude   Angle(radians)\n');
fprintf(fid,'%-10.3f %-10.3f\n\n',abs(fault_current),angle(fault_current));
fprintf(fid,'Post Fault Voltages\n\n');
fprintf(fid,'Bus-number  Magnitude  Angle(radians)\n');
for i=1:length(vf)
    fprintf(fid,'%-10.3f  %-10.3f  %-10.3f\n',i,vf_mag(i),vf_phase(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 Magnitude Phase\n');
fprintf(fid,'This represents the I to J flows\n\n');
for i=1:length(posflow_mag)
    fprintf(fid,'%-10.3f %-10.3f %-10.3f %-10.3f\n',start_bus(i),end_bus(i),posflow_mag(i),posflow_phase(i));
end
fprintf(fid,'\n\n');
fprintf(fid,'This represents the J to I flows\n\n');
fprintf(fid,'Start-bus End-bus Magnitude Phase\n');
kk=find(end_bus==(nb+1));
for i=1:length(negflow_mag)
    test=find(kk==i);
    if length(test)==0
        fprintf(fid,'%-10.3f %-10.3f %-10.3f %-10.3f\n',end_bus(i),start_bus(i),negflow_mag(i),negflow_phase(i));
    end
end
fprintf(fid,'\n\n');
fclose(fid);
 
done=1;