clc
clear
Fs=1; % Fault Starting Time
Fc=1.06; %Fault clearing time
%for graph plotting
h=0.02; Time=7; t=linspace(0,Time,Time*(1/h));
%for no controller
Vb=1.0; Vt=1; Xd=0.81; Xdd=0.3; Xq=0.6; Xt=0.1;
P1=0.8; Hz=50;
Xl=0.15; Th0=0.0;
Td0d=5.044; w0=2*pi*Hz; M=8; D1=0;
%initial values calculation
Xla=Xl+Xt;
Th1=asin((P1*Xla)/(Vb*Vt));
Vtd=Vt*cos(Th1);
Vtq=Vt*sin(Th1);
Vbd=Vb*cos(Th0);
Vbq=Vb*sin(Th0);
Iqq=(Vtq-Vbq)/Xla;
Idd=(Vbd-Vtd)/Xla;
%Vdd=Vtd-Xt*Iqq
%Vqq=Vtq+Idd*Xt
Edd=Vtd+Xd*Idd;
Eqq=Vtq+Iqq*Xq;
del1=atan(Eqq/Edd);
Id=Iqq*sin(del1)-Idd*cos(del1);
Iq=Iqq*cos(del1)+Idd*sin(del1);
Vd=Vtd*sin(del1)-Vtq*cos(del1);
Vq=Vtd*cos(del1)+Vtq*sin(del1);
Edm=Vd-(Iq*Xdd)-(Xq-Xdd)*Iq;
Eqm=Vq+Id*Xdd;
Ef=Eqm+(Xd-Xdd)*Id;
%P1=Vd*Id+Vq*Iq
w1=w0;
Dw1=(w1-w0);
Ic1=2*M/w0;
Pm=P1;
%RK method
%initial values for RK method
x1=del1;
x2=w1;
w0=2*pi*Hz;
x3=Eqm;
T=0;
for m=1:Time/h
if T<=Fc & T>=Fs
P1=0;
% k1 for delta, k2 for change in speed, k3 for Eq' of M/c
%First iteration########################################
k11=(x2-w0)*h;
k21=((1/Ic1)*(Pm-P1-D1*Dw1))*h;
k31=((1/Td0d)*(Ef-(Xd-Xdd)*Id-x3))*h;
%Second iteration#####################################
k12=((x2+k11/2)-w0)*h;
k22=((1/Ic1)*(Pm-P1-D1*Dw1))*h;
k32=(1/Td0d)*(Ef-(Xd-Xdd)*Id-(x3+k31/2))*h;
%Third iteration#######################################
k13=((x2+k12/2)-w0)*h;
k23=((1/Ic1)*(Pm-P1-D1*Dw1))*h;
k33=(1/Td0d)*(Ef-(Xd-Xdd)*Id-(x3+k32/2))*h;
%Fourth iteration#######################################
k14=((x2+k13)-w0)*h;
k24=((1/Ic1)*(Pm-P1-D1*Dw1))*h;
k34=(1/Td0d)*(Ef-(Xd-Xdd)*Id-(x3+k33))*h;
x1=x1+(1/6)*(k11+k12*2+2*k13+k14);
x2=x2+(1/6)*(k21+k22*2+2*k23+k24);
x3=x3+(1/6)*(k31+k32*2+2*k33+k34);
del1=x1;
w1=x2;
Eqm=x3;
Dw1=(w1-w0);
Vbd=Vb*cos(Th0);
Vbq=Vb*sin(Th0);
Ebd=Vbd*sin(del1)-Vbq*cos(del1);
Ebq=Vbd*cos(del1)-Vbq*sin(del1);
Id=(Eqm-Ebq)/(Xdd+Xt+Xl);
Iq=(Ebd-Edm)/(Xq+Xt+Xl);
Vd=(Edm+Iq*Xq);
Vq=(Eqm-Xdd*Id);
Em=sqrt(Vd*Vd+Vq*Vq);
%for Plotting graphs of different parameters######################
d11(m)=x1*180/pi;
w11(m)=x2-w0;
pg11(m)=P1;
Vt11(m)=Em;
T=T+h; %iteration step length
else
% k1 for delta, k2 for change in speed, k3 for Eq' of M/c
%First iteration########################################
k11=(x2-w0)*h;
k21=((1/Ic1)*(Pm-P1-D1*Dw1))*h;
k31=((1/Td0d)*(Ef-(Xd-Xdd)*Id-x3))*h;
%Second iteration#####################################
k12=((x2+k11/2)-w0)*h;
k22=((1/Ic1)*(Pm-P1-D1*Dw1))*h;
k32=(1/Td0d)*(Ef-(Xd-Xdd)*Id-(x3+k31/2))*h;
%Third iteration#######################################
k13=((x2+k12/2)-w0)*h;
k23=((1/Ic1)*(Pm-P1-D1*Dw1))*h;
k33=(1/Td0d)*(Ef-(Xd-Xdd)*Id-(x3+k32/2))*h;
%Fourth iteration#######################################
k14=((x2+k13)-w0)*h;
k24=((1/Ic1)*(Pm-P1-D1*Dw1))*h;
k34=(1/Td0d)*(Ef-(Xd-Xdd)*Id-(x3+k33))*h;
x1=x1+(1/6)*(k11+k12*2+2*k13+k14);
x2=x2+(1/6)*(k21+k22*2+2*k23+k24);
x3=x3+(1/6)*(k31+k32*2+2*k33+k34);
del1=x1;
w1=x2;
Eqm=x3;
Dw1=(w1-w0);
Vbd=Vb*cos(Th0);
Vbq=Vb*sin(Th0);
Ebd=Vbd*sin(del1)-Vbq*cos(del1);
Ebq=Vbd*cos(del1)-Vbq*sin(del1);
Id=(Eqm-Ebq)/(Xdd+Xt+Xl);
Iq=(Ebd-Edm)/(Xq+Xt+Xl);
Vd=(Edm+Iq*Xq);
Vq=(Eqm-Xdd*Id);
Em=sqrt(Vd*Vd+Vq*Vq);
P1=(Vd*Id+Vq*Iq);
%for Plotting graphs of different parameters######################
d11(m)=x1*180/pi;
w11(m)=x2-w0;
pg11(m)=P1;
Vt11(m)=Em;
T=T+h; %iteration step length
end
end;
%for pss controller
%PSS
KS=.08; TW=10.41; T1=0.254;
T2=0.033;
Vb=1.0; Vt=1; Xd=0.81; Xdd=0.3; Xq=0.6; Xt=0.1;
P1=0.8; Hz=50;
Xl=0.15; Th0=0.0;
Td0d=5.044; w0=2*pi*Hz; M=8; D1=0;
%Excitation system parameters AVR
VSMAX=0.2; VSMIN=-0.2; EFMAX=7.0;
EFMIN=-6.4; TR=0.076; KA=10.0;
%initial values calculation
Xla=Xl+Xt;
Th1=asin((P1*Xla)/(Vb*Vt));
Vtd=Vt*cos(Th1);
Vtq=Vt*sin(Th1);
Vbd=Vb*cos(Th0);
Vbq=Vb*sin(Th0);
Iqq=(Vtq-Vbq)/Xla;
Idd=(Vbd-Vtd)/Xla;
%Vdd=Vtd-Xt*Iqq
%Vqq=Vtq+Idd*Xt
Edd=Vtd+Xd*Idd;
Eqq=Vtq+Iqq*Xq;
del1=atan(Eqq/Edd);
Id=Iqq*sin(del1)-Idd*cos(del1);
Iq=Iqq*cos(del1)+Idd*sin(del1);
Vd=Vtd*sin(del1)-Vtq*cos(del1);
Vq=Vtd*cos(del1)+Vtq*sin(del1);
Em=sqrt(Vd*Vd+Vq*Vq);
Edm=Vd-(Iq*Xdd)-(Xq-Xdd)*Iq;
Eqm=Vq+Id*Xdd;
Ef=Eqm+(Xd-Xdd)*Id;
%P1=Vd*Id+Vq*Iq
w1=w0;
Dw1=(w1-w0);
Ic1=2*M/w0;
Pm=P1;
Ef0=Ef;
Vref=Em;
%RK method
%initial values for RK method
x1=del1;
x2=w1;
w0=2*pi*Hz;
x3=Eqm;
x4=Em;
x5=0.0;
x6=0.0;
T=0;
for m=1:Time/h
if T<=Fc & T>=Fs
p1=0;
% k1 for delta, k2 for change in speed, k3 for Eq' of M/c
%k4 for v1 of AVR, k5 for v2 of pss, k6 for v3 of pss
%First iteration########################################
k11=(x2-w0)*h;
k21=((1/Ic1)*(Pm-P1-D1*Dw1))*h;
k31=((1/Td0d)*(Ef-(Xd-Xdd)*Id-x3))*h;
k41=((Em-x4)/TR)*h;
k51=(KS*(k21/h)-x5/TW)*h;
k61=((T1*(k51/h)+x5-x6)/T2)*h;
%Second iteration#####################################
k12=((x2+k11/2)-w0)*h;
k22=((1/Ic1)*(Pm-P1-D1*Dw1))*h;
k32=(1/Td0d)*(Ef-(Xd-Xdd)*Id-(x3+k31/2))*h;
k42=((Em-(x4+k41/2))/TR)*h;
k52=(KS*(k22/h)-(x5+k51/2)/TW)*h;
k62=((T1*(k52/h)+(x5+k51/2)-(x6+k61/2))/T2)*h;
%Third iteration#######################################
k13=((x2+k12/2)-w0)*h;
k23=((1/Ic1)*(Pm-P1-D1*Dw1))*h;
k33=(1/Td0d)*(Ef-(Xd-Xdd)*Id-(x3+k32/2))*h;
k43=((Em-(x4+k42/2))/TR)*h;
k53=(KS*(k23/h)-(x5+k52/2)/TW)*h;
k63=((T1*(k53/h)+(x5+k52/2)-(x6+k62/2))/T2)*h;
%Fourth iteration#######################################
k14=((x2+k13)-w0)*h;
k24=((1/Ic1)*(Pm-P1-D1*Dw1))*h;
k34=(1/Td0d)*(Ef-(Xd-Xdd)*Id-(x3+k33))*h;
k44=((Em-(x4+k43))/TR)*h;
k54=(KS*(k24/h)-(x5+k53)/TW)*h;
k64=((T1*(k54/h)+(x5+k53)-(x6+k63))/T2)*h;
x1=x1+(1/6)*(k11+k12*2+2*k13+k14);
x2=x2+(1/6)*(k21+k22*2+2*k23+k24);
x3=x3+(1/6)*(k31+k32*2+2*k33+k34);
x4=x4+(1/6)*(k41+k42*2+2*k43+k44);
x5=x5+(1/6)*(k51+k52*2+2*k53+k54);
x6=x6+(1/6)*(k61+k62*2+2*k63+k64);
del1=x1;
w1=x2;
Eqm=x3;
Dw1=(w1-w0);
VS=x6;
if VS>VSMAX
VS=VSMAX;
end
if VS<VSMIN
VS=VSMIN;
end
x6=VS;
Ef=Ef0+KA*(Vref-x4+x6);
if Ef>EFMAX
Ef=EFMAX;
end
if VS<VSMIN
Ef=EFMIN;
end
Ebd=Vbd*sin(del1)-Vbq*cos(del1);
Ebq=Vbd*cos(del1)-Vbq*sin(del1);
Id=(Eqm-Ebq)/(Xdd+Xt+Xl);
Iq=(Ebd-Edm)/(Xq+Xt+Xl);
Vd=(Edm+Iq*Xq);
Vq=(Eqm-Xdd*Id);
Em=sqrt(Vd*Vd+Vq*Vq);
P1=0;
Q1=0;
%for Plotting graphs of different parameters######################
d41(m)=x1*180/pi;
w41(m)=x2-w0;
pg41(m)=P1;
Vt41(m)=Em;
T=T+h; %iteration step length
else
% k1 for delta, k2 for change in speed, k3 for Eq' of M/c
%k4 for v1 of AVR, k5 for v2 of pss, k6 for v3 of pss
%First iteration########################################
k11=(x2-w0)*h;
k21=((1/Ic1)*(Pm-P1-D1*Dw1))*h;
k31=((1/Td0d)*(Ef-(Xd-Xdd)*Id-x3))*h;
k41=((Em-x4)/TR)*h;
k51=(KS*(k21/h)-x5/TW)*h;
k61=((T1*(k51/h)+x5-x6)/T2)*h;
%Second iteration#####################################
k12=((x2+k11/2)-w0)*h;
k22=((1/Ic1)*(Pm-P1-D1*Dw1))*h;
k32=(1/Td0d)*(Ef-(Xd-Xdd)*Id-(x3+k31/2))*h;
k42=((Em-(x4+k41/2))/TR)*h;
k52=(KS*(k22/h)-(x5+k51/2)/TW)*h;
k62=((T1*(k52/h)+(x5+k51/2)-(x6+k61/2))/T2)*h;
%Third iteration#######################################
k13=((x2+k12/2)-w0)*h;
k23=((1/Ic1)*(Pm-P1-D1*Dw1))*h;
k33=(1/Td0d)*(Ef-(Xd-Xdd)*Id-(x3+k32/2))*h;
k43=((Em-(x4+k42/2))/TR)*h;
k53=(KS*(k23/h)-(x5+k52/2)/TW)*h;
k63=((T1*(k53/h)+(x5+k52/2)-(x6