Modal-based stability analysis of

power systems using energy

functions
K Khorasani, M A Pal and P W Sauer
Dept. of Electrical and Computer Engineering, University
of Illinois, Urbana, I L 61801, USA

approximate method of calculating Ecr by computing E~nax

In this paper a physically based decomposition technique is on the faulted trajectory. The method of decomposition
exploited to perform direct stability analysis o f power employs the slow coherency based grouping technique 11-13
systems using the energy function method. The slow- based on a linearized analysis of the system. This grouping
coherency analysis decomposes the system into r areas is used to develop (a) slow energies associated with the
associated with (r -- 1) slow modes o f the linearized system. centres of inertia of the areas, (b) fast energies associated
The centres o f inertia (COIs) o f these areas form the slow with intermachine oscillations within an area, and (c)
subsystem and the rest o f the fast modes are associated slow-fast energies associated with both the slow and the
with different areas. The system transient energy function fast variables. It is shown that depending on the fault
E is decomposed into Eslow associated with the slow location, only a few areas are disturbed and it is possible to
variables o f the slow subsystem, Efast associated with the approximate the system energy function by the fast energies
fast variables o f the r-areas and Efastoslowwhich is a func- of the disturbed areas, global slow energy and the slow-fast
tion o f both fast and slow variables. Depending on the fault energy. The method is illustrated using a 10-machine 39-bus
location and strength o f connection between areas, the sum system.
o f Eslow , Efast_slow and Efast o f the disturbed areas consti-
tutes a good approximation to the system energy function.
Using this partial E-function and the potential energy II. Modal energy decomposition in linear systems
boundary surface (PEBS) method, both Ecr and tcr are Consider a linearized multimachine power system based on
computed. Numerical results on a l O-machine 39-bus the classical model representation. The mathematical model
system are presented in support o f the technique. is given by

Keywords: network analysis, transient stability, power = Ax (1)

system element modelling
where

I. Introduction A= M:qK ,, , x=
Since 1979 there has been renewed interest in indirect
methods of stability analysis of power systems because of M = inertia matrix and K = synchronizing coefficient
the close agreement of the critical clearing times using the matrix. The slow-coherency method exploits the eigen-
transient energy function with those obtained using step- value spread of M-tK to divide the system into areas. Let
by-step methods 1~. The transient energy function in the
kl, ~,2, . . . , kr be the slow eigenvalues and ~,r+l . . . . . ~'n the
absence of transfer conductances is a Lyapunov function, fast ones. The condition to be satisfied is that
but is not proven to be so otherwise 1'6. Hence we use the
symbol E instead of V. The success is mainly due to the use
of fault-dependent computation of the critical energy Ecr
using either the relevant UEP method or the potential I XX-~-r-~
~ =e
energy boundary surface (PEBS) method. Stability analysis
of large-scale power systems by decomposition 7-1° met with be a small number < 1. The eigenvalues of A are the square
limited success both because of the purely mathematical roots of the eigenvalues of M" 1K, the latter being real and
decomposition technique employed and the conservative ~< 0. Hence the eigenvalues of A lie on the imaginary axis in
nature of the critical energy computation. In this paper a the s-plane and represent modes of oscillation of the
new method of decomposition is used and is combined with systemll, 12
the PEBS method to compute the critical clearing time
quickly and with good accuracy. The PEBS method is an The technique discussed in References 12 and 13 decomposes
the system into r areas whose centre of inertia (COl) angles
and speeds, representing the slow variables, constitute the
Received: 26 June 1984, Revised: 2 January 1985 slow subsystem variables. The differences variables of rotor

Vol 8 No 1 January 1986 0142-0615/86/010011-06 © 1986 Butterworth & Co (Publishers) Ltd 11

