%電力系統(tǒng)故障程序
%*******************************************************************
                     %讀入潮流結(jié)果數(shù)據(jù)
%*******************************************************************
fidflow = fopen('e:\故障程序\inputdata\flowresult.txt','r');

dataflow=fscanf(fidflow,'%d  %f  %f',[3,inf]);

fclose(fidflow);
dataflow=dataflow'
aa=dataflow(:,2)
theta=dataflow(:,3)*3.1416/180
U=aa.*(cos(theta)+sin(theta)*i)


%*******************************************************************
                     %讀入零序網(wǎng)絡(luò)數(shù)據(jù)
%*******************************************************************

n=5%系統(tǒng)節(jié)點總數(shù)
malfunction_nodenumber=1%短路點
malfunction_type=11%故障類型

a=cos(120*3.1416/180)+sin(120*3.1416/180)*i;

fid0 = fopen('e:\故障程序\inputdata\f0.txt','r');

data0=fscanf(fid0,'%d  %d  %d  %f',[4,inf]);

fclose(fid0);
data0=data0';
%*******************************************************************
Fr_i=data0(:,2);
To_j=data0(:,3);
x0=data0(:,4)*i;
branch=size(x0);

b0=ones(branch)./x0;
m=zeros(n,n);

for nb=1:branch
    if Fr_i(nb)~=0
        m(Fr_i(nb),Fr_i(nb))=m(Fr_i(nb),Fr_i(nb))+b0(nb);
    end
    if To_j(nb)~=0 
        m(To_j(nb),To_j(nb))=m(To_j(nb),To_j(nb))+b0(nb);
    end
    if Fr_i(nb)~=0&To_j(nb)~=0
        m(Fr_i(nb),To_j(nb))=-b0(nb);
        m(To_j(nb),Fr_i(nb))=-b0(nb);    
    end
end
y0=m;%零序網(wǎng)節(jié)點導納陣
z0=inv(y0);
%*******************************************************************
                    %讀入正序和負序網(wǎng)絡(luò)數(shù)據(jù)
%*******************************************************************
fid = fopen('e:\故障程序\inputdata\f1,2.txt','r');


data1=fscanf(fid,'%d  %d  %d  %f',[4,inf]);

fclose(fid);
data1=data1';
%*******************************************************************
Fr_i=data1(:,2);
To_j=data1(:,3);
x1=data1(:,4)*i;
branch=size(x1);

b1=ones(branch)./x1;
m=zeros(n,n);

for nb=1:branch
    if Fr_i(nb)~=0
        m(Fr_i(nb),Fr_i(nb))=m(Fr_i(nb),Fr_i(nb))+b1(nb);%計算自導納
    end
    if To_j(nb)~=0 
        m(To_j(nb),To_j(nb))=m(To_j(nb),To_j(nb))+b1(nb);
    end
    if Fr_i(nb)~=0&To_j(nb)~=0
        m(Fr_i(nb),To_j(nb))=-b1(nb);%計算互導納
        m(To_j(nb),Fr_i(nb))=-b1(nb);    
    end
end
y1=m;%正序網(wǎng)節(jié)點導納陣
z1=inv(y1);


%*******************************************************************計算短路點電流
U=ones(n,1);%先設(shè)為1，實際上需要從潮流計算求得
Z1=diag(z1);%正（負）序網(wǎng)節(jié)點入端阻抗
Z0=diag(z0);%零序網(wǎng)節(jié)點入端阻抗

%*********************
%給If1，If2，If0賦值
%*********************

if malfunction_type==3%發(fā)生三相短路時短路點的正序短路電流
    If1_three=U./Z1;
    If2_three=zeros(n,1);
    If0_three=zeros(n,1);

    If1=If1_three;
    If2=zeros(n,1);
    If0=zeros(n,1);
end
if malfunction_type==1%發(fā)生單相短路接地時短路點的正負零序短路電流
    If1_one0=U./(Z1+Z1+Z0);
    If2_one0=If1_one0;
    If0_one0=If1_one0;

    If1=If1_one0;
    If2=If2_one0;
    If0=If0_one0;
end
if malfunction_type==2%發(fā)生兩相短路時短路點的正負序短路電流
    If1_two=U./(Z1+Z1);
    If2_two=-If1_two;
    If0_two=zeros(n,1);

    If1=If1_two;
    If2=If2_two;
    If0=zeros(n,1);
end
if malfunction_type==11%發(fā)生兩相短路接地時短路點的正負零序短路電流
    If1_two0=U./(Z1+(Z1.*Z0./(Z1+Z0)));
    If2_two0=-If1_two0.*(Z0./(Z1+Z0));
    If0_two0=-If1_two0.*(Z1./(Z1+Z0));

    If1=If1_two0;
    If2=If2_two0;
    If0=If0_two0;
