短路计算matlab程序.doc

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Power System Analysis' Project: Matlab Program
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
close all
clc

a=exp(j*2*pi/3);
A=[1 1 1 ;
1 a^2 a ;
1 a a^2];
Sref=900e6;
Uref=300e3;
Iref=Sref/(sqrt(3)*Uref);

% data from the different documents:
Positive;
Negative;
Zero;

% run power flow:
opt = mpoption('ENFORCE_Q_LIMS', 1);
results=runpf(case9gs,opt);
Vp=[results.bus(:,8).*(cos(pi/180*results.bus(:,9))+j*sin(pi/180*results.bus(:,9)))];

clc;

% Menu:

Bucle=1;
while (Bucle==1)
% Menu:
disp('________________________________________________________________________________________')
disp('%%%%%% Program to calculate the symetrical and unsymetrical short circuit faults %%%%%%');
disp(' 1) Symetrical fault (three phase short circuit fault)');
disp(' 2) Unsymetrical fault (A single phase to ground fault)');
disp(' 0) For exit the program');
obtion=input('Choose the option you want to be calculated by typing its number:','s');
disp('________________________________________________________________________________________')

% Calculation of the fault
switch (obtion)
%%
case '1'
disp('You Have choose: SYMETRICAL FAULT');
disp(' ');
disp('Please follow the menu');

% Ask for the fault's bus number and the Rf value:
disp(' ');
disp('*******************************************************************************')
disp('* DATA *')
disp('*******************************************************************************')
bus=input ('Please, number the bus where the three phase short circuit fault occur. Bus number:');
%bus=bus+1; % Since node zero is also analyzed, the bus number X will have a X+1 index
Rf=input('Please write the value of the arc resistance Rf per phase. Rf=');

% Calculation:
If=zeros(8,1);
If(bus)=Vp(bus)/(Zpos(bus,bus)+Rf);
Vs=Vp-Zpos*If;

% Showing the Results:
disp(' ');
disp('*******************************************************************************')
disp('* RESULTS *')
disp('*******************************************************************************')