angles and speeds in each area are the fast variables. If e is
small enough, there is a good decoupling o f the slow and
the fast subsystems, which is an advantage both in simula-
tion and energy computation. For energy function analysis
this is preferable though not essential. In what follows, we
will assume that an area decomposition has been performed
based on the pre-fault linearized A matrix. Moreover, to
facilitate the use of the transient energy function, we define
the fast variables o f each area to be the rotor angles and
speed referenced to the local COI. It is important to note
that the linearized model is used only for the decomposition
into areas.

III. Modal-based decomposition of nonlinear

energies
Consider the post-fault multimachine swing equations in
the synchronous frame o f reference:

Migi = P m i - - E ~ G i i - ~, (Ci] sin6q + D o cos6q) (2)

]=1

where M i = inertia constant

Ci! = EiEIBq; Di] = E i E i G q , etc. Figure 1. Area representation and interconnections

Letting Pmi - - E ~ G u = Pi, the transient energy function 1'3 is

given by
cally the connections between machines. For subsequent
n discussion we denote two typical areas a and ~ whose COI
angles are y~ and y0 respectively. Define the rotor angles
e ( 5 , oa) " = 21i =E1 M,,,4 2
with respect to the local COIs as E i = 6 i --Yc, i E a.
a = 1, 2, . . . , r and those of local COIs with respect to
1 n-1 n
the global COl 6o a s i a = ya -- 6o, a ~- 1. . . . . r. It is also
E
M T i=1 ] = / + 1
-- (6;j --
easy to verify that 0 i = x i + Yc~, i E a. Also we define

n-I n
- Ei=1 Y Cq(cosaq-cos ) MTe = ~ Mi
i@u
j=i+l
With the above definitions we now express the transient
n-1 n 6i+ 26o81--
q
energy function (3) in terms o f the slow variables y and
+ E E
i=1 j = i + l s
jS fast variables ~. It is noted that
6i+6I - 25os
(3) 6i--6]=Xi--E j
where
when i a n d j belong to the same area and that
1 n n
80 = M'-TT i=IEMi6i, MT = i=,EMi (4) 6i --61 = ~i --Ya - - E l +Y0 = Ei] - - Y ~

Og0=~ 0 when i a n d j belong to different areas. We will use the

symbol Z7 to imply summation over all machines in area
Defining for convenience COl referenced angles Oi = 6i--8o, a. The transient energy function (3) in terms of the variables
the path dependent term is numerically evaluated 1 and is x i , y c ~ , i @ a , a = 1. . . . . r is given by
approximated by
E - - E KE -F E PE (6)
I ( k =) i ( k - 1)
where
Jr" ~Dij[Cos(O~ k) -- O}k)) "4- COS(0~k- 1) __ 0}k--1))]
1 r +L r ~-2
x + o(/%, - (o(.k-l), + (5)
= i 2 ~--1MT~Y~
Let there be r areas so that we have ( r - - 1) slow modes.
Consider Figure 1, indicating the r areas. Each area has one
or more machines in it and with the internal node model = E, E~ast,a
KE + f KE ~
E~low,
~=1 e~=l
representation there are lines from each node to every
other node. The thick lines between areas show symboli- KE KE
= E/ast + Eslow (7)

12 Electrical Power and Energy Systems

 PE [ 1 '~
EPE = I - - E (eiM~-e~Mi) Efast'a= [---MTT E E
-1 i fl=a+l J i j=i+l
Ot OL
1 +Z Z
i /=i+1

(P~M! --PjMi) ('xO __~S

xii) -E a Ea S , . k) ]
× E l 1=i+1 d
1=i+1
r
+ E E (eiMi--PIMt)(xq--Y'~ (13)
+E E E (PiM] -- PjMi) (xii -- x~') ) ~¢a I
i 3=a+l j

a a & 7
The first term in the square brackets in equation (13) is
+ the internal PE and the second term is the PE along the cut-
E E E E
i ]=i+1 i ]=i+1 set around the area. Note the summation Z~ e a in equation
(13) instead of E#=a +1 as in the terms in equation (8). It
+ rE__
1 i 3=a+l ]
is not helpful to associate Efast_slow with each area since in
general the slow variables of all areas get disturbed and all
the terms in equation (11) should be computed. However,
in a very large-scale system weak coupling between some of
-- cos(if/} +yaS~)) + Z E the areas could be exploited to reduce the computation as
i /=i+I will be explained in the next section.