end
If1_malfunction=If1(malfunction_nodenumber)
If2_malfunction=If2(malfunction_nodenumber)
If0_malfunction=If0(malfunction_nodenumber)

If1_abs=abs(If1);
If1_angle=angle(If1)*180/3.1416;
If2_abs=abs(If2);
If2_angle=angle(If2)*180/3.1416;
If0_abs=abs(If0);
If0_angle=angle(If0)*180/3.1416;

If1_abs_malfunction=If1_abs(malfunction_nodenumber)
If1_angle_malfunction=If1_angle(malfunction_nodenumber)
If2_abs_malfunction=If2_abs(malfunction_nodenumber)
If2_angle_malfunction=If2_angle(malfunction_nodenumber)
If0_abs_malfunction=If0_abs(malfunction_nodenumber)
If0_angle_malfunction=If0_angle(malfunction_nodenumber)



%a相
Ifa=If1+If2+If0;
Ifa_abs=abs(Ifa)
Ifa_angle=angle(Ifa)*180/3.1416

%b相
Ifb=a*a*If1+a*If2+If0;
Ifb_abs=abs(Ifb)
Ifb_angle=angle(Ifb)*180/3.1416

%c相
Ifc=a*If1+a*a*If2+If0;
Ifc_abs=abs(Ifc)
Ifc_angle=angle(Ifc)*180/3.1416

%*******************************************************************計算發(fā)生短路時各節(jié)點正負零序電壓

    I=zeros(n,1);
    I(malfunction_nodenumber)=-If1(malfunction_nodenumber);
    U1=z1*I;
    U1=U1+U

    I=zeros(n,1);
    I(malfunction_nodenumber)=-If2(malfunction_nodenumber);
    U2=z1*I

    I=zeros(n,1);
    I(malfunction_nodenumber)=-If0(malfunction_nodenumber);
    U0=z0*I

U1_abs=abs(U1)
U1_angle=angle(U1)*180/3.1416 
U2_abs=abs(U2)
U2_angle=angle(U2)*180/3.1416
U0_abs=abs(U0)
U0_angle=angle(U0)*180/3.1416
    
    

%a相
Ua=U1+U2+U0;
Ua_abs=abs(Ua)
Ua_angle=angle(Ua)*180/3.1416

%b相
Ub=a*a*U1+a*U2+U0;
Ub_abs=abs(Ub)
Ub_angle=angle(Ub)*180/3.1416

%c相
Uc=a*U1+a*a*U2+U0;
Uc_abs=abs(Uc)
Uc_angle=angle(Uc)*180/3.1416

%*******************************************************************計算支路電流
%求各序支路電流
%*******************************************************************
branch=size(x1);
I_branch1=zeros(branch,1);
I_branch2=zeros(branch,1);
I_branch0=zeros(branch,1);

for k=1:branch
    if Fr_i(k)~=0&To_j(k)~=0
        I_branch1(k)=(U1(Fr_i(k))-U1(To_j(k)))*b1(k);
        I_branch2(k)=(U2(Fr_i(k))-U2(To_j(k)))*b1(k);
        if Fr_i(k)<=3&To_j(k)<=3
            I_branch0(k)=(U0(Fr_i(k))-U0(To_j(k)))*b0(k);
        end
    end
end
I_branchmatrix1=[Fr_i,To_j,I_branch1]
I_branchmatrix2=[Fr_i,To_j,I_branch2]
I_branchmatrix0=[Fr_i,To_j,I_branch0]

I_branch_1_abs=abs(I_branch1)
I_branch_1_angle=angle(I_branch1)*180/3.1416
I_branch_2_abs=abs(I_branch2)
I_branch_2_angle=angle(I_branch2)*180/3.1416
I_branch_0_abs=abs(I_branch0)
I_branch_0_angle=angle(I_branch0)*180/3.1416

%a相
I_branch_a=[Fr_i,To_j,I_branch1+I_branch2+I_branch0]
I_branch_a_abs=abs(I_branch1+I_branch2+I_branch0)
I_branch_a_angle=angle(I_branch1+I_branch2+I_branch0)*180/3.1416

%b相
I_branch_b=[Fr_i,To_j,a*a*I_branch1+a*I_branch2+I_branch0]
I_branch_b_abs=abs(a*a*I_branch1+a*I_branch2+I_branch0)
I_branch_b_angle=angle(a*a*I_branch1+a*I_branch2+I_branch0)*180/3.1416

%c相
I_branch_c=[Fr_i,To_j,a*I_branch1+a*a*I_branch2+I_branch0]
I_branch_c_abs=abs(a*I_branch1+a*a*I_branch2+I_branch0)
I_branch_c_angle=angle(a*I_branch1+a*a*I_branch2+I_branch0)*180/3.1416