Power System Stability and Control

by

 

small signal analysis of a power system

kundur_12o3_v4.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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.

Contact us