The slow and fast subsystem models are given in the

+ E E E IFS,0 Appendix.
i 3=~+1 ] _1 (8)
PE PE PE
= Eslow + Efast + gfast, slow IV. Approximate transient energy function Epart
As in the case of the relevant UEP concept t' 2,6, we postulate
where the concept of relevant disturbed areas (RDAs) depending
on the fault location. This in turn enables us to approximate
i(k)
i ,7 + ~.nl r.os.v. (k) + cos~i~k-O) the system energy function E by Epart. The following
× (:~i(k) + :~j(/O _ ~i(k - 1) _ ~i(k - 1)) (9) assumptions result in an overall computational saving both
in terms of integrating the faulted equations and computing
/(& . --'FS,
_ .0, -,) , -oo
q + Dq(cosx q + the transient energy function.

× O~(k)__fi(k-1)) (10) (a) If the fault is inside an area, the fast modes of that area
are excited while the fast modes of other areas are not.
](k),i ] _-- i(k-1)
"FS,i] q-
12DijE This justifies approximating the fast energies in equa-
+ I tions (7) and (12) to those corresponding to the
2nqE cos [~i(k) + .~(k)
_ (~](k) + fi(k))] disturbed area only.
+ COS [ ~ i ( k - - 1 ) + ~ ( k - 1 ) _ _ ( . ~ j ( k - 1 ) + ~ ( k - - 1 ) ) ]
(b) From singular perturbation theory 12, the fast variables
× + + + y(]) in the slow subsystem are at their post-fault steady-state
--(Xi(k-1) + y(k-i) + x/(k-1) + y(k-i))] (11) equilibrium values. Similarly, the slow variables in the
fast subsystems are at their initial values. Hence we
obtain the decoupled equations (24) and (25) of the
Although ~f + Ya = Oi, the expressions (9)-(1 l) are written
Appendix, which are approximations to the exact
in terms of ~i,.~ a because the differential equations of the
equations but faster to compute for the faulted
decoupled fast and slow subsystems are expressed in terms
trajectory.
o f ~ i and Ya respectively (see equations (22) and (23) in
the Appendix). The computation o f ~ from ya is trivial. (c) Because of (a) and (b) the computation of Efast_slow
PE
becomes simplified because it is a function of the
The decomposition of KE in terms of fast and slow com- global slow variables and fast variables of the disturbed
ponents is clearly exhibited in equation (7). In the case of area only.
PE in equation (8), the expressions in the three outer
brackets on the right-hand side of equation (8) are respec- We now define the approximate transient energy function
tively in terms of slow, fast and mixed fast-slow variables. Epart to be
We can therefore denote E vE as
KE ~KE . ~PE
E PE __ PE PE PE Epart = gslow + E ~fast,¢~ 7-JZslow
-- Eslow + Efast + Efast, slow (12)
The potential energy terms due to the interconnections + ~ PE r~ -[-
Efast, PE
Efast_slow
(14)
between areas cannot be associated with any particular aER
areas. However, if we consider only one area then it is
= E KE VE
possible to define Efast,a
PE for that area as part + Epart

Vol 8 No 1 January 1986 13

where R = set of relevant disturbed areas.
I a9
The notion of approximating the system for a specific 24 c~)38
disturbance by means of a two-machine system is con- IB 17 27
sidered in Reference 14. The notion of RDAs is, however,
more general. For a fault on the boundary between areas,
the connected areas are the RDAs for the purpose of com-
3g
I
puting Epart •

In view of the assumptions made above, the theory applic-

