- Email: [email protected]
Singular perturbation analysis of large-scale power systems
Singular perturbation analysis of large-scale power systems J H Chow Electrical, Computer and Systems Engineering Department, Rensselaer Polytechnic institute, Troy, NY 12180, USA J R Winkelman Ford Motor Company, Dearborn, MI 48t21, USA M A Pai a n d P W S a u e r
Department of Electrical and Computer Engineering, University of Illinois. Urbana, IL 61801. USA
This paper reviews some recent results in applying singular perturbation theory to obtain simplified power system models for stability analysis and control design. The topics include synchronous machine modelling, slow coherency and dynamic equivalence of large power networks, and transient stability analysis using direct methods. The objective is to introduce power system engineers to singular perturbation methods as a tool for providing additional insights into power system dynamics of different time scales and for analytically deriving reduced models. Keywords: singular perturbations, power system dynamics, model reduction, power system stability
I. I n t r o d u c t i o n Power system dynamic analysis encompasses a wide time span of responses, ranging from lightning phenomena in microseconds to automatic generation control over periods of minutes 1. Within the time span of stability analysis, there are also time scales arising from the various speeds of responses of different devices such as synchronous machines and excitation systems 2, as well as from the interconnections within large power systems 3. Including all time scale aspects in a stability analysis of a large power system not only results in a large computational burden but also produces a voluminous amount of output, making the interpretation of the outcome time-consuming. The objective of this paper is to review some recent results from singular perturbations as a means to obtain reduced models, and to discuss their potential applications to power system problems. Singular perturbations are most conveniently performed on two-time-scale systems in the 'explicit' form 4'5
/,=f(x, z, e),
X(to)=Xo
(1) e~=g(x, z, E),
Z(to)=Zo
Received: May 1987
Vol 1 2 Number 2 April 1 990
where x is the n-dimensional vector of slow variables and e > 0 is the singular perturbation parameter multiplying the time derivative of the m-dimensional vector of fast variables, z. Not all models of power systems with time scales are in the form of (1). Thus the first task of time scale modelling is to identify e, which could be due to ratios of small and large time constants, stray and linkage inductances, weak and strong connections, etc., and to formulate physical transformations to obtain the slow and fast variables. Once the two-time-scale system is in the form of (1), an n-dimensional slow subsystem can be obtained from (1) by setting e = 0, while an m-dimensional fast subsystem can be obtained to model the fast transients of z from its 'quasi-steady' state. The subsystems can then be used for analysis in separate time scales. In this paper the above two-time-scale modelling technique is applied to several modelling problems occurring in power systems. A non-ideal transformer example illustrates how a model in the non-explicit form is transformed to an explicit singularly perturbed model, and a voltage regulator example illustrates an iterative technique to improve transfer function approximations. These simple examples serve both as an introduction to the technique and as non-trivial applications. For synchronous machine models we use the slow manifold approach to recover damping coefficients in the slow subsystem. For the study of large power systems the slow coherency approach proposes the lumping of machines which are coherent with respect to the low-frequency modes. The theory allows for various levels of approximations, depending on the study requirements. Simulation results demonstrating the approximations are included. For direct stability analysis the singular perturbation approach is to decompose the system energy into its slow and fast components. This decomposition simplifies the detection of instability between slow coherent groups of machines. The organization of the paper is as follows. Section II illustrates the singular perturbation model reduction technique on a non-ideal transformer and a voltage
0142-0615/90/020117-10/S03.00 © Butterworth & Co (Publishers) Ltd
1 17
regulator example. Section III discusses reduced modelling of synchronous machines. Section IV considers model reduction of large power systems using the slow coherency approach. The use of the reduced models from Section IV for direct stability analysis is discussed in Section V. Some subjects in these sections are discussed only briefly. Readers interested in more details are referred to the cited references.
where
A°=
I - - R 1/'L x -R~/~/(LIL2)
As a first example of the presence of time scales in individual power system components and as an illustration of the singular perturbation method, consider the power transformer. The single-phase-balanced circuit for the transformer is shown in Figure 1, where L~, L 2 are the self-inductances of the coils, M is the mutual inductance, i 1, i 2 a r e the currents through the coils and v is a slowly varying voltage. Using x~=i~, x 2 = i 2 a s the states, the equations of the transformer model are
-R2/L2
l
(6)
The presence of time scales is not obvious from (4)-(6). To exhibit the time scale behaviour explicitly, we introduce a change of coordinates y=Px,
II. Introductory singular perturbation model reduction examples
-- R 2 / ' ~ / ( L 1L 2 )~
z=Qx
(7)
where P and Q form bases for the left null and row spaces of A o respectively (see Reference 6, Chapter 3). The left null space of A 0, P=[L~
- ~/(L~L2) ]
(6)
and the row space of A o, Q = [R, R z ~ / ( L 1 / L 2 ) ]
(9)
define, according to (7), the slow variable dx 1/dt = - [ R 1L 2 / ( L 1
L2 -- M 2
)Ix 1 y= Llxl-x/(L1Lz)x2
-- [ M R 2 / ( L I L 2 -- M 2 ) ] x 2 + [L 2fiL 1L 2 -- M 2 )] v (2) dx2/dt = _ [ M R ~ / ( L 1 L 2 - M2)]xa - [R2L 1/(LIL2
+ [ M / ( L I L 2 - M2)]v In the case of an ideal transformer, L 1 L 2 = M 2. For a non-ideal transformer, L 1 L z - ~ M 2 and this leakage is modelled as M2=(1--g,,)L1L
and the fast variable Z= R1x 1 + R2N/(L 1/L2)X 2
-- M 2 ) ] x 2
2
--RjL 1
A(~:)=(l/e)[ __X,r(l _e)R1/v/(L1L2)
-V''(I -e')R2/v/(L1L2)] - R E/Lz J (4)
which, with the approximation x/(1 - e ) rewritten as
1 - e / 2 , can be
A(e) = (1/e)(A o + eA 1 + O(e 2))
(5)
([1)
Note that the slow variable y has dimensions of a flux linkage and the fast variable z has dimensions of a voltage. In the new coordinates y, z the model (2) becomes, using equations (5), (10) and (11), dy/dt = --[1/(T~ + T2)]y + [ ( T , - T2)/2(T , + T2)]z + v/2 +O(~ 2 )
(3)
where e is a small dimensionless parameter. Using (3), the system dynamic matrix of (2) becomes
(lO)
dz/dt =
1/T~)/2(T1 + T2)]y - [ l / L + 1/r2-e/(rl + Tg]z + ( l / T 1 + 1/T2)vq-O(g2)
e[(1/T 2 -
(12)
where T~ = L ~ / R 1 and T 2 = L 2 / R 2 a r e time constants of the primary and secondary R L circuits respectively. From an examination of the fast state equation in (12), we see that there is an order-of-one presence of the slow forcing voltage v. Thus the fast state z will be composed of a fast part superimposed on a slowly varying part. The slow part of z will influence the slow subsystem dynamics because of the order-of-one coupling into the slow state equations. The model for the slow subsystem
d~/dt = [-- 1/(T, + T2)]~ + [T,/(T 1 + T2)]v
(13)
M R1 •
•
+ "----~ i 2
+
v_()
v1
L1
L2
v2
•R •
2
is obtained by setting e = 0 in the second equation of (12) and substituting the quasi-steady-state of the fast subsystem, £ = v, into the first equation. To obtain the model for the fast subsystem in the fast time scale, we scale time such that t=ev. Changing the time scale in the second equation, the fast model
d~/dz = Figure 1. A non-ideal transformer 1 18
-
(1/T 1 + I/T2)~
(14)
is obtained by setting e = O. Let us give a physical interpretation to the new state
Electrical Power & Energy Systems
variables y, z. For the slow variable y we can write Efd y = L I X 1 -- 4 ( L 1 L 2 ) x 2 = L I X 1 -- M x 2 + 0 ( ~ )
= kIJ11 -F kI/12 -~- O(e )
(15)
(
We see that y is, to order e, the total flux linkage in the primary coil, where W~I is due to the current in the primary coil and W~2 is due to the current in the secondary coil. This total flux linkage is a slowly varying variable. Using (11 ) and writing the fast variable in terms of its fast and slowly varying parts, we obtain Z=2 + ~=V+'2=R1x
(16)
a+ R2~/(L1/L2)x 2
Rearranging (16) and solving for the secondary voltage v 2, we obtain
Figure 2. Linearized IEEE type 1 exciter or
v2 = R z x 2 = x / ( L 2 / L l ) ( v - R l X l +2)
= (N2/Nx)(v a +_~)
(17)
where N1 and N 2 are the numbers of turns in the primary and secondary coils, respectively. After ~, owing to initial conditions, has decayed, v2=(N2/N1)v 1, which is the voltage relationship for an ideal transformer. We illustrate with this analysis how simplified models may be systematically obtained for components which exhibit two time scales. Interested readers may wish to contrast this approach with an empirical two-time-scale analysis of the non-ideal transformer in Reference 7. As a second example of how singular perturbation techniques may be used in analysis of power system equipment, we consider the simplification of the transfer function of an IEEE type 1 exciter s. The small-signal model, neglecting saturation, is shown in Figure 2, where VR,Efd and Rf are the states associated with the regulator, field voltage and the rate feedback circuit respectively. Based on typical values of parameters, 1 / K , = e and Ta = e/c~are small, where e is a scaling constant. Using this fact, we write the state equations as
"~"
-1 ]~f t
e~
~ Efd/ 'I G/T~
1/T~
I I I [--KE/TEK ~ I
0
I _O~KF/TF
I
o-
AET c( I D
Substituting (20) into the slow equation in (19), we get the following slow subsystem model: x=(All -A12A2~(O)Azl)2-A12A2)(O)B2
J
Vol 12 Number 2 April 1990
AET
(21)
Using (20) and (21) to solve for the field voltage in terms of the terminal voltage error, we obtain the following approximation to the exciter transfer function:
Erd -- 1K,( 1\ + sr "} aET
(22)
From (22) we see that the low-frequency characteristic of the transfer function of the full model is poorly approximated. Equation (22) exhibits an integral characteristic at low frequency while the full-order model exhibits a lead/lag-type characteristic. To improve the transfer function approximation, we make a first-order correction to the slow subsystem; that is, we include an e, correction term in the slow subsystem. Expanding
Rf I
A221= ( A 2o2 ( I + . ( A 2 2o)
(23)
1A22)) 1 -1
= (I + e.(A°2) - 'A~2)-'(A°2) -1 -~ (I - e(A°2) 1A~2)- 1(A2°2)-1
1/T E --0~
(24)
Making the appropriate substitutions into (19), we have, after some algebraic manipulations, the improved approximation to the exciter transfer function
.VR/
T~
OI
(20)
we obtain, to a first order of e,
m
-07
+
5 = - A221(0)(A 2,2 + B 2 AET)
A22(~:) = A°2 +/3A12+ O(~2)
~VR%_I -
where e is a small parameter. If we let e --* 0 in (19), we get for the quasi-steady-state value of the fast subsystem
(18)
+-K#~vs / K E} AE T
(25)
We see from (25) that the first-order correction of the transfer function approximation exhibits the lead/lag characteristic of the full-order model transfer function characteristic.
119
To test the accuracy of (25), we assume the following typical values for the linearized model parameters:
for i
TF= 1.0 s Kv=0.16 p.u./p.u. TE =0.5 s K. = 100 p.u./p.u.
1. . . . . n,
6ji=c~j--3 i
(32)
Traditional approaches to reduced-order modelling recognize that T~oiis large while Tqol is small compared to the dynamics of the mechanical system 3~, ~oi. Following this reasoning, the reduced multi-machine electromechanical model with constant field flux (rEdoi ~ ,-f~) and no damper windings (Tqoi ~ O) gives
Ta = 0.05 s
K E= 1.0 p.u./p.u. The slow eigenvalue for the full-order system is -0.057 rad s 1, the transfer function zero is - 1.0 rad sand the DC gain is 100p.u./p.u. From (25) the slow eigenvalue is - 0 . 0 6 2 rad s- 1, the transfer function zero is -1.0 rad s 1 and the DC gain is 94 p.u./p.u. We see from this example how reduced-order models of varying degrees of approximation may be developed using series expansions. This basic concept is utilized in the next section to develop non-linear reduced-order models of synchronous machines.
Efdi=.~)((~ 1. . . . .
0n)
(33)
and d3i
~-=
(Db ((D i - - 1)
(34)
dcD i
(35)
2H i ~[-=gi(31 . . . . . 3,)
where the n fi functions of6 a. . . . . 6, are the solution of the n linear (in Eh) algebraic equations
III. R e d u c e d - o r d e r
machine modelling
While the synchronous machine has been widely studied, its rich time scale properties have yet to be fully exploited. The natural time scale separations which exist between various electrical and mechanical subsystems allow for systematic model reduction. In this section we present some recent results on model reduction of synchronous machines. For the non-linear synchronous machine models, the theory of integral manifolds 9 seems best suited for their analysis. While integral manifolds are not easily found (and may not exist) for general non-linear systems, they frequently can be found systematically in multi-time-scale systems. This is illustrated in a multi-machine example and used to justify the assumption of inter-machine damping terms in reduced-order models. In commonly used notations, an n-machine power system model including the effect of field flux decay and damper windings (two-axis model of Anderson and Fouad 1°) is , dE'qi
T~toi ~ -dE'di _
7q°' dt
E'q,--(La,-- Ld,)Id, + Efd ~ ,
Ed,+(Lq,-Lqi)lq,
(26)
E'di =(Lqi-Eql) ~ Yq[E'q°jcos(6j, +Vo) +Ehj sin(fji +71)] j=l
(36) This multi-machine reduced-order model is undamped since the effect of taking T~o~ as zero is equivalent to omitting the damper windings. This is also clear from the absence of speed terms in (35). The earlier result of single-machine model reduction using an integral manifold approach 11 is extended to show how the inter-machine damping in a multi-machine system can be systematically restored without resorting to the addition of the damper-winding differential equations. Considering each T~oias a constant (ci) times a small parameter E, we seek the multi-dimensional integral manifold 9
E~ = hi(f1 . . . . . 6,, o9t . . . . . ~o,, ~)
(37)
Each hi must be found from the partial differential equations obtained by substituting (37) into (27) as k = l ~ k (,0b((Ok - l ) - ' ] - ~ k ~ g k
(27)
=
--hi+(Lqi--Cqi) (38)
d8 i
~ - = ~)b(~)i-- 1)
• Yij[E'q°j cos(6ji + 70)
(28)
+
hj sin(6ji + 71j)]
j=l
d°Ji _ 2H~ dt
,
,
(Edildi+Eqilql)+ Tmi=gi
(29)
For hi represented in a power series in e as
hi=hio + ehil +8,2hi2 + ' ' "
where
ldi= ~ j=l
Yij[E',xjcos(fji+?ii)-E',dsin(3ji+?q)]
(39)
evaluation at e = 0 gives
(30) 0
=
-- hlo +
(Lq,-
Lqi) ~'. Yo[E'q~cos(6ji + ?q) j=l
Iqi= ~ j=l
120
Yij[E'q,iCOS(,Sji+?ii)+E'
(31)
+ hjo sin(3ji + 7q)]
(40)
Electrical Power & Energy Systems
Solving these n linear (in hjo) algebraic equations gives n non-linear functions of only 31 . . . . . 6. as
hio = f i ( 6 , . . . . . 6,)
(41)
In order to find higher-order terms for the integral manifold (37), we keep the e terms but neglect the e 2 and higher-order terms to obtain
connections between machines in the same area. The parameter ~ is small when the external connections are high-impedance lines 6 or when the external connections are sparse compared to the internal connections t*. Let 6i be the rotor angle of the ith machine with respect to the synchronous frame of reference, J , be the index set of all machines i in area c~, and ~ ' = summation over all machines iE J~
(k=~1 (')hio = -hi~ +(Lq~--Lq~) ~, Yuhj~ sin(6ji+7 u) j=l
i
(42)
Again, solving these n linear (in hj~) algebraic equations gives n non-linear functions of 3~ . . . . . 3, and co~. . . . . co,, written symbolically as
hii =.fil (bl . . . . . 6,, col . . . . . co,)
d3, - ' = cob(coi- 1) dt
(44)
2Hi ~
(45)
while the actual use of the integral manifold non-linear expressions may be impractical for multi-machine dynamic studies, the damping contribution may be approximated by a linear term such as Duo) j where the Dq are constants evaluated from a linearization about an operating point. This illustration shows how the time scale separation due to small Tqo~ can be exploited to improve on traditional reduced-order models and give a basis for previously heuristic techniques of accounting for inter-machine damping.
IV. Slow coherency and dynamic equivalencing In realistic large-scale power systems, groups of machines tend to oscillate at frequencies slower than the oscillations between individual close-by machines. Reduced models constructed on the basis of this physical consideration would allow the study of the slow dynamics due to inter-area power transfer. In this section we present a slow coherency approach to obtain two-time-scale reduced models of power systems 6'~ 2,13. The approach consists of two parts: the first is the identification of the slow coherent machines, i.e. which machines are to be grouped as an aggregate; the second is the construction of reduced models. IV.1 Identification of s l o w coherent m a c h i n e s To illustrate slow coherency, we consider an n-machine power system consisting o f r areas, and denote by e > 0 the ratio of the strength of the external connections between machines in different areas to the strength of internal
Vol 12 Number 2 April 1990
mi d2 31/dt 2 = P m i - Vi2Gii -- 2 ~ Vi VjBij sin 6ij J
- ~ ~PeV~VjBijsin3~
(43)
Substitution of (37) using (39), (41) and (43) into (28)-(32) gives the improved reduced-order model which reflects the fundamental inter-machine damping terms visible by the presence of COgin Gi:
= Gi(61 . . . . . 6,, (Ol . . . . . co.)
Normalizing the external connections by separating out the small parameter e, the swing equations of the system can be written as
[~-1 j
for all iEd=, ~ t = l , 2 . . . . . r
(46)
where mi is the inertia, Pmi the input mechanical power, V~ the voltage behind the transient reactance of machine i, Bi~ is the transfer admittance between machines i and j, G, represents the load and 6~j=6~-6j. For simplicity, damping and transfer conductance terms are neglected in (46). To show that (46) has two time scales, we linearize (46) about an equilibrium point &e to obtain in matrix form M~ = (K; + e,KE)x
(47)
where x = 6 - 6 e, M is the diagonal inertia matrix and K ~ and K E are the internal and external connection matrices of synchronizing coefficients V~ViBi j cos 6~j between machines in the same areas and in different areas respectively. An important property is that entries in every row of K t and K E sum to zero. Since the time scales are not explicit in (47), a transformation based on the null space and the row space of K ~is used 6. We introduce the centre of angle variables
m i x i / ~ mi,
y,=~ i
'
~ = 1,2 . . . . . r
(48a)
i
which are the inertia-weighted angles, and the difference variables
z~, = x i - x j,
k = 1,2,..., n,- 1 for all i ~ J,, i =#j and ct = 1,2 . . . . . r
(48b)
where j is the reference machine in area ct and n, is the number of machines in area cc Denoting U as an n x r grouping matrix whose (i, ~) entry is 1 if i~ J , and 0 otherwise, the transformation (48) becomes
121
where z is comprised of the subvectors z ", thus defining G, and Ca = (UTMU) 1UTM = M a 1UTM. The inverse of T is T I=[U
G +]
(50)
where G + = M - 1GT(GM -
1G T) - 1 =
M - 1 GTMd 1
Applying the transformation (49) to the linearized model (47), and using the property K~U=0, we obtain
IMaMd
lI:l
01[~]=F~Ka ~Sad Le K d a K d + ~ K d d
(51)
where K.
= UTKEU
Kad = UTKEM - 1GTMd Kda -- Kad m
T
K d = M d G M - 1KIM- 1GTMd Kdd = MdGM- 1KEM - IGTMd For e small, (51) is a two-time-scale system where the centres of angles y are the slow variables and the difference variables z are the fast variables. In addition, the grouping matrix U is an O(e) approximation of the slow eigensubspace of (47). Since the rows of U are equal for machines in the same area, i.e. the mode shapes are identical, the motions of these machines with respect to the slow modes are the same. This phenomenon is called 'slow coherency'. This slow coherency behaviour, which is due to the weak and strong connections, is also present in the non-linear model (46). In realistic power systems, slow coherency is not exact and the slow eigensubspace is only an approximate grouping matrix. This property can be used to identify the slow coherent machines if they are not known a priori. The grouping algorithm 12 computed the slow eigensubspace and performs column and row operations on the eigenvector matrix to obtain an approximation to a grouping matrix, which assigns the machines into coherent groups. For large power systems the Lanczos algorithm is used to exploit sparsity by including the load buses in computing the slow eigensubspace 15. The grouping algorithm has been applied to power systems as large as 4000 buses and 600 machines. As an illustration, the grouping algorithm was applied to a 411-machine/1750-bus model of the US West Coast System (WSCC). A partition into 11 areas, obtained by using the 11 lowest-frequency mode shapes, is shown in Figure 3. The dots in Figure 3 denote the locations of some of the machines. The partition consists of eight large coherent groups, each corresponding to a load and generation centre. The three small groups are machines connected to the main system through relatively weak ties. The lowest-frequency mode is computed to be at 0.34 Hz, which has been measured on the system 3. Using this mode and the system mode, the WSCC system was partitioned into a North System and a South System (Figure 4), which agrees with the observed north-south swing pattern.
122
Figure 3. An 11-area partition of WSCC
Figure 4. A two-area partition of WSCC
Electrical Power & Energy Systems
IV.2 Reduced-order models To exhibit the time scales in the non-linear model (46), we apply the same slow coherency transformation (48)
" /Z= ya-= ~ mi&i i
~
100
," ~ E
V
mi
i
(52)
z~=6~-@
k = l , 2 . . . . . n ~ - l , i~J~,, i # j
oop
for all ct= 1,2 . . . . . r, to obtain a non-linear singularly perturbed system in explicit form
it
z')
el = ol(z ~) + @ ~ ( y , , z ~. . . . . z ' )
(53)
40
I'~
(54)
where f,, g~, 0~ are non-linear functions of the sine type. In (53), (54) the area variables y, are slow while the local difference variables z" are fast. We note that the coupling of the slow variables into the fast variables and the coupling between the fast variables from different areas are weak, which offers many possibilities for obtaining reduced models for power system dynamic simulations and direct stability analysis. The latter subject will be discussed in the next section. In power system dynamic simulations, if only the slow inter-area dynamics is of interest, then the fast variables z" are assumed to be at their quasi-steady states and solved from gkl ~ , ~ ~+ ~-~ _ ~1 . . . . . 2'r) = 0 gk(Y=,
%
ff,,,f /,,,c
) 3~: = ~ f : ( y : , z ~. . . . .
\
20 1
'~.~.~.."s\t "x~l= I
0
I
2 Time
I
I
3
4
5
(s)
Figure 5. Machine angles--exact and approximation A1
100[
(55)
Substituting U, which is only a function of .~, into (53), the slow variables are approximated by the zero-order model
(56) We shall denote (55), (56) as approximation A1. Equation (55) can be interpreted as a load flow equation yielding the slowly varying angles between machines in the same areas. Modelling these time-varying angles as phase shifters, the slow variable ~, can be regarded as the motion of an equivalent single machine connected to individual coherent machine buses through these phase shifters. If, in addition to the slow dynamics, the fast dynamics of the study areas are also of interest, then only the fast variables outside of the study area are assumed to be at their quasi-steady state and solved from (55) (approximation A2). A further approximation is to assume that the fast variables outside of the study area are constant, that is, the phase shifters are fixed (approximation A3). In this case, no algebraic equations of the type (55) are needed. These three simplifications are illustrated with the simulation of a 48-machine system 13. Figures 5-7 compare the approximations with the exact response (E) for several machines in the coherent group containing the disturbance. Approximation A1 captures the slow behaviour of the machine response whereas A2 corrects for the fast behaviour as well. The loss of accuracy of approximation A3 is small. IV.3 Large-power-system application For large power systems, approximation A3 can be implemented using the EPRI dynamic equivalencing
Vol 12 Number 2 April 1990
? °°if'l,,;
x
! ,.J X
20
\?. ,,,, A2
0
I
1
I
I
2
3
I
4
5
Time (s) Figure 6. Machine angles--exact and approximation A2
programs 16. These programs generate reduced models with physical parameters suitable for use in conventional power system dynamic simulation programs. To use this approach, the slow coherent areas are identified using the grouping algorithm. A coherent area that does not need to be studied in detail is aggregated into a single equivalent machine using the EPRI programs. If only the inter-area dynamics is of interest, then all the areas are represented by equivalent machines. For extremely large systems the procedure may have to be repeated several times to obtain
123
100
---,
E\r"
For convenience we define the difference between the angle of machine i in area ~ and 3.,~as (57)
zi = 6~ - 3'~ 80
such that 2ij : 2 i
-
if
2 j = (~ij
-
i,j~J=
(58)
6O
Z i j = t~ij
'/
O3 t-
<
40
/,'
!,\\%, ~, ~\
,
IIII/,.,'/ ~,~ '.,..\ \\~/ex\,?~_/,//
/
",',
6 /
d2)?==P=-e k ~ = ~
20
\x
0
if i ~ J = , j c J ~ , O~#fl
Y~t[3
-
V.1 S l o w global s u b s y s t e m The slow subsystem representing the inter-area power transfer is obtained from (46) by using (48a) and (57) and then setting the fast variables to their initial condition:
,,yl/Y
,,
-
#=l
iI
I
1
I
I
I
2
3
4
Time
Figure 7. Machine angles A3
i
ViV~B,jsin(~=~+z~j)
(59)
j
where
=E (Pml- Vi2Gii)
{s)
P=
exact and approximation
(60)
i
and the superscript 'o' denotes initial values of the variables. a reduced model of appropriate size. For instance, a coherent area may contain many machines, not all of which are of interest. Then an interim reduced model is obtained by aggregating only the external coherent areas. The interim reduced model allows the use of the grouping algorithm to identify smaller coherent groups. The reducing procedure using slow coherent areas has been applied to several 2000-bus power system models. An example is a stability investigation of several units in a 211-machine/1078-bus system. For this problem the external coherent areas were reduced to single-machine equivalents. The number of machines in the study area was determined to be too large, and thus the coherent area identification process is applied to the study area also, allowing for the reduction of the study area. Since local stability is of interest, the system around the study units was kept in detail, as well as major transmission facilities within the study area. The resulting reduced model consists of 18 machines and 67 buses. A comparison of the response of the full model (dashed curve) and the reduced model (solid curve) is shown in Figure 8. In the machine angle and speed plots the dashed curves are the response of the machines in a coherent group within the study area. The solid curve is the equivalent machine of the coherent group, which tracks the slow behaviour and does not contain any fast inter-machine-type oscillations, as predicted by theory.
V.2 Fast local s u b s y s t e m s The fast subsystem representing the dynamics of machines within an area is obtained from (46) by using
100
E
=~
," " ' - ~iL~-'='~--_=-.
50
""
cxJ
~~'---~..I--~"~-'------
-------'--- ~- -===
0 f
-50 I .015 "O
g
1.olo
o.
II
II
z~ :~
,e..~ ~
• •
1.000 0.995 1000
750 O~
V. Transient stability analysis In direct or energy methods of transient stability analysis it is crucial to know the energy threshold for instability in order to obtain an accurate stability assessment. For larger power systems with slow inter-area dynamics the system energy can be decomposed into a slow global energy and fast local energies, thus allowing for more accurate energy monitoring ~7. In this section we present the ideas using the zeroth-order models from (53), (54).
124
~ -~_o - ~
500 250
>~-
0 -250 0
20
40
60 Time
80
100
120
(cycles)
Figure 8. The 211 -machine/1078-bus example
Electrical Power & Energy Systems
conditions(57) and :setting the slow variables at their initial d2-~i
~
~--,a
m i ~ i ~ = l ' m i -- ~.
J
Vi V j B i j
sin
V~ast, =PE =
at
_~
p~(~_~o)
+~.= ~ ViVjBij(cos~ij-cos2i°j)
~ij
i
j=i+
(67)
1
It can readily be shown that
~, ~ ViVjBijsin(yia+ ~.,j)
-e
r KE Vsll°w -~- 2 ~=1
fl=l j
i = t , 2 . . . . . n~,c~= 1,2 . . . . . r
V.3 Slow and fast energies The energy function associated with the slow global subsystem (59) is called the slow transient energy and is given by (62)
Here the slow kinetic energy
V KE _1 ~ ~2 -~low-9 /.2. M~yat
(63)
~CC=I
and the slow potential energy PE ~ 7,0 V~low = - P a t ( 9 ~ - - yat) at=l
--F, ~ ZatEfl ViVjBij ffl~ot ill j x [cos( ~,~ + z°fl - cos(~°a +
z°j)]
(64)
where 9, = ~, -.Vo with r
t r
Note that Yo is the global centre of angle. When the area variables are thus defined relative to this global centre of angle, expressions for kinetic and potential energy similar to those found in Athay et al. as will be obtained. Similarly, the energy function associated with the fast subsystem of area c~is called the fast transient energy and is given by Vfast PE . , at = V ; KE a s t , a t + Vlast,a
(65)
where V~ast,:~
1
at
M,~ 2
Vol 12 Number 2 April 1990
KE
(68)
(61)
Since subsystems (60) and (61) are decoupled, we can decompose the system transient energy function into a slow component and r fast components. The advantage of denoting the fast variables relative to the centre of an area, instead of the difference variables used in Section IV, is that the transient energy function for both slow and fast subsystems has the same form as the system-wide transient energy function ~8,~ 9
Vslow= VslKo~w+ V~orw
KE
v;.s,.= = V~ys,om
(66)
is the total transient kinetic energy of the system. However, because of the decoupling of the fast and slow subsystems due to the 'freezing' of the variables in the slow and fast subsystems (59) and (61) respectively, we cannot make a similar statement regarding the potential energy. The following argument will be used to motivate the use of either the slow or the fast energy to determine the critical clearing time for a given disturbance. Suppose for a given fault on the system that the slow and fast dynamics get disturbed, but that the instability of the global system is due to the separation of the areas (inter-area instability mode) and not due to instability inside one of the areas (local instability mode). The total system energy function will not differentiate between energies due to separation of the areas (which are relevant in this case) and energies due to oscillations of the machines inside each area (which are irrelevant in this case), i.e. the total system energy will add the two energies as if both energies were relevant for this disturbance. Consequently, this may lead to erroneous conclusions about the system stability. However, the slow energy is the inter-area energy and each fast energy is the energy within an area. Hence, for the class of faults where the instability is due to the disturbance of the inter-area dynamics, only the slow energy is relevant and should be used to determine the critical clearing time. The above concept has been rigorously justified via the concept of integral manifolds, where corrections to the zeroth-order slow energy function are successively made 2°. Similarly, when a disturbance is such that the fast subsystem severely impacts the slow subsystem (e.g. a fault near a generator), we can introduce a coordinate transformation which reflects the effect of the boundary layer correction term on the slow subsystem as in the case of the linear dynamic system 21. Thus for a large class of disturbances the slow subsystem can be used to predict critical clearing times by using any of the techniques in direct methods. Numerical results are contained in References 17 and 20. It is conjectured that the kinetic and potential energy corrections to the energy function found in the literature 22'23 are equivalent to applying this methodology. Research in testing the system on large power networks is continuing.
Vl. Conclusions The results described in this paper illustrate the role of singular perturbations as an analytical tool in power system modelling at both the device and system level. They also illustrate how additional insights into power system dynamics can be obtained in the process. Research is continuing on improved reduced-order models of machines and excitation systems, more accurate aggregate models of groups of machines, and the use of reduced models for transient stability analysis. Additional
125
effort is needed to develop procedures and computer programs ready for use by the power industry.
11
VII. A c k n o w l e d g e m e n t s
12
This work is supported in part by the US Department of Energy under contract DE-AC0!-84CE76249 and in part by the National Science Foundation under Grant ECS-84-14677.
VIII. References 1
2
3 4 5 6 7 8
9 10
126
Concordia, C and Schulz, R P "Appropriate component representation for the simulation of power system dynamics' Syrup. on Adequacy and Philosophy of Modeling." Dynamic System Performance, IEEE Pub/. 75CHO970-4-PWR (1 975) Cate, E G, H e m m a p l a r d h , K, M a n k e , J W and Gelopulos, D P 'Time frame notion and time response of the models in transient, mid-term and long-term stability programs" /EEE Trans. PowerAppar. & Syst. Vol PAS-103 (1984) pp 143-151 Cresap, R Land Hauer, J F'Emergenceofanewswingmode in the Western Power System' IEEE Trans. Power Appar. & Syst. Vol PAS-100 (1981 ) pp 2037 2045 K o k o t o v i c , PV, O ' M a l l e y , R EJrand Sannuti, P'Singular perturbations and order reduction in control theory--an overview" Automatica Vol 12 (1 976) pp 123-1 32 Saksena, V K, O'Reilly, J and K o k o t o v i c , P V 'Singular perturbations and time-scale methods in control theory: survey 1976-1983" Automatica Vol 20 (1984) pp 273-293 C h o w , J H (ed.) Time-Scale Modeling of Dynamic Networks with Applications to Power Systems," Lecture Notes in Control and Information Sciences Vol 46, Springer, New York (1 982) Fitzgerald, A E and Kingsley, C Electric Machinery 2nd edn, McGraw-Hill, New York (1961) IEEE Committee Report Excitation system models for power system stability studies" IEEE Trans. Power Appar. & Syst. Vol 100 (1981 ) pp 494-509 Carr, J Applications of Center Manifold Theory Springer, New York (1981) Anderson, P i and Fouad, A A Power System Control and Stability Iowa State University Press, Ames (1 977)
13 14 15
16 17
18 19
20 21 22
23
K o k o t o v i c , P V and Sauer, P W 'Integral manifold as a tool for reduced order modeling of nonlinear systems: a case study' Proc. 1986 IEEE Int. Symp. on Circuits and Systems (1986) pp 908 911 K o k o t o v i c , P V, A v r a m o v i c , B, C h o w , J H and W i n k e l m a n , J R 'Coherency-based decomposition and aggregation' Automatica Vol 18 (1 982) pp 47-56 W i n k e l m a n , J R et al. 'An analysis of inter-area dynamics of multi-machine systems' IEEE Trans. Power Appar. & Syst. Vol PAS-100 (t 981 ) pp 754-763 C h o w , J H and K o k o t o v i c , P V 'Time scale modeling of sparse dynamic networks' IEEE Trans. Automatic Control Vol AC-30 (1985) pp714-722 C h o w , J H, Cullum, J and Willoughby, R A "A sparsity-based technique for identifying slow-coherent areas in large power systems' IEEE Trans. Power Appar. & Syst. Vol PAS-103 (1984) pp463-473 'Development of dynamic equivalents for transient stability studies' EPRI Report EL-456 (1 977) Pal, i A, O t h m a n , H, C h o w , J H and W i n k e l m a n , J R 'Time scale energies in power systems' Proc. 4th IFAC/IFORS Syrup. on Large Scale Systems Theory and Applications, Zurich (26 29 August 1986) A t h a y , T, Podmore, R and V i r m a n i , S A practical method for direct analysis of transient stability' IEEE Trans. Power Appar. & Syst. Vol PAS-102 (1979) pp 573-594 Ribbens-Pavella, M and Evans, F J 'Direct methods in the study of the dynamics of large-scale electric power systems--a survey' Proc. 8th IFAC World Cong., Kyoto (24-28 August 1981 ) Pal, a . A. et al. 'Decoupled stability analysis of large scale power systems using integral manifold approach' Proc. 25th IEEE Conf on Decision and Control, Athens (1 986) Chang, K W 'A diagonalization technique for singular perturbation problems' Proc. 12th AIlerton Conf. on Circuit and System Theory, University of Illinois, Urbana (1 974) Fouad, A A and V i t t a l , A 'Power system response to a large disturbance: energy associated with system separation' IEEE Trans. Power Appar. & Syst. Vol PAS-102 (1 983) pp 3534-3540 Y i - x i n Ni, Bun Kee Wee and Fouad, A A 'Accuracy of transient energy functions' Proc. 1986 North American Power Symp., Ithaca, New York (1 3-14 October 1986)
Electrical Power & Energy Systems
Department of Electrical and Computer Engineering, University of Illinois. Urbana, IL 61801. USA
This paper reviews some recent results in applying singular perturbation theory to obtain simplified power system models for stability analysis and control design. The topics include synchronous machine modelling, slow coherency and dynamic equivalence of large power networks, and transient stability analysis using direct methods. The objective is to introduce power system engineers to singular perturbation methods as a tool for providing additional insights into power system dynamics of different time scales and for analytically deriving reduced models. Keywords: singular perturbations, power system dynamics, model reduction, power system stability
I. I n t r o d u c t i o n Power system dynamic analysis encompasses a wide time span of responses, ranging from lightning phenomena in microseconds to automatic generation control over periods of minutes 1. Within the time span of stability analysis, there are also time scales arising from the various speeds of responses of different devices such as synchronous machines and excitation systems 2, as well as from the interconnections within large power systems 3. Including all time scale aspects in a stability analysis of a large power system not only results in a large computational burden but also produces a voluminous amount of output, making the interpretation of the outcome time-consuming. The objective of this paper is to review some recent results from singular perturbations as a means to obtain reduced models, and to discuss their potential applications to power system problems. Singular perturbations are most conveniently performed on two-time-scale systems in the 'explicit' form 4'5
/,=f(x, z, e),
X(to)=Xo
(1) e~=g(x, z, E),
Z(to)=Zo
Received: May 1987
Vol 1 2 Number 2 April 1 990
where x is the n-dimensional vector of slow variables and e > 0 is the singular perturbation parameter multiplying the time derivative of the m-dimensional vector of fast variables, z. Not all models of power systems with time scales are in the form of (1). Thus the first task of time scale modelling is to identify e, which could be due to ratios of small and large time constants, stray and linkage inductances, weak and strong connections, etc., and to formulate physical transformations to obtain the slow and fast variables. Once the two-time-scale system is in the form of (1), an n-dimensional slow subsystem can be obtained from (1) by setting e = 0, while an m-dimensional fast subsystem can be obtained to model the fast transients of z from its 'quasi-steady' state. The subsystems can then be used for analysis in separate time scales. In this paper the above two-time-scale modelling technique is applied to several modelling problems occurring in power systems. A non-ideal transformer example illustrates how a model in the non-explicit form is transformed to an explicit singularly perturbed model, and a voltage regulator example illustrates an iterative technique to improve transfer function approximations. These simple examples serve both as an introduction to the technique and as non-trivial applications. For synchronous machine models we use the slow manifold approach to recover damping coefficients in the slow subsystem. For the study of large power systems the slow coherency approach proposes the lumping of machines which are coherent with respect to the low-frequency modes. The theory allows for various levels of approximations, depending on the study requirements. Simulation results demonstrating the approximations are included. For direct stability analysis the singular perturbation approach is to decompose the system energy into its slow and fast components. This decomposition simplifies the detection of instability between slow coherent groups of machines. The organization of the paper is as follows. Section II illustrates the singular perturbation model reduction technique on a non-ideal transformer and a voltage
0142-0615/90/020117-10/S03.00 © Butterworth & Co (Publishers) Ltd
1 17
regulator example. Section III discusses reduced modelling of synchronous machines. Section IV considers model reduction of large power systems using the slow coherency approach. The use of the reduced models from Section IV for direct stability analysis is discussed in Section V. Some subjects in these sections are discussed only briefly. Readers interested in more details are referred to the cited references.
where
A°=
I - - R 1/'L x -R~/~/(LIL2)
As a first example of the presence of time scales in individual power system components and as an illustration of the singular perturbation method, consider the power transformer. The single-phase-balanced circuit for the transformer is shown in Figure 1, where L~, L 2 are the self-inductances of the coils, M is the mutual inductance, i 1, i 2 a r e the currents through the coils and v is a slowly varying voltage. Using x~=i~, x 2 = i 2 a s the states, the equations of the transformer model are
-R2/L2
l
(6)
The presence of time scales is not obvious from (4)-(6). To exhibit the time scale behaviour explicitly, we introduce a change of coordinates y=Px,
II. Introductory singular perturbation model reduction examples
-- R 2 / ' ~ / ( L 1L 2 )~
z=Qx
(7)
where P and Q form bases for the left null and row spaces of A o respectively (see Reference 6, Chapter 3). The left null space of A 0, P=[L~
- ~/(L~L2) ]
(6)
and the row space of A o, Q = [R, R z ~ / ( L 1 / L 2 ) ]
(9)
define, according to (7), the slow variable dx 1/dt = - [ R 1L 2 / ( L 1
L2 -- M 2
)Ix 1 y= Llxl-x/(L1Lz)x2
-- [ M R 2 / ( L I L 2 -- M 2 ) ] x 2 + [L 2fiL 1L 2 -- M 2 )] v (2) dx2/dt = _ [ M R ~ / ( L 1 L 2 - M2)]xa - [R2L 1/(LIL2
+ [ M / ( L I L 2 - M2)]v In the case of an ideal transformer, L 1 L 2 = M 2. For a non-ideal transformer, L 1 L z - ~ M 2 and this leakage is modelled as M2=(1--g,,)L1L
and the fast variable Z= R1x 1 + R2N/(L 1/L2)X 2
-- M 2 ) ] x 2
2
--RjL 1
A(~:)=(l/e)[ __X,r(l _e)R1/v/(L1L2)
-V''(I -e')R2/v/(L1L2)] - R E/Lz J (4)
which, with the approximation x/(1 - e ) rewritten as
1 - e / 2 , can be
A(e) = (1/e)(A o + eA 1 + O(e 2))
(5)
([1)
Note that the slow variable y has dimensions of a flux linkage and the fast variable z has dimensions of a voltage. In the new coordinates y, z the model (2) becomes, using equations (5), (10) and (11), dy/dt = --[1/(T~ + T2)]y + [ ( T , - T2)/2(T , + T2)]z + v/2 +O(~ 2 )
(3)
where e is a small dimensionless parameter. Using (3), the system dynamic matrix of (2) becomes
(lO)
dz/dt =
1/T~)/2(T1 + T2)]y - [ l / L + 1/r2-e/(rl + Tg]z + ( l / T 1 + 1/T2)vq-O(g2)
e[(1/T 2 -
(12)
where T~ = L ~ / R 1 and T 2 = L 2 / R 2 a r e time constants of the primary and secondary R L circuits respectively. From an examination of the fast state equation in (12), we see that there is an order-of-one presence of the slow forcing voltage v. Thus the fast state z will be composed of a fast part superimposed on a slowly varying part. The slow part of z will influence the slow subsystem dynamics because of the order-of-one coupling into the slow state equations. The model for the slow subsystem
d~/dt = [-- 1/(T, + T2)]~ + [T,/(T 1 + T2)]v
(13)
M R1 •
•
+ "----~ i 2
+
v_()
v1
L1
L2
v2
•R •
2
is obtained by setting e = 0 in the second equation of (12) and substituting the quasi-steady-state of the fast subsystem, £ = v, into the first equation. To obtain the model for the fast subsystem in the fast time scale, we scale time such that t=ev. Changing the time scale in the second equation, the fast model
d~/dz = Figure 1. A non-ideal transformer 1 18
-
(1/T 1 + I/T2)~
(14)
is obtained by setting e = O. Let us give a physical interpretation to the new state
Electrical Power & Energy Systems
variables y, z. For the slow variable y we can write Efd y = L I X 1 -- 4 ( L 1 L 2 ) x 2 = L I X 1 -- M x 2 + 0 ( ~ )
= kIJ11 -F kI/12 -~- O(e )
(15)
(
We see that y is, to order e, the total flux linkage in the primary coil, where W~I is due to the current in the primary coil and W~2 is due to the current in the secondary coil. This total flux linkage is a slowly varying variable. Using (11 ) and writing the fast variable in terms of its fast and slowly varying parts, we obtain Z=2 + ~=V+'2=R1x
(16)
a+ R2~/(L1/L2)x 2
Rearranging (16) and solving for the secondary voltage v 2, we obtain
Figure 2. Linearized IEEE type 1 exciter or
v2 = R z x 2 = x / ( L 2 / L l ) ( v - R l X l +2)
= (N2/Nx)(v a +_~)
(17)
where N1 and N 2 are the numbers of turns in the primary and secondary coils, respectively. After ~, owing to initial conditions, has decayed, v2=(N2/N1)v 1, which is the voltage relationship for an ideal transformer. We illustrate with this analysis how simplified models may be systematically obtained for components which exhibit two time scales. Interested readers may wish to contrast this approach with an empirical two-time-scale analysis of the non-ideal transformer in Reference 7. As a second example of how singular perturbation techniques may be used in analysis of power system equipment, we consider the simplification of the transfer function of an IEEE type 1 exciter s. The small-signal model, neglecting saturation, is shown in Figure 2, where VR,Efd and Rf are the states associated with the regulator, field voltage and the rate feedback circuit respectively. Based on typical values of parameters, 1 / K , = e and Ta = e/c~are small, where e is a scaling constant. Using this fact, we write the state equations as
"~"
-1 ]~f t
e~
~ Efd/ 'I G/T~
1/T~
I I I [--KE/TEK ~ I
0
I _O~KF/TF
I
o-
AET c( I D
Substituting (20) into the slow equation in (19), we get the following slow subsystem model: x=(All -A12A2~(O)Azl)2-A12A2)(O)B2
J
Vol 12 Number 2 April 1990
AET
(21)
Using (20) and (21) to solve for the field voltage in terms of the terminal voltage error, we obtain the following approximation to the exciter transfer function:
Erd -- 1K,( 1\ + sr "} aET
(22)
From (22) we see that the low-frequency characteristic of the transfer function of the full model is poorly approximated. Equation (22) exhibits an integral characteristic at low frequency while the full-order model exhibits a lead/lag-type characteristic. To improve the transfer function approximation, we make a first-order correction to the slow subsystem; that is, we include an e, correction term in the slow subsystem. Expanding
Rf I
A221= ( A 2o2 ( I + . ( A 2 2o)
(23)
1A22)) 1 -1
= (I + e.(A°2) - 'A~2)-'(A°2) -1 -~ (I - e(A°2) 1A~2)- 1(A2°2)-1
1/T E --0~
(24)
Making the appropriate substitutions into (19), we have, after some algebraic manipulations, the improved approximation to the exciter transfer function
.VR/
T~
OI
(20)
we obtain, to a first order of e,
m
-07
+
5 = - A221(0)(A 2,2 + B 2 AET)
A22(~:) = A°2 +/3A12+ O(~2)
~VR%_I -
where e is a small parameter. If we let e --* 0 in (19), we get for the quasi-steady-state value of the fast subsystem
(18)
+-K#~vs / K E} AE T
(25)
We see from (25) that the first-order correction of the transfer function approximation exhibits the lead/lag characteristic of the full-order model transfer function characteristic.
119
To test the accuracy of (25), we assume the following typical values for the linearized model parameters:
for i
TF= 1.0 s Kv=0.16 p.u./p.u. TE =0.5 s K. = 100 p.u./p.u.
1. . . . . n,
6ji=c~j--3 i
(32)
Traditional approaches to reduced-order modelling recognize that T~oiis large while Tqol is small compared to the dynamics of the mechanical system 3~, ~oi. Following this reasoning, the reduced multi-machine electromechanical model with constant field flux (rEdoi ~ ,-f~) and no damper windings (Tqoi ~ O) gives
Ta = 0.05 s
K E= 1.0 p.u./p.u. The slow eigenvalue for the full-order system is -0.057 rad s 1, the transfer function zero is - 1.0 rad sand the DC gain is 100p.u./p.u. From (25) the slow eigenvalue is - 0 . 0 6 2 rad s- 1, the transfer function zero is -1.0 rad s 1 and the DC gain is 94 p.u./p.u. We see from this example how reduced-order models of varying degrees of approximation may be developed using series expansions. This basic concept is utilized in the next section to develop non-linear reduced-order models of synchronous machines.
Efdi=.~)((~ 1. . . . .
0n)
(33)
and d3i
~-=
(Db ((D i - - 1)
(34)
dcD i
(35)
2H i ~[-=gi(31 . . . . . 3,)
where the n fi functions of6 a. . . . . 6, are the solution of the n linear (in Eh) algebraic equations
III. R e d u c e d - o r d e r
machine modelling
While the synchronous machine has been widely studied, its rich time scale properties have yet to be fully exploited. The natural time scale separations which exist between various electrical and mechanical subsystems allow for systematic model reduction. In this section we present some recent results on model reduction of synchronous machines. For the non-linear synchronous machine models, the theory of integral manifolds 9 seems best suited for their analysis. While integral manifolds are not easily found (and may not exist) for general non-linear systems, they frequently can be found systematically in multi-time-scale systems. This is illustrated in a multi-machine example and used to justify the assumption of inter-machine damping terms in reduced-order models. In commonly used notations, an n-machine power system model including the effect of field flux decay and damper windings (two-axis model of Anderson and Fouad 1°) is , dE'qi
T~toi ~ -dE'di _
7q°' dt
E'q,--(La,-- Ld,)Id, + Efd ~ ,
Ed,+(Lq,-Lqi)lq,
(26)
E'di =(Lqi-Eql) ~ Yq[E'q°jcos(6j, +Vo) +Ehj sin(fji +71)] j=l
(36) This multi-machine reduced-order model is undamped since the effect of taking T~o~ as zero is equivalent to omitting the damper windings. This is also clear from the absence of speed terms in (35). The earlier result of single-machine model reduction using an integral manifold approach 11 is extended to show how the inter-machine damping in a multi-machine system can be systematically restored without resorting to the addition of the damper-winding differential equations. Considering each T~oias a constant (ci) times a small parameter E, we seek the multi-dimensional integral manifold 9
E~ = hi(f1 . . . . . 6,, o9t . . . . . ~o,, ~)
(37)
Each hi must be found from the partial differential equations obtained by substituting (37) into (27) as k = l ~ k (,0b((Ok - l ) - ' ] - ~ k ~ g k
(27)
=
--hi+(Lqi--Cqi) (38)
d8 i
~ - = ~)b(~)i-- 1)
• Yij[E'q°j cos(6ji + 70)
(28)
+
hj sin(6ji + 71j)]
j=l
d°Ji _ 2H~ dt
,
,
(Edildi+Eqilql)+ Tmi=gi
(29)
For hi represented in a power series in e as
hi=hio + ehil +8,2hi2 + ' ' "
where
ldi= ~ j=l
Yij[E',xjcos(fji+?ii)-E',dsin(3ji+?q)]
(39)
evaluation at e = 0 gives
(30) 0
=
-- hlo +
(Lq,-
Lqi) ~'. Yo[E'q~cos(6ji + ?q) j=l
Iqi= ~ j=l
120
Yij[E'q,iCOS(,Sji+?ii)+E'
(31)
+ hjo sin(3ji + 7q)]
(40)
Electrical Power & Energy Systems
Solving these n linear (in hjo) algebraic equations gives n non-linear functions of only 31 . . . . . 6. as
hio = f i ( 6 , . . . . . 6,)
(41)
In order to find higher-order terms for the integral manifold (37), we keep the e terms but neglect the e 2 and higher-order terms to obtain
connections between machines in the same area. The parameter ~ is small when the external connections are high-impedance lines 6 or when the external connections are sparse compared to the internal connections t*. Let 6i be the rotor angle of the ith machine with respect to the synchronous frame of reference, J , be the index set of all machines i in area c~, and ~ ' = summation over all machines iE J~
(k=~1 (')hio = -hi~ +(Lq~--Lq~) ~, Yuhj~ sin(6ji+7 u) j=l
i
(42)
Again, solving these n linear (in hj~) algebraic equations gives n non-linear functions of 3~ . . . . . 3, and co~. . . . . co,, written symbolically as
hii =.fil (bl . . . . . 6,, col . . . . . co,)
d3, - ' = cob(coi- 1) dt
(44)
2Hi ~
(45)
while the actual use of the integral manifold non-linear expressions may be impractical for multi-machine dynamic studies, the damping contribution may be approximated by a linear term such as Duo) j where the Dq are constants evaluated from a linearization about an operating point. This illustration shows how the time scale separation due to small Tqo~ can be exploited to improve on traditional reduced-order models and give a basis for previously heuristic techniques of accounting for inter-machine damping.
IV. Slow coherency and dynamic equivalencing In realistic large-scale power systems, groups of machines tend to oscillate at frequencies slower than the oscillations between individual close-by machines. Reduced models constructed on the basis of this physical consideration would allow the study of the slow dynamics due to inter-area power transfer. In this section we present a slow coherency approach to obtain two-time-scale reduced models of power systems 6'~ 2,13. The approach consists of two parts: the first is the identification of the slow coherent machines, i.e. which machines are to be grouped as an aggregate; the second is the construction of reduced models. IV.1 Identification of s l o w coherent m a c h i n e s To illustrate slow coherency, we consider an n-machine power system consisting o f r areas, and denote by e > 0 the ratio of the strength of the external connections between machines in different areas to the strength of internal
Vol 12 Number 2 April 1990
mi d2 31/dt 2 = P m i - Vi2Gii -- 2 ~ Vi VjBij sin 6ij J
- ~ ~PeV~VjBijsin3~
(43)
Substitution of (37) using (39), (41) and (43) into (28)-(32) gives the improved reduced-order model which reflects the fundamental inter-machine damping terms visible by the presence of COgin Gi:
= Gi(61 . . . . . 6,, (Ol . . . . . co.)
Normalizing the external connections by separating out the small parameter e, the swing equations of the system can be written as
[~-1 j
for all iEd=, ~ t = l , 2 . . . . . r
(46)
where mi is the inertia, Pmi the input mechanical power, V~ the voltage behind the transient reactance of machine i, Bi~ is the transfer admittance between machines i and j, G, represents the load and 6~j=6~-6j. For simplicity, damping and transfer conductance terms are neglected in (46). To show that (46) has two time scales, we linearize (46) about an equilibrium point &e to obtain in matrix form M~ = (K; + e,KE)x
(47)
where x = 6 - 6 e, M is the diagonal inertia matrix and K ~ and K E are the internal and external connection matrices of synchronizing coefficients V~ViBi j cos 6~j between machines in the same areas and in different areas respectively. An important property is that entries in every row of K t and K E sum to zero. Since the time scales are not explicit in (47), a transformation based on the null space and the row space of K ~is used 6. We introduce the centre of angle variables
m i x i / ~ mi,
y,=~ i
'
~ = 1,2 . . . . . r
(48a)
i
which are the inertia-weighted angles, and the difference variables
z~, = x i - x j,
k = 1,2,..., n,- 1 for all i ~ J,, i =#j and ct = 1,2 . . . . . r
(48b)
where j is the reference machine in area ct and n, is the number of machines in area cc Denoting U as an n x r grouping matrix whose (i, ~) entry is 1 if i~ J , and 0 otherwise, the transformation (48) becomes
121
where z is comprised of the subvectors z ", thus defining G, and Ca = (UTMU) 1UTM = M a 1UTM. The inverse of T is T I=[U
G +]
(50)
where G + = M - 1GT(GM -
1G T) - 1 =
M - 1 GTMd 1
Applying the transformation (49) to the linearized model (47), and using the property K~U=0, we obtain
IMaMd
lI:l
01[~]=F~Ka ~Sad Le K d a K d + ~ K d d
(51)
where K.
= UTKEU
Kad = UTKEM - 1GTMd Kda -- Kad m
T
K d = M d G M - 1KIM- 1GTMd Kdd = MdGM- 1KEM - IGTMd For e small, (51) is a two-time-scale system where the centres of angles y are the slow variables and the difference variables z are the fast variables. In addition, the grouping matrix U is an O(e) approximation of the slow eigensubspace of (47). Since the rows of U are equal for machines in the same area, i.e. the mode shapes are identical, the motions of these machines with respect to the slow modes are the same. This phenomenon is called 'slow coherency'. This slow coherency behaviour, which is due to the weak and strong connections, is also present in the non-linear model (46). In realistic power systems, slow coherency is not exact and the slow eigensubspace is only an approximate grouping matrix. This property can be used to identify the slow coherent machines if they are not known a priori. The grouping algorithm 12 computed the slow eigensubspace and performs column and row operations on the eigenvector matrix to obtain an approximation to a grouping matrix, which assigns the machines into coherent groups. For large power systems the Lanczos algorithm is used to exploit sparsity by including the load buses in computing the slow eigensubspace 15. The grouping algorithm has been applied to power systems as large as 4000 buses and 600 machines. As an illustration, the grouping algorithm was applied to a 411-machine/1750-bus model of the US West Coast System (WSCC). A partition into 11 areas, obtained by using the 11 lowest-frequency mode shapes, is shown in Figure 3. The dots in Figure 3 denote the locations of some of the machines. The partition consists of eight large coherent groups, each corresponding to a load and generation centre. The three small groups are machines connected to the main system through relatively weak ties. The lowest-frequency mode is computed to be at 0.34 Hz, which has been measured on the system 3. Using this mode and the system mode, the WSCC system was partitioned into a North System and a South System (Figure 4), which agrees with the observed north-south swing pattern.
122
Figure 3. An 11-area partition of WSCC
Figure 4. A two-area partition of WSCC
Electrical Power & Energy Systems
IV.2 Reduced-order models To exhibit the time scales in the non-linear model (46), we apply the same slow coherency transformation (48)
" /Z= ya-= ~ mi&i i
~
100
," ~ E
V
mi
i
(52)
z~=6~-@
k = l , 2 . . . . . n ~ - l , i~J~,, i # j
oop
for all ct= 1,2 . . . . . r, to obtain a non-linear singularly perturbed system in explicit form
it
z')
el = ol(z ~) + @ ~ ( y , , z ~. . . . . z ' )
(53)
40
I'~
(54)
where f,, g~, 0~ are non-linear functions of the sine type. In (53), (54) the area variables y, are slow while the local difference variables z" are fast. We note that the coupling of the slow variables into the fast variables and the coupling between the fast variables from different areas are weak, which offers many possibilities for obtaining reduced models for power system dynamic simulations and direct stability analysis. The latter subject will be discussed in the next section. In power system dynamic simulations, if only the slow inter-area dynamics is of interest, then the fast variables z" are assumed to be at their quasi-steady states and solved from gkl ~ , ~ ~+ ~-~ _ ~1 . . . . . 2'r) = 0 gk(Y=,
%
ff,,,f /,,,c
) 3~: = ~ f : ( y : , z ~. . . . .
\
20 1
'~.~.~.."s\t "x~l= I
0
I
2 Time
I
I
3
4
5
(s)
Figure 5. Machine angles--exact and approximation A1
100[
(55)
Substituting U, which is only a function of .~, into (53), the slow variables are approximated by the zero-order model
(56) We shall denote (55), (56) as approximation A1. Equation (55) can be interpreted as a load flow equation yielding the slowly varying angles between machines in the same areas. Modelling these time-varying angles as phase shifters, the slow variable ~, can be regarded as the motion of an equivalent single machine connected to individual coherent machine buses through these phase shifters. If, in addition to the slow dynamics, the fast dynamics of the study areas are also of interest, then only the fast variables outside of the study area are assumed to be at their quasi-steady state and solved from (55) (approximation A2). A further approximation is to assume that the fast variables outside of the study area are constant, that is, the phase shifters are fixed (approximation A3). In this case, no algebraic equations of the type (55) are needed. These three simplifications are illustrated with the simulation of a 48-machine system 13. Figures 5-7 compare the approximations with the exact response (E) for several machines in the coherent group containing the disturbance. Approximation A1 captures the slow behaviour of the machine response whereas A2 corrects for the fast behaviour as well. The loss of accuracy of approximation A3 is small. IV.3 Large-power-system application For large power systems, approximation A3 can be implemented using the EPRI dynamic equivalencing
Vol 12 Number 2 April 1990
? °°if'l,,;
x
! ,.J X
20
\?. ,,,, A2
0
I
1
I
I
2
3
I
4
5
Time (s) Figure 6. Machine angles--exact and approximation A2
programs 16. These programs generate reduced models with physical parameters suitable for use in conventional power system dynamic simulation programs. To use this approach, the slow coherent areas are identified using the grouping algorithm. A coherent area that does not need to be studied in detail is aggregated into a single equivalent machine using the EPRI programs. If only the inter-area dynamics is of interest, then all the areas are represented by equivalent machines. For extremely large systems the procedure may have to be repeated several times to obtain
123
100
---,
E\r"
For convenience we define the difference between the angle of machine i in area ~ and 3.,~as (57)
zi = 6~ - 3'~ 80
such that 2ij : 2 i
-
if
2 j = (~ij
-
i,j~J=
(58)
6O
Z i j = t~ij
'/
O3 t-
<
40
/,'
!,\\%, ~, ~\
,
IIII/,.,'/ ~,~ '.,..\ \\~/ex\,?~_/,//
/
",',
6 /
d2)?==P=-e k ~ = ~
20
\x
0
if i ~ J = , j c J ~ , O~#fl
Y~t[3
-
V.1 S l o w global s u b s y s t e m The slow subsystem representing the inter-area power transfer is obtained from (46) by using (48a) and (57) and then setting the fast variables to their initial condition:
,,yl/Y
,,
-
#=l
iI
I
1
I
I
I
2
3
4
Time
Figure 7. Machine angles A3
i
ViV~B,jsin(~=~+z~j)
(59)
j
where
=E (Pml- Vi2Gii)
{s)
P=
exact and approximation
(60)
i
and the superscript 'o' denotes initial values of the variables. a reduced model of appropriate size. For instance, a coherent area may contain many machines, not all of which are of interest. Then an interim reduced model is obtained by aggregating only the external coherent areas. The interim reduced model allows the use of the grouping algorithm to identify smaller coherent groups. The reducing procedure using slow coherent areas has been applied to several 2000-bus power system models. An example is a stability investigation of several units in a 211-machine/1078-bus system. For this problem the external coherent areas were reduced to single-machine equivalents. The number of machines in the study area was determined to be too large, and thus the coherent area identification process is applied to the study area also, allowing for the reduction of the study area. Since local stability is of interest, the system around the study units was kept in detail, as well as major transmission facilities within the study area. The resulting reduced model consists of 18 machines and 67 buses. A comparison of the response of the full model (dashed curve) and the reduced model (solid curve) is shown in Figure 8. In the machine angle and speed plots the dashed curves are the response of the machines in a coherent group within the study area. The solid curve is the equivalent machine of the coherent group, which tracks the slow behaviour and does not contain any fast inter-machine-type oscillations, as predicted by theory.
V.2 Fast local s u b s y s t e m s The fast subsystem representing the dynamics of machines within an area is obtained from (46) by using
100
E
=~
," " ' - ~iL~-'='~--_=-.
50
""
cxJ
~~'---~..I--~"~-'------
-------'--- ~- -===
0 f
-50 I .015 "O
g
1.olo
o.
II
II
z~ :~
,e..~ ~
• •
1.000 0.995 1000
750 O~
V. Transient stability analysis In direct or energy methods of transient stability analysis it is crucial to know the energy threshold for instability in order to obtain an accurate stability assessment. For larger power systems with slow inter-area dynamics the system energy can be decomposed into a slow global energy and fast local energies, thus allowing for more accurate energy monitoring ~7. In this section we present the ideas using the zeroth-order models from (53), (54).
124
~ -~_o - ~
500 250
>~-
0 -250 0
20
40
60 Time
80
100
120
(cycles)
Figure 8. The 211 -machine/1078-bus example
Electrical Power & Energy Systems
conditions(57) and :setting the slow variables at their initial d2-~i
~
~--,a
m i ~ i ~ = l ' m i -- ~.
J
Vi V j B i j
sin
V~ast, =PE =
at
_~
p~(~_~o)
+~.= ~ ViVjBij(cos~ij-cos2i°j)
~ij
i
j=i+
(67)
1
It can readily be shown that
~, ~ ViVjBijsin(yia+ ~.,j)
-e
r KE Vsll°w -~- 2 ~=1
fl=l j
i = t , 2 . . . . . n~,c~= 1,2 . . . . . r
V.3 Slow and fast energies The energy function associated with the slow global subsystem (59) is called the slow transient energy and is given by (62)
Here the slow kinetic energy
V KE _1 ~ ~2 -~low-9 /.2. M~yat
(63)
~CC=I
and the slow potential energy PE ~ 7,0 V~low = - P a t ( 9 ~ - - yat) at=l
--F, ~ ZatEfl ViVjBij ffl~ot ill j x [cos( ~,~ + z°fl - cos(~°a +
z°j)]
(64)
where 9, = ~, -.Vo with r
t r
Note that Yo is the global centre of angle. When the area variables are thus defined relative to this global centre of angle, expressions for kinetic and potential energy similar to those found in Athay et al. as will be obtained. Similarly, the energy function associated with the fast subsystem of area c~is called the fast transient energy and is given by Vfast PE . , at = V ; KE a s t , a t + Vlast,a
(65)
where V~ast,:~
1
at
M,~ 2
Vol 12 Number 2 April 1990
KE
(68)
(61)
Since subsystems (60) and (61) are decoupled, we can decompose the system transient energy function into a slow component and r fast components. The advantage of denoting the fast variables relative to the centre of an area, instead of the difference variables used in Section IV, is that the transient energy function for both slow and fast subsystems has the same form as the system-wide transient energy function ~8,~ 9
Vslow= VslKo~w+ V~orw
KE
v;.s,.= = V~ys,om
(66)
is the total transient kinetic energy of the system. However, because of the decoupling of the fast and slow subsystems due to the 'freezing' of the variables in the slow and fast subsystems (59) and (61) respectively, we cannot make a similar statement regarding the potential energy. The following argument will be used to motivate the use of either the slow or the fast energy to determine the critical clearing time for a given disturbance. Suppose for a given fault on the system that the slow and fast dynamics get disturbed, but that the instability of the global system is due to the separation of the areas (inter-area instability mode) and not due to instability inside one of the areas (local instability mode). The total system energy function will not differentiate between energies due to separation of the areas (which are relevant in this case) and energies due to oscillations of the machines inside each area (which are irrelevant in this case), i.e. the total system energy will add the two energies as if both energies were relevant for this disturbance. Consequently, this may lead to erroneous conclusions about the system stability. However, the slow energy is the inter-area energy and each fast energy is the energy within an area. Hence, for the class of faults where the instability is due to the disturbance of the inter-area dynamics, only the slow energy is relevant and should be used to determine the critical clearing time. The above concept has been rigorously justified via the concept of integral manifolds, where corrections to the zeroth-order slow energy function are successively made 2°. Similarly, when a disturbance is such that the fast subsystem severely impacts the slow subsystem (e.g. a fault near a generator), we can introduce a coordinate transformation which reflects the effect of the boundary layer correction term on the slow subsystem as in the case of the linear dynamic system 21. Thus for a large class of disturbances the slow subsystem can be used to predict critical clearing times by using any of the techniques in direct methods. Numerical results are contained in References 17 and 20. It is conjectured that the kinetic and potential energy corrections to the energy function found in the literature 22'23 are equivalent to applying this methodology. Research in testing the system on large power networks is continuing.
Vl. Conclusions The results described in this paper illustrate the role of singular perturbations as an analytical tool in power system modelling at both the device and system level. They also illustrate how additional insights into power system dynamics can be obtained in the process. Research is continuing on improved reduced-order models of machines and excitation systems, more accurate aggregate models of groups of machines, and the use of reduced models for transient stability analysis. Additional
125
effort is needed to develop procedures and computer programs ready for use by the power industry.
11
VII. A c k n o w l e d g e m e n t s
12
This work is supported in part by the US Department of Energy under contract DE-AC0!-84CE76249 and in part by the National Science Foundation under Grant ECS-84-14677.
VIII. References 1
2
3 4 5 6 7 8
9 10
126
Concordia, C and Schulz, R P "Appropriate component representation for the simulation of power system dynamics' Syrup. on Adequacy and Philosophy of Modeling." Dynamic System Performance, IEEE Pub/. 75CHO970-4-PWR (1 975) Cate, E G, H e m m a p l a r d h , K, M a n k e , J W and Gelopulos, D P 'Time frame notion and time response of the models in transient, mid-term and long-term stability programs" /EEE Trans. PowerAppar. & Syst. Vol PAS-103 (1984) pp 143-151 Cresap, R Land Hauer, J F'Emergenceofanewswingmode in the Western Power System' IEEE Trans. Power Appar. & Syst. Vol PAS-100 (1981 ) pp 2037 2045 K o k o t o v i c , PV, O ' M a l l e y , R EJrand Sannuti, P'Singular perturbations and order reduction in control theory--an overview" Automatica Vol 12 (1 976) pp 123-1 32 Saksena, V K, O'Reilly, J and K o k o t o v i c , P V 'Singular perturbations and time-scale methods in control theory: survey 1976-1983" Automatica Vol 20 (1984) pp 273-293 C h o w , J H (ed.) Time-Scale Modeling of Dynamic Networks with Applications to Power Systems," Lecture Notes in Control and Information Sciences Vol 46, Springer, New York (1 982) Fitzgerald, A E and Kingsley, C Electric Machinery 2nd edn, McGraw-Hill, New York (1961) IEEE Committee Report Excitation system models for power system stability studies" IEEE Trans. Power Appar. & Syst. Vol 100 (1981 ) pp 494-509 Carr, J Applications of Center Manifold Theory Springer, New York (1981) Anderson, P i and Fouad, A A Power System Control and Stability Iowa State University Press, Ames (1 977)
13 14 15
16 17
18 19
20 21 22
23
K o k o t o v i c , P V and Sauer, P W 'Integral manifold as a tool for reduced order modeling of nonlinear systems: a case study' Proc. 1986 IEEE Int. Symp. on Circuits and Systems (1986) pp 908 911 K o k o t o v i c , P V, A v r a m o v i c , B, C h o w , J H and W i n k e l m a n , J R 'Coherency-based decomposition and aggregation' Automatica Vol 18 (1 982) pp 47-56 W i n k e l m a n , J R et al. 'An analysis of inter-area dynamics of multi-machine systems' IEEE Trans. Power Appar. & Syst. Vol PAS-100 (t 981 ) pp 754-763 C h o w , J H and K o k o t o v i c , P V 'Time scale modeling of sparse dynamic networks' IEEE Trans. Automatic Control Vol AC-30 (1985) pp714-722 C h o w , J H, Cullum, J and Willoughby, R A "A sparsity-based technique for identifying slow-coherent areas in large power systems' IEEE Trans. Power Appar. & Syst. Vol PAS-103 (1984) pp463-473 'Development of dynamic equivalents for transient stability studies' EPRI Report EL-456 (1 977) Pal, i A, O t h m a n , H, C h o w , J H and W i n k e l m a n , J R 'Time scale energies in power systems' Proc. 4th IFAC/IFORS Syrup. on Large Scale Systems Theory and Applications, Zurich (26 29 August 1986) A t h a y , T, Podmore, R and V i r m a n i , S A practical method for direct analysis of transient stability' IEEE Trans. Power Appar. & Syst. Vol PAS-102 (1979) pp 573-594 Ribbens-Pavella, M and Evans, F J 'Direct methods in the study of the dynamics of large-scale electric power systems--a survey' Proc. 8th IFAC World Cong., Kyoto (24-28 August 1981 ) Pal, a . A. et al. 'Decoupled stability analysis of large scale power systems using integral manifold approach' Proc. 25th IEEE Conf on Decision and Control, Athens (1 986) Chang, K W 'A diagonalization technique for singular perturbation problems' Proc. 12th AIlerton Conf. on Circuit and System Theory, University of Illinois, Urbana (1 974) Fouad, A A and V i t t a l , A 'Power system response to a large disturbance: energy associated with system separation' IEEE Trans. Power Appar. & Syst. Vol PAS-102 (1 983) pp 3534-3540 Y i - x i n Ni, Bun Kee Wee and Fouad, A A 'Accuracy of transient energy functions' Proc. 1986 North American Power Symp., Ithaca, New York (1 3-14 October 1986)
Electrical Power & Energy Systems