FaultCurrent=sprintf('The short circuit currents is: %gA. with angle %g?, Iref*abs(If(bus)),180/pi*angle(If(bus)));
disp(FaultCurrent);

figure();
compass(If(bus),'blue');
hold on
compass(If(bus)*a,'red');
compass(If(bus)*a^2,'green');
text(-abs(If(bus)),-abs(If(bus)),'Short Circuit current');

disp(' ')
disp ('The Voltages at all the buses in the system are:')
disp (' _____________________________________________')
disp (' Bus V Angle[Degree]')
disp (' _____________________________________________')
VoltageTokke= sprintf(' Tokke %g V %g?,Uref*abs(Vs(1)),180/pi*angle(Vs(1)));disp(VoltageTokke)
VoltageVinje= sprintf(' Vinje %g V %g?,Uref*abs(Vs(2)),180/pi*angle(Vs(2)));disp(VoltageVinje)
VoltageSonga= sprintf(' Songa %g V %g?,Uref*abs(Vs(3)),180/pi*angle(Vs(3)));disp(VoltageSonga)
VoltageVemork= sprintf(' Vemork %g V %g?,Uref*abs(Vs(4)),180/pi*angle(Vs(4)));disp(VoltageVemork)
VoltageRjukan= sprintf(' Rjunkan %g V %g?,Uref*abs(Vs(5)),180/pi*angle(Vs(5)));disp(VoltageRjukan)
VoltageFlesaker=sprintf(' Flesaker %g V %g?,Uref*abs(Vs(6)),180/pi*angle(Vs(6)));disp(VoltageFlesaker)
VoltageTveiten= sprintf(' Tveiten %g V %g?,Uref*abs(Vs(7)),180/pi*angle(Vs(7)));disp(VoltageTveiten)
VoltageRod= sprintf(' Rod %g V %g?,Uref*abs(Vs(8)),180/pi*angle(Vs(8)));disp(VoltageRod)
disp (' _____________________________________________')

figure()
subplot(2,4,1)
compass(Vs(1),'b');
hold on
compass(Vs(1)*a,'r');
compass(Vs(1)*a^2,'g');
text(-abs(Vs(1)),-abs(Vs(1)),'Tokke`s Voltage vectors');

subplot(2,4,2)
compass(Vs(2));
hold on
compass(Vs(2)*a,'r');
compass(Vs(2)*a^2,'g');
text(-abs(Vs(2)),-abs(Vs(2)),'Vinje`s Voltage vectors');

subplot(2,4,3)
compass(Vs(3));
hold on
compass(Vs(3)*a,'r');
compass(Vs(3)*a^2,'g');
text(-abs(Vs(3)),-abs(Vs(3)),'Songa`s Voltage vectors');

subplot(2,4,4)
compass(Vs(4));
hold on
compass(Vs(4)*a,'r');
compass(Vs(4)*a^2,'g');
text(-abs(Vs(4)),-abs(Vs(4)),'Vemork`s Voltage vectors');

subplot(2,4,5)
compass(Vs(5));
hold on
compass(Vs(5)*a,'r');
compass(Vs(5)*a^2,'g');
text(-abs(Vs(5)),-abs(Vs(5)),'Rjukan`s Voltage vectors');

subplot(2,4,6)
compass(Vs(6));
hold on
compass(Vs(6)*a,'r');
compass(Vs(6)*a^2,'g');
text(-abs(Vs(6)),-abs(Vs(6)),'Flesaker`s Voltage vectors');

subplot(2,4,7)
compass(Vs(7));
hold on
compass(Vs(7)*a,'r');
compass(Vs(7)*a^2,'g');
text(-abs(Vs(7)),-abs(Vs(7)),'Tveiten`s Voltage vectors');

subplot(2,4,8)
compass(Vs(8));
hold on
compass(Vs(8)*a,'r');
compass(Vs(8)*a^2,'g');
text(-abs(Vs(8)),-abs(Vs(8)),'Rod`s Voltage vectors');

disp(' ')
disp('*******************************************************************************')
disp('* CURRENTS *')
disp('*******************************************************************************')
Q=input ('Do you want to know any of the currents through two connected points? Y/N:','s');
k=0;
while Q=='Y'
k=k+1;
B1(k)=input('Write the first bus number:');
B2(k)=input ('Write the second bus number:');
if Ypos(B1(k),B2(k))==0
disp(' ')
disp('Warning: *************There is no connection between*************')
disp(' ')
k=k-1;
else
Ibranch(k)=(Vs(B1(k))-Vs(B2(k)))*(-Ypos(B1(k),B2(k)));

end
Q=input ('Do you want to know any more of the currents? Y/N:','s');
end

disp(' ')
disp('*******************************************************************************')
disp('* CURRENTS RESULTS *')
disp('*******************************************************************************')
disp('')
if k==0
disp('There are no currents results');
end
for K=1:k
Current=sprintf ('From bus no %g to bus no %g the current flowing is: %gA %g?,B1(K),B2(K),Iref*abs(Ibranch(K)),180/pi*angle(Ibranch(K)));disp (Current)

figure();
compass(Ibranch(K),'blue');
hold on
compass(Ibranch(K)*a,'red');
compass(Ibranch(K)*a^2,'green');
text(-abs(Ibranch(K)),-abs(Ibranch(K)),'Current between selected nodes');

end

%%
case '2'
disp('You Have choose: UNSYMETRICAL FAULT');
disp('The program only calcutes single line to ground faults, it is consider to happen in phase "a"')
disp(' ')
disp('Please follow the menu');

% Ask for the fault's bus number and the Rf value:
disp(' ');
disp('*******************************************************************************')
disp('* DATA *')
disp('*******************************************************************************')
bus=input ('Please, number the bus where the single phase to ground fault occur. Bus number:');
%bus=bus+1; % Since node zero is also analyzed, the bus number X will have a X+1 index
Rf=input('Please write the value of the arc resistance Rf per phase. Rf=');

% Calculation of the currents:
Ifphaseground=3*Vp(bus)/(Zpos(bus,bus)+Zneg(bus,bus)+Zzero(bus,bus)+3*Rf);

%Fase a
Ifa=zeros(8,1);
Ifa(bus)=Ifphaseground;

Ifa1=Ifa/3;
Ifa2=Ifa/3;
Ifa0=Ifa/3;

for k=1:8
Ifault=A*[Ifa0(k);Ifa1(k);Ifa2(k)];
IF(k,:)=Ifault';
end

% Voltages;

%Fase a:
V1=Vp-Zpos*Ifa1;
V2=-Zneg*Ifa2;
V0=-Zzero*Ifa0;

for k=1:8
U=A*[V0(k);V1(k);V2(k)]
V(k,:)=transpose(U)
end

% Showing the Results:
disp(' ');
disp('*******************************************************************************')
disp('* RESULTS *')
disp('*******************************************************************************')
disp(' ')
FaultCurrent=sprintf('The short circuit currents is: %g A. with angle %g ?, Iref*abs(Ifa(bus)),180/pi*angle(Ifa(bus)));
disp(FaultCurrent);

figure();
compass(IF(bus,1),'red');
hold on
compass(IF(bus,2),'blue');
compass(IF(bus,3),'green');
text(-abs(IF(bus,1)),-abs(IF(bus,1)),'Short Circuit current');

disp(' ')
disp ('The Voltages at all the buses in the system are:')

disp('__________________________________________________________________________________________________')
disp(' Bus Va Vb Vc');
disp('__________________________________________________________________________________________________')
VoltageTokke= sprintf('Tokke %gV and %g? %gV and %g? %gV and %g? ',Uref/sqrt(3)*abs(V(1,1)),180/pi*angle(V(1,1)),Uref/sqrt(3)*abs(V(1,2)),180/pi*angle(V(1,2)),Uref/sqrt(3)*abs(V(1,3)),180/pi*angle(V(1,3)));disp(VoltageTokke)
VoltageVinje= sprintf('Vinje %gV and %g? %gV and %g? %gV and %g? ',Uref/sqrt(3)*abs(V(1,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(2,2)),180/pi*angle(V(2,2)),Uref/sqrt(3)*abs(V(2,3)),180/pi*angle(V(2,3)));disp(VoltageVinje)
VoltageSonga= sprintf('Songa %gV and %g? %gV and %g? %gV and %g? ',Uref/sqrt(3)*abs(V(3,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(3,2)),180/pi*angle(V(3,2)),Uref/sqrt(3)*abs(V(3,3)),180/pi*angle(V(3,3)));disp(VoltageSonga)
VoltageVemork= sprintf('Vemork %gV and %g? %gV and %g? %gV and %g? ',Uref/sqrt(3)*abs(V(4,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(4,2)),180/pi*angle(V(4,2)),Uref/sqrt(3)*abs(V(4,3)),180/pi*angle(V(4,3)));disp(VoltageVemork)
VoltageRjukan= sprintf('Rjukan %gV and %g? %gV and %g? %gV and %g? ',Uref/sqrt(3)*abs(V(5,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(5,2)),180/pi*angle(V(5,2)),Uref/sqrt(3)*abs(V(5,3)),180/pi*angle(V(5,3)));disp(VoltageRjukan)
VoltageFlesaker=sprintf('Flesaker %gV and %g? %gV and %g? %gV and %g? ',Uref/sqrt(3)*abs(V(6,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(6,2)),180/pi*angle(V(6,2)),Uref/sqrt(3)*abs(V(6,3)),180/pi*angle(V(6,3)));disp(VoltageFlesaker)
VoltageTveiten= sprintf('Tveiten %gV and %g? %gV and %g? %gV and %g? ',Uref/sqrt(3)*abs(V(7,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(7,2)),180/pi*angle(V(7,2)),Uref/sqrt(3)*abs(V(7,3)),180/pi*angle(V(7,3)));disp(VoltageTveiten)
VoltageRod= sprintf('Rod %gV and %g? %gV and %g? %gV and %g? ',Uref/sqrt(3)*abs(V(8,1)),180/pi*angle(V(2,1)),Uref/sqrt(3)*abs(V(8,2)),180/pi*angle(V(8,2)),Uref/sqrt(3)*abs(V(8,3)),180/pi*angle(V(8,3)));disp(VoltageRod)
disp('__________________________________________________________________________________________________')

figure()
subplot(2,4,1)
compass(V(1,2),'b');
hold on
compass(V(1,1),'r');
compass(V(1,3),'g');
text(-abs(V(1,2)),-abs(V(1,2)),'Tokke`s Voltage vectors');

subplot(2,4,2)
compass(V(2,2));
hold on
compass(V(2,1),'r');
compass(V(2,3),'g');
text(-abs(V(2,2)),-abs(V(2,2)),'Vinje`s Voltage vectors');

subplot(2,4,3)
compass(V(3,2));
hold on
compass(V(3,1),'r');
compass(V(4,3),'g');
text(-abs(V(3,2)),-abs(V(3,2)),'Songa`s Voltage vectors');

subplot(2,4,4)
compass(V(4,2));
hold on
compass(V(4,1),'r');
compass(V(4,3),'g');
text(-abs(V(4,2)),-abs(V(4,2)),'Vemork`s Voltage vectors');

subplot(2,4,5)
compass(V(5,2));
hold on
compass(V(5,1),'r');
compass(V(5,3),'g');
text(-abs(V(5,2)),-abs(V(5,2)),'Rjukan`s Voltage vectors');

subplot(2,4,6)
compass(V(6,2));
hold on
compass(V(6,1),'r');
compass(V(6,3),'g');
text(-abs(V(6,2)),-abs(V(6,2)),'Flesaker`s Voltage vectors');

subplot(2,4,7)
compass(V(7,2));
hold on
compass(V(7,1),'r');
compass(V(7,3),'g');
text(-abs(V(7,2)),-abs(V(7,3)),'Tveiten`s Voltage vectors');

subplot(2,4,8)
compass(V(8,2));
hold on
compass(V(8,1),'r');
compass(V(8,3),'g');
text(-abs(V(8,2)),-abs(V(8,2)),'Rod`s Voltage vectors');

disp(' ')
disp('*******************************************************************************')