able to the system energy function is applicable to Epa~
also. In particular, the potential energy boundary surface
method 1'4 can be readily implemented.

Algorithm Figure 2. The New England system line diagram

1 Assume that the areas and the set R of disturbed areas
for a given fault are known. It is best to store this infor-
mation in the database from offline studies.

Integrate the approximate differential equations (24) Area I __ AreaE

and (25) of the RDAs and the slow subsystem in the
faulted state, assuming that in the non-disturbed areas
the fast variables are at their steady-state value. ( r--r--g, 4 )
U .stag r£part, r£part
PE(max)
along the faulted trajectory is
computed. This is a good approximation for Ecr.
-,4t- s s ~

4 Using this Ecr, find when Epart = Ecr to yield tcr.

V. Numericalresults ~ A r e a
The 10-machine 39-bus New England system shown in
Figure 2 is considered. The eigenvatues, in rad/s, are 0,
Area 3
+j 3.916, +j 5.95, +j 6.344, -+j 7.037, +j 7.845, +j 8.05,
+--j 8.918, +-j 9.618, + j 9.722. The zero eigenvalue is a
repeated one. Although there is no clear-cut separation of
slow and fast modes, the first five pairs of eigenvalues were Figure 3. G r o u p i n g and eutsets f o r various faults on the
taken as slow and the rest as fast. system. - - - Fault bus 29 - disturbed areas 2 and 4;
+ + + Fault bus 31 - disturbed areas 1 and 2; x - - - x
We thus have a 5-area grouping with the largest machine Fault line 1 5 - 1 6 ' - disturbed area 2
associated with the zero eigenvalue. The slow coherency
based grouping TM13 resulted in the following areas:

Area 1 Machines 2, 3 Table 1

Area 2 Machines l, 4, 6, 7, 8
Results using Results using
Area 3 Machine 5
Area 4 Machine 9 Fault Vpart Vsyst
Area 5 Machine 10 location Disturbed
and area areas Vcr ter (s) Ver tcr (s)
This grouping is shown in a simple manner in Figure 3.
Bus29 (4) 2,4 2.71 0.12 2.819 0.124
Some of the connections between machines are weak
Bus31 (1) 1,2 5.75 0.21 7.03 0.22
while others are not. The areas are therefore semi-
15-16"(2) 2 6.23 0.14 6.57 0.142
4"-14(2) 2 10.23 0.228 10.74 0.231
Four faults were considered - two on the boundary between
areas and the other two well inside an area. These faults
are:

1 Fault on bus 29 (no line switching) disturb only area 2. Using the algorithm, Eer and tcr are
2 Fault on bus 31 (no line switching) computed. They agree well with the values obtained by using
3 Fault 15-16" followed by line switching the system energy function.
4 Fault 4 " - 1 4 followed by line switching
(* indicates the point of application of the three-phase fault) V.1 Discussion
For faults inside an area on the transmission system, it is
Fault at bus 29 disturbs both areas 2 and 4 while fault at sufficient to compute only EFast, u of that area in addition
bus 31 disturbs areas 1 and 2. Faults 15-16" and 4 " - 1 4 t o E s and EF, S. For faults on the boundary, at most 2 or 3

