%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Title: Example 12.3 (page 752)
% By: Sonny Lloyd
% Copyright: 2009
%
% Text Reference: Power System Stability and Control by Prabha Kundur
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear,clc,close all; % clear memory and workspace
s = tf('s'); % specifies the transfer function H(s) = s (Laplace variable)
digits(4); % define precision, useful for symbolic math operations
deg2rad = pi/180; % conversion factor from degs to radians
rad2deg = 1/deg2rad; % conversion factor from radians to deg
%% Objectives:
% The objective of this example is to analyze the small-signal stability
% characterics of the system about the steady-state opeating condition
% following the loss of circuit 2.
%% Assumptions:
% The effects of amortisseurs may be neglected
% Neglect saliency of rotor
%% Operating Conditions:
Pt = 0.9; %[w] active power supplied by generator
Qt = 0.3; %[var] reactive power supplied by generator, overexcited)
Et_mag = 1; %[Vll] magnitude, line-to-line terminal voltage
Et_angle = 36; %[deg]
Et=Et_mag*exp(j*deg2rad*(Et_angle)); % polar notation
% NOTE to get Et in rectangular form: abs(Et), rad2deg*(angle(Et))
EB_mag = 0.995; %{Vll] magnitudeline-to-line bus voltage
EB_angle = 0; %[deg]
EB=EB_mag*exp(j*deg2rad*(EB_angle)); % polar notation
%% A 555MVA, 60Hz turbine generator has the following parameters:
% GIVEN SPECIFICAITON:
S = 555e6; %[VA] total power base
f0 = 60; %[Hz] frequency
Xd=1.81; % note that in PU system, Ld = Xd
Xq=1.76;
Xp_d=0.3; %X'd
Xl=0.16;
Ra = 0.003;
Tp_d0 = 8; %[s]
H=3.5;
KD=0;
% The above parameters are unsaturated values
% The effect of saturation is to be represented by assuming that d and q
% axes have similar saturation characteristics with:
Asat =0.031; Bsat = 6.93; YTI = 0.8; % based on definitions of section 3.8.2
% NOTE: giving the saturation characteristics in this manor means you have
% to CALCULATE saturation coefficients Ksd, Ksq.
%% Network Connection
% Generator is connected to an infinite bus through a step-up transformer
% (Xtrans = j0.15) and two parallel transmission lines (Xl1 = j0.5, Xl2 = 0.093)
% assuming lossless conditions.
Xtrans = 0.15; %[ohm] impedance of step-up transformer
Xl1 = 0.5; %[ohm] impedance of line 1
Xl2 = 0.93; %[ohm] impedance of line 2
% Compute the equivalent network (single line)
Xe = Xtrans + Xl1; %[ohm] equivalent network impedance, recall line 2 is out of service
Re = 0; %[ohm] assume no losses
%% The per unit fundamental parameters
% elements of d- and q-axis equivalent circuits of the equivalent generator
%
% First compute the unsaturated mutual inductances (pg. 154)
Ladu = Xd - Xl; % mutual inductance, unsaturated value
Laqu = Xq - Xl; % mutual inductance, unsaturated value
% Compute the rotor leakage inductances from the expressions for transient
% and subtransient inductances (pg. 154).
% Based on equation 4.29, the expression for L'd based on the classical
% definition is: L'd = Ll + Ladu*Lfd/(Ladu+Lfd);
% note L'd =X'dm & Ll = Xl -- L'd is a GIVEN VALUE (and is the unsaturated value)
% we need to solve for Lfd (unsaturated rotor field inductance)
% therefore: L'd(Ladu+Lfd) = Ll*(Ladu + Lfd)+Ladu*Lfd
% Lfd(L'd-Ll-Ladu) = -L'd*Ladu + Ll*Ladu
Lfd = (-Xp_d*Ladu + Xl*Ladu)/(Xp_d-Xl-Ladu); % unsaturated value, rotor field inductance
Lp_adu = 1/((1/Ladu)+1/(Lfd)); % L'ad = X'ad mutual unsaturated value %%%%%%%% is this even used anywhere? %%%%%%%
%
% Next we compute the rotor resistances from the expressions for the
% open-circuit time constants: Equ. 4.15 and 4.23
% T'd0 = (Ladu+Lfd)/Rfd [pu]
% NOTE: The time constant in per unit is equal to 377 times the time
% constant in seconds (pg 155)
Rfd = (Ladu + Lfd)/(377*Tp_d0); %[pu]
%% Initial Steady-State Values:
St = Pt + j*Qt; %[VA] total power at generator output in rect. notation
St_mag = abs(St); %[VA] magnitude of total power
St_angle = rad2deg*(angle(St)); %[deg] angle
It = conj(St)/conj(Et); % [A] phase current. Take conjugate to recognize that current is positive out of generator
It_mag = abs(It);
It_angle = rad2deg*(angle(It)); % [deg] phase current angle
%
% SOLVE for saturated coefficients of Ksq, Ksd
Yat0 = abs(Et + It*(Ra+j*Xl));
YIO = Asat*exp(Bsat*(Yat0-YTI));
Ksd = Yat0/(Yat0+YIO); % saturation coefficient, same as value given on pg. 753
Ksq = Ksd; % saturation coefficient, same as value given on pg. 753
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Find the saturated values for parameters
Lads = Ksd * Ladu; % Ld mutual inductance, saturated value
Laqs = Ksq * Laqu; % Lq mutual inductance, saturated value
Lp_ads = 1/((1/Lads)+(1/Lfd)); % saturated value inductance of L'd (see eq. 12.91)
Lds = Lads+Xl; % the SATURATED euiv. of Ld or Xd (previously defined)
Lqs = Laqs+Xl; % the SATURATED euiv. of Lq or Xq (previously defined)
%%%Lp_ds = 1/((1/Lds)+(1/Lfd)); ???????????
Xads = Lads; % equivalenet when in per unit
Xds = Lds; % equivalenet when in per unit
Xqs = Lqs; % equivalenet when in per unit
%Xp_ds = Lp_ds; ???????????
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The network constraint equation
% As a check, we know Ebus (EB) as it was given, 0.9950 <0.
% Since, Et = EB + (RE + jXE)*It
EB_check = Et - (Re + j*Xe)*It;
EB_mag_check = abs(EB_check);
EB_angle_check = rad2deg*(angle(EB_check));
%
% Calculate the internal rotor angle (Si). See Figure 3.23, pg 101.
%Check calculation of Si_angle using test data from text page 103
%Xqs = 1.494; It_mag = 1; Et_mag = 1; St_angle = 26 % test #'s page 103
Si = rad2deg*(atan((Xqs*It_mag*cos(deg2rad*(St_angle))-Ra*It_mag*sin(deg2rad*(St_angle)))/(Et_mag+Ra*It_mag*cos(deg2rad*(St_angle))+Xqs*It_mag*sin(deg2rad*(St_angle)))));
% Another way of calculating the internal rotor angle:
% Decompose into d-q components, see figure 12.6
Etq = Et + (Ra + j*Xqs)*It;
Etq_mag = abs(Etq);
Etq_angle = rad2deg*(angle(Etq));
Si_angle_verify = Etq_angle - Et_angle; % [deg] internal angle (see Section 3.6.3), between Et and q-axis
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate Initial Conditions, as outlined on page 746, solutions given on pg. 753
ed0 = real(Et_mag*sin(deg2rad*(Si))); % real part of |Et|*sin(Si)
eq0 = real(Et_mag*cos(deg2rad*(Si))); % real part of |Et|*cos(Si)
id0 = real(It_mag*sin(deg2rad*(Si+St_angle))); % real part of |It|*sin(Si+phi)
iq0 = real(It_mag*cos(deg2rad*(Si+St_angle))); % real part of |It|*cos(Si+phi)
%
% Solve for S0, equations on page 746, solutions given on pg. 753
EBd0 = ed0 - Re*id0+Xe*iq0;
EBq0 = eq0 - Re*iq0-Xe*id0;
S0 = rad2deg*(atan(EBd0/EBq0));
%
% Solve for Efd0, equations on page 746, solutions given on pg. 753
%EB2 = sqrt((EBd0^2+EBq0^2)); % NOTE: this is same value as given for EB
ifd0 = (eq0+Ra*iq0+Lds*id0)/Lads; % as per pg 746
Efd0 = Ladu*ifd0; % as per pg 746
Yad0 = Lads*(-id0+ifd0); % as per pg 746
Yaq0 = -Laqs*iq0; % as per pg 746
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Small Signal Analysis
% "Since we are expressing small-signal performance in terms of perturbed
% values of flux linkages and currents, a distinction has to be made etween
% TOTAL saturation "s" and INCREMENTAL saturation "i". Incremental
% saturation is associated with perturbed values of flux linkages and
% currents."
% Solve for Incremental values, Ksd(incr) & Ksq(incr)
Ksdi = 1/(1+Bsat*Asat*exp(Bsat*(Yat0 - YTI))); % incremental saturation, verify on page 753
Ksqi = Ksdi;
Ladi = Ksdi*Ladu; % Incremental mutual inductances
Laqi = Ksdi*Laqu; % Incremental mutual inductances
Ldi = Ladi + Xl; % Incremental inductance
Lqi = Laqi + Xl; % Incremental inductance
Lp_adi = 1/((1/Ladi)+(1/Lfd)); % Incremental mutual inductance of L'd
Xadi = Ladi; % Equivalent in per unit system
Xaqi = Laqi; % Equivalent in per unit system
Xdi = Ldi; % Equivalent in per unit system
Xqi = Lqi; % Equivalent in per unit system
Xp_di = Lp_adi; % Equivalent in per unit system
%
%
% Preparing to linearize.. see page 742
% NOTE: use of "incremental" values, as per pg.745
RT = Ra+Re;
XTq = Xe+Xqi; % Text shows Xqs .. but here we use "i", the incremental value
XTd = Xe+Xp_di+Xl; % Text shows Xqs .. but here we use "i", the incremental value
D = RT^2+XTq*XTd;
% Linearized system equations (pg. 742)
m1 = (EB_mag*(XTq*sin(deg2rad*(S0)) - RT*cos(deg2rad*(S0))))/D;
n1 = (EB_mag*(RT*sin(deg2rad*(S0)) + XTd*cos(deg2rad*(S0))))/D;
m2 = XTq/D * Ladi/(Ladi + Lfd);
n2 = RT/D * Ladi/(Ladi+Lfd);
% %%%% SUMMARY: Initial Steady-state Values of the System Varaibles
% disp('SUMMARY: Initial Steady-state Values of the System Varaibles:')
% Ksd
% Ksq
% Si
% ed0
% eq0
% id0
% iq0
% S0
% Efd0
% Ksdi
% Ksqi
% m1
% m2
% n1
% n2
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The Constants associated with the block diagram (fig. 12.9)
K1 = n1*(Yad0+Laqi*id0)-m1*(Yaq0+Lp_adi*iq0); % eq. 12.113 NOTE: incremental values used
K2 = n2*(Yad0 + Laqi*id0)-m2*(Yaq0+Lp_adi*iq0)+Lp_adi/Lfd*iq0; % eq. 12.114 NOTE: incremental values used
%% Construct State-Space Matrix
a11 = -KD/(2*H);
a12 = -K1/(2*H);
a13 = -K2/(2*H);
a21 = 2*pi*f0;
a22 = 0;
a23 = 0;
a31 = 0;
a32 = -((2*pi*f0)*Rfd/Lfd)*m1*Lp_adi; % NOTE: use of incremental value
a33 = -((2*pi*f0)*Rfd/Lfd)*(1-Lp_adi/Lfd + m2*Lp_adi); % NOTE: use of incremental value
% Construct the A-matrix
A = [a11 a12 a13; a21 a22 a22; a31 a32 a33]; % A-matrix of SMIB
%disp('The A-matrix:'); disp(A); %Displays the A-matrix to workspace
%
b11 = 1/(2*H);
b12 = 0;
b21 = 0;
b22 = 0;
b31 = 0;
b32 = (2*pi*f0)*Rfd/Ladu; %%%%% Why is this unsaturated mutual inductance Ld?
% The Constants associated with the block diagram (fig. 12.9)
K3 = -b32/a33;
K4 = -a32/b32;
T3 = -1/a33;
%% Eigenvalues
% %%%%%%%%%%%%%%%%%%%%%%%%% TEST from pg 734 %%%%%%%%%%%%%%%%%%%%%%%%
A = [(-1.43) -0.108; 377 0]; % with KD = 10 %%%%%%%%%%%%%%%%%%%%%%%%%%%
%
[eigen_Vec, eigen_Val] = eig(A); % eigen vectors and eigenvalues, where the eigen values are in a diagonal matrix
disp('The Eignvalues:'); disp(diag(eigen_Val)); % Displays eigenvalues to workspace
%disp('Characteristic Equation:'); disp(poly(A)); % Displays characteristic equation
syms h; % define a symbol for lamda% At = [(-1.43) -0.108; 377 0];
Dh = diag(h)*eye(size(A)); % form a diagonal matrix of lambda, same size as A-matrix
dec= vpa((A-Dh)); % convert to decimal values
eq1 = sort(det(dec)); % characteristic equation is the determinant of (A-D), sorted
disp('Characteristic Equation:')
pretty(eq1) % displays result in a "pretty" way
%
% Define 2nd order characteristic equation variables
% syms phi wn;
% char_eq = h^2+2*phi*wn*h+wn^2;
% disp('of the form:')
% pretty(char_eq)
%
% Look at the characteristic equation enter in wn
% wn = sqrt(8.318); % [rad] undamped natural frequency
% fn = wn/(2*pi); % [Hz] undamped natural frequency
% sigma = 1.430/(2*wn); % Damping ratio
% wd = wn*sqrt(1-sigma^2); %[rad] damped frequency
% fd = wd/(2*pi); % [Hz] damped frequency
% Find Right Eigenvector and eigenvalue, check with page 735
syms phi_11 phi_21 phi_31 phi_12 phi_22 phi_32 phi_31 phi_32 phi_33; % manually determined
[right,ev]=eig(A);
disp('The Right Eigenvector Matrix:')
disp(right)
% Method 2
% For the first eigenvalue
eig1= diag(eigen_Val(1))*eye(size(A)); % (lamda * I)
eig1a= vpa((A-eig1)); % (A-lamda*I) answer convert to decimal values
phi_m1 = [phi_11; phi_21;]; % first matrix of phi variables
eig1b = eig1a * phi_m1;
% one of the eigenvecotors corresesponding to an eigenvalue has to be set
% arbitrarily, therefore let phi_21 = 1, and solve for phi_11.
%phi_21 = 1;
[phi_11]=solve(subs(eig1b(1), phi_21, 1))
%
% Similarily, for the second eigenvalue
eig2= diag(eigen_Val(2))*eye(size(A)); % (lamda * I)
eig2a= vpa((A-eig2)); % (A-lamda*I) answer convert to decimal values
phi_m2 = [phi_12; phi_22;]; % first matrix of phi variables
eig2b = eig2a * phi_m2;
% one of the eigenvecotors corresesponding to an eigenvalue has to be set
% arbitrarily, therefore let phi_22 = 1, and solve for phi_12.
%phi_22 = 1;
[phi_12]=solve(subs(eig1b(2), phi_22, 1))
%
% Find left eigenvector, check with page 736
left=inv(right);
disp('The Left Eigenvector Matrix:')
disp(left)
%% Plotting
% Plot Eigenvalues:
% figure(1)
% real_eigen=real(eig(A));
% imag_eigen=imag(eig(A));
% plot(real_eigen,imag_eigen,'*')
% xlabel ('Real')
% ylabel ('Imaginary')
% title ('System Eigen Values')
% grid
break
Identity = eye(size(A)); % returns an identity matrix the same size as A
Use [W,D] = eig(A.'); W = conj(W) to compute the left eigenvectors.