14 Electrical PowerandEnergySystems
areas will get disturbed. There is a considerable saving in Scale Power Systems, New Delhi, India (16-18 August 1979)
terms of less DEs to be integrated and a smaller E-function Pergamon Press, UK
4 Kakimoto, N, Ohsawa, Y and Hayashi, M 'Transient stability
to be evaluated. Once the area grouping is assumed and the analysis of electric power system via Lur~ type Lyapunov
disturbed areas are known in the case of a large-scale power function, Parts I and II' Trans. lEE Japan Vol 98 No 5/6 (May/
system, the technique is very efficient and time saving. June 1978)
Since the values of tcr using Esystem agree well with the 5 Fouad, A A and Stanton, S E 'Transient stability of multi-
machine power system, Parts I and I1' IEEE Trans. Power
simulation, the latter are not listed separately. The concept
Appar. & Syst. Vol PAS-1 O0 (July 1981) pp 3408-3424
of cutsets for different faults in Figure 3 is useful in the 6 Pal, M A Power System Stability North Holland, New York
following sense. The PE terms in equation (14) need only (1981)
be computed along the cutset and internal to the disturbed 7 Pal, MA and Narayana, C L 'Stability of large scale power
areas. The slow KE terms must be evaluated globally while systems' Prec. 6th IFAC World Congress, Boston (1975)
8 Jocic, Lj B, Ribbens-Pavella, M and Siljak, D D 'Multi-machine
the fast KE are computed only for the disturbed areas. power systems: stability, decomposition and aggregation'
IEEE Trans. Automatic Control Vol AC-23 (1978) pp 325-332
9 Araki, M, Metwally, M M and Siljak, D D 'Generalized de-
compositions for transient stability stability analysis of multi-
VI. Conclusion machine power systems' Prec. JACC, San Francisco (August
1980)
The application of modal-based stability analysis by the
10 Pal, M A and Vittal, V 'Multi-machine stability analysis using
direct method has been demonstrated. It involves the vector Lyapunov functions with inertial center decomposition'
following steps: Int. J. Elec. Power & Energy Syst. Vol 5 No 3 (July 1983)
11 Anderson, P M and Fouad, A A Power System Control and
1 Definition of areas. The modal based slow-coherency Stability Iowa State Univ. Press, Ames, Iowa (1977)
12 Chow, J H Time Scale Modelling of Dynamic Networks with
grouping is used. Applications to Power Systems Springer-Verlag, New York
(1982)
For each fault, we associate relevant disturbed areas 13 Winkelman, J R, Chow, J H, Wheeler, B C, Avramovic, B and
denoted by the set R which may consist of a single area Kokotovic, P V 'An analysis of inter-area dynamics of multi-
only or more than one area depending on fault location. machine systems' IEEE Trans. Power Appar. & Syst. Vol PAS-
100 (1981) pp 754-763
For a previously specified list of contingencies, this 14 Fouad, A A and Vittal, V 'Power system response to a large
information can be stored in the database. disturbance: energy associated with system separation' IEEE
Trans. Power Appar. & Syst. Vol PAS-102 (November 1983)
For a given fault, the system transient energy function 15 Zaborsky, J, Whang, K W, Huang, G M, Chiang, L J and Lin, S Y
'A clustered dynamic model for a class of linear autonomous
is approximated as systems using simple enumerative sorting' IEEE Trans. Circuits
& Syst. Vol CAS-29 No 11 (November 1982)
Ep~= E pKE
~+ E pPE
~ 16 Michel, A N, Fouad, A A and Vittal, V 'Power systems transi-
ent stability using individual machine energy functions' IEEE
Trans. Circuits & Syst. Vol CAS-30 No 5 (May 1983)
4 Using the PEBS method, Ecr and ter are computed with 17 Bergen, A R and Hill, D J 'A structure preserving model for
Epart function. power system stability analysis' IEEE Trans. Power Appar. &
Syst. Vol PAS-IO0 (January 1981) pp 25-35
The decentralized nature of computation of Ecr and tcr 18 Athay, T and Sun, D I 'An improved energy function for transi-
depending on the fault location makes the method quite ent stability analysis' Prec. /nt. Syrup. Circuits and Systems,
Chicago (April 1981 )
attractive for large-scale power system applications. New
areas of research are quite evident from this work. It would
be interesting to establish connection between this technique
Appendix
The exact and approximate differential equations for the
and the more specialized case of individual machine energy
slow and fast subsystems for the nonlinear model of equa-
function 16. Extension to sparse formulation is another
tion (2) are presented in this Appendix.
area of future investigation 17' as. Extensions to larger
systems are in progress.
Let U be an n x r partition matrix whose rows represent the
n machines and the columns represent the areas 1 , . . . , r.
Let the machines be ordered by area. The i--a entry in U is
VII. Acknowledgements 1 if machine i belongs to area a, and 0 otherwise. Hence
The authors wish to thank Prof. P Kokotovic of University each row contains only one non-zero entry which is equal
of Illinois and Bill Price of General Electric Co. Schenectady to 1. Let (2 be an r × n aggregation matrix 12
for useful discussions. This research was supported in part
by the General Electric Co. and by the Power Affiliates = [ U ' M U]-aU' M (15)
Program of the University of Illinois.
The slow and fast variables are defined as

VIII. References y = C6 (16)

1 Athay, T, Sherket, R, Podmore, R, Virmani, S and Puech, C
'Transient energy stability analysis' Prec. Conf. System Engin-
eering for Power, Dares, Switzerland Also US Department of
= (I--UC) 6 (17)
Energy Publ. No Conf-790904-PL 1980
2 Athay, T, Podmore, R and Virmani, S 'A practical method for 6 =Uy+K
direct analysis of transient stability' IEEE Trans. Power Appar.
& Syst. Vol PAS-98 No 2 (March/April 1979) pp 573-584 where I = identity matrix. With the above definitions
3 Ribbens-Pavella, M, Gruijc, Lj T, Sabatel, J and Bouffioux, A
'Direct methods for stability analysis of large-scale power Y = [Yl...Y~...Yr]' and ~ = [ ~ ' . . . ~ . . . ~ r ] ' where ~
systems' Prec. IFAC Symp. Computer Applications in Large is the state vector corresponding to fast variables of area ct

Vol 8 No 1 January 1 9 8 6 15
 with.~i = 6i --Ya, i E a, a = 1, 2 . . . . . r. We define the With these notations the swing equations (2) become
following n x n matrices:
= -- [I --U(2] M -1 [[sin6] I + [sin ~5]E + [cos~5] 1
[sin 6 ] = Cii sin 6q i ~j
(18) + [cos6]El En+ [ I - - U C ] M-zPF.n (22)
=0 i=j
3 ; = - - C M - 1 [ [ s i n 3 ] E + [cos3] E] En + C.M-1PF.n (23)
[cos6]=Dijcos6 U i4=j
(19) since CM -1 [sin 3] I = CM -1 [cos6] I = 0.
=0 i=j
Functionally [ ]l are functions of.~ only and [ ]E are func-
Also we define tions of both x and y. Following the thinking in the linear
case 11, the decoupled state-space model assumes that (a) the
M = diag(M/) i = 1, 2 . . . . , n fast variables in equation (23) are at their post-fault steady-
state equilibrium value ~s and in equation (22) we neglect
P = diag(Pi) i = 1, 2 . . . . . n the terms [ ]E involving external connections. This results
in the following r subsystems and one slow subsystem all
En = n-vector with all entries equal to 1 decoupled from each other:

Furthermore we decompose [sin 6 ] and [cos 6 ] respectively .~ = -- [I -- UC] M -1 [[sin 6 ]I + [cos ~]I En
as the sum of two matrices, one representing internal con- + [ I - UC] M-IPEn (24)
nections of an area and the other pertaining to external
connections as ~; = --C-M-1 [[sinS] E + [C°s6]E] F'n + CM-1PF-n (25)

[sin6] = sin6] I + [sin6] E (20) Equations (24) and (25) are used for integrating the faulted
equations in a fast manner. Note that only the disturbed
[cos~ l = [cos~] I + [cos~] E (21) area equations in equation (24) (fast subsystem) need to be
integrated. The fast variables of the non-disturbed areas are
The [ ]I is block diagonal and [ ]E has 0 entries along the held constant as explained in Section IV. The 6ij terms are
diagonal and non-zero entries otherwise. expressed in terms of ~ and y as explained in Section III.

16 Electrical Power and Energy Systems