RETURNING MATERIALS:
IVIESI_J Place in book drop to ,
remove this checkout from E
4::::;:Eii_ your record. FINES will "
_ be charged if book is
returned after the date ;
stamped be10w. \\

 

 

 

 

 

LYAPUNOV STABILITY BASED NETWORK CONDITIONS.
FOR CHARACTERIZING RETENTION AND
LOSS OF POWER SYSTEM TRANSIENT STABILITY

BY

Gholamhossein Sigari

A DISSERTATION

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

DOCTOR OF PHILOSOPHY

Department of Electrical Engineering
and Systems Science

1985

 

 

Copyright by
Gholamhossein Sigari

1985

ABSTRACT
LYAPUNOV STABILITY BASED NETWORK CONDITIONS

FOR CHARACTERIZING RETENTION AND
LOSS OF POWER SYSTEM TRANSIENT STABILITY

BY

Gholamhossein Sigari

A topological energy function has been derived from the
first principles that retain the network without aggregation
back to generator internal buses and includes a general
description of real and reactive power load models as a
function of voltage. For the special case of a constant
real power, constant current reactive load model and the one
axis generator model a topological Lyapunov energy function
has been derived.

Based on the Popov stability criterion a characteriza-
tion of the region of stability and a definition of loss of
transient stability have been presented. Then three differ-
ent theorems stating necessary and sufficient conditions for
the loss of transient stability have been established. The
conditions stated in these theorems describe an estimate of
the region of instability. Although a precise boundary
between the region of stability and the region of
instability is not established a relationship between the

regions of stability and instability has been established.

..;

sh

A
-..I

V."
O..-

(H

u’ )

,‘find'

iv “

H‘AL

-
5..

Q .5.
Ui~°i

(0
GD

Gholamhossein Sigari

The necessary condition for instability provides a condition
that if satisfied, will ensure that a loss of stability will
occur. Thus the necessary and sufficient conditions provide
easily tested conditions on the entire states of the system
which can be clearly identified with the loss of synchronism
inherent in the loss of transient stability. These
necessary and sufficient conditions for loss of transient
stability permit the proper definition of stability margin
that measures the relative security of the system for a
particular fault cleared at a particular clearing time with
a given fault clearing action.

Stability criteria are proposed based on the Popov
stability criterion and these necessary and sufficient
conditions for loss of stability. These stability criteria
are then related to the PEBS method developed for the single
machine energy function and the PEBS method based on the
outset energy function.

The above theoretical results were tested on the 39 bus
New England system. These tests confirmed the properties of
the energy function and Lyapunov energy functions. The
stability criteria tested clearly and accurately determined
the critical outset and critical clearing time and verified
the accuracy and validity on the stability criteria that

were developed based on the theory developed.

To Mojtaba and Horiyeh

iii

aivis
ccurs

flfl‘uv
4&5!

E u
.x “we:

of 31's
SIaRR‘
.. w.

i5 315

I

U
ire-t:
ta,
bvie

ACKNOWLEDGEMENTS

 

I would like to express my appreciation to my thesis
advisor, Dr. R.A. Schlueter, for his guidance throughout the
course of my doctoral program. Dr. Schlueter served not
only as advisor but also a friend and his friendship will
never be forgotten.

Special acknowledgement is made for the contributions
of Drs. G. Park and H. Khalil. The participation of Drs. M.
Shanblatt and D. Yen and N. Guven on my guidance committee
is also appreciated and acknowledged.

Dr. Schlueter and the Division of Engineering research
are thanked for providing financial support over the final
two years of my doctoral program.

Special thanks to ”Sam" Miller with the Wizards of
Words for her word processing talents.

Most of all, I would like to thank my father, for his
encouragement and financial support during the course of my

graduate studies.

iv

TABLE OF CONTENTS

Page
LIST OF FIGURES 000.000.000.000...00000000000000000000o Vii

CHAPTERl-INTRODUCTION oooooooooooooooooooooooooooooo 1

CHAPTER 2 ‘ REVIEW OF LYAPUNOV STABILITY
THEORY APPLIED TO POWER SYSTEMS ........... 7

2.1 Equal Area Criterion for the Single Machine
Infinite Bus "Odel .OOOOOOOOOOOOOOOOOOOOOOO0000. 7

2.2 Potential Energy Boundary Surface, (PEBS)
MethOd 0.0.000.........OOOOOOOOOOOOOO0.0.0.00... 13

2.3 Review of Literature on Energy and Lyapunov
Functions OOOOOOOOOOOOOOOOOOOOOOOO0.00.00.00.00. 16

2.4 Lyapunov's Direct Method and Present Energy
Functions ......OOOOOOOOOOOOOOOOOO00.000.00.000. 20

2.4.1 Lyapunov Energy Function for the
Aggregated Network Model ................. 21

2.4.2 Individual Machine Energy Function ....... 24

2.4.3 A TOpological Lyapunov Energy Function ... 26

2.5 Summary ........................................ 28
CHAPTER 3 - TOPOLOGICAL ENERGY FUNCTION ............... 30
3.1 Derivation of the Energy Function .............. 32

3.2 Proof of Conservativeness of the Energy
Function ......OOOOOOOOOOOOOOO..OOOOOOOOOOOOOOOO 36

3.3 Energy Function with Different Load Models ..... 40
3.3.1 A General Real Power Load Model .......... 40

3.3.2 The Energy Function With Different
Reactive Load "Odels ......OOOOOOOOOOOOOO. 41

a"
q
\n.

..n

m...

[J-

(’1

LI.

‘1

3.4Simu1ationResu1ts ......OOOOOOOOOOOOOOOOO......

Appendix 30A .........OOOOOOOOOOOOOOOO......OOOOOOOO

CHAPTER 4 - A LYAPUNOV ENERGY FUNCTION ................

4.1 Theoretical Background .........................

4.2 Lyapunov Energy Function Derivation ............

4.2.1 System

Medeling OOOOOOOOOOOOOOOOO.........

4.2.2 Satisfaction of Moore-Anderson Theorem
conditions 0.00.0000...OOOOOOOOOOOOOOOOOOO

4.3 Construction

of the Lyapunov Function ..........

4.4Simu1ationResu1ts ......OOOOOOOOOOOOO0.0.0.0...

Appendix 40A .....OOOOOOOOOOOOOO.......OOOOOOOOOOOO.

Appendix ‘08 ....OOOOOOOOO0.000000000000000000000000

Appendix 40C ..OOOOOOOOOOOOOOOOOOO......OOOOOOOOOOOO

CHAPTER 5 - STABILITY CRITERIA BASED ON REGION
OF STABILITY AND INSTABILITY ..............

5.1 Introduction

......IOOOOOOOOOOOO0.00.0.0...CO...

5.2 Regions of Stability and Instability for
Transient Stability Models .....................

5.3 A Necessary and Sufficient Condition for
Loss of Transient Stability ....................

5.4 Stability Criteria Based on the Boundaries
of the Region of Instability and Stability .....

CHAPTER 6 - VERIFICATION OF STABILITY
CRITERIA VIA SIMULATION ......OOOOOOOOOOOO.

CHAPTER 7 - OVERVIEW AND FUTURE RESEARCH o o o o o o o o o o o o o o

7.10verView 00...........OOOOOOOOOOOOOOOO0....0....

7.2 Future ResearCh ......OOOOOOO......OOOOOOOOOO...

LIST OF REFERENCES

vi

46
86
88
89
92

93

95
101
116
134
135
140

142
142

144

150

165

169
189
189
192

194

Fig.

Fig.

Fig.

Fig.

Fig.

Fig.

Fig.

Fig.

Fig.

Fig.

(1.a)

(1.b)

(2.a)

(2.b)

(3)

(4)
(5.a)

(5.b)

(5.c)

(5.d)

LIST OF FIGURES

Power-angle curves showing the
critical clearing angles 6 and
c
areas A1 and A3 .....................OOOCCO

Equal area criterion net acceleration
energy E(t) versus time ...................

Maximum values of generator angles
for different clearing times ..............

Potential energy stored in the

system for different clearing times .......

Power system model. For the
classical model all the elements
inside the dotted box are aggregated ......

New England study system ..................

Variations of angles of generators

7, 8 and 9 for stable case where

the fault on line 26-27 is cleared

in 0.35 seconds by removing the

faulted line ..............................

Variations of potential. kinetic
and total energy for the stable case
where the fault is on line 26-27 ..........

Variations of real load energy and
transmission line energy and energy

due to angle displacement of

generators rotors for the stable case

when the fault is on line 26-27 ..........

Variation of reactive load energy for
the stable case where fault is on
line 21:22 ......OOOOOOOOOO00.00.000.000...

vii

11

14

15

22

47

54

55

56

57

RI.

HI.

T i-

RIL

Fig.

Fig.

Fig.

Fig.

Fig.

Fig.

Fig.

Fig.

Fig.

(6.a)

(6.b)

(6.c)

(6.d)

(7.3)

(7.b)

(7.c)

(7.d)

(8.a)

Variation of angles of generators
7, 8 and 9 for unstable case where
the fault on line 21-22 is cleared

in 0.36
removed

Changes

seconds and line 26-27 is

of potential, kinetic and

potential energy for the unstable
case when the fault is on line 26-27 ......

Variation of real load and trans-

mission line energies and energy due

to angle displacement of generator

rotors for the unstable case when

the fault is on line 26-27 ................

Changes of reactive load energy for
the unstable case where fault is on

line 26

-27 0..........OOOOOOOOOOOIOO0......

Variation of generators angles for
stable case that fault occurs on

line 26
at 0.30
removed
model)

-27. The fault is cleared

seconds and line 26-27 is
. (Constant impedance load

Changes of kinetic, potential and

total energies for the stable case

when the fault is on line 26-27.

(Constant impedance load model) ...........

Changes of energy due to generator

angle displacement and loads real

power energy and energy stored in
transmission grid for the stable

case where fault on line 26-27 is

cleared in 0.3 seconds. (Constant
impedance load model) .....................

Changes of loads reactive power

energy for the stable case when

fault is on line 26-27. (Constant
impedance load model) .....................

Variation of generators angles for
unstable case when the fault occurs
on line 26-27. The fault is cleared
at 0.31 seconds and line 26-27 is
removed. (Constant impedance load

model)

viii

58

59

60

61

62

63

64

65

66

‘Q
‘I.

.\Q
'A

»4
Qt.

Fig.

Fig.

Fig.

Fig.

Fig.

Fig.

Fig.

Fig.

Fig.

(8.b) Variations of kinetic, potential and
total energies for the unstable case
when the fault is on line 26-27.
(Constant impedance load model) ...........

(8.c) Variations of real load energy,
enery due to generator angle
displacement and energy stored in
transmission line for unstable
case when fault is on line 26-27.

(Constant impedance load model) ...........

(8.d) Variations of reactive load
energies for the unstable case when
the fault is on line 26-27.
(Constant impedance load model) ...........

(9.a) Variation of angles and generators
7, 8 and 6 for the stable case when
the fault on line 21-22 is cleared
at 0.15 seconds and line 21-22 is
removed. (Constant current load

medal) O0.0......OOOOOOOCOOOOOOOOOO0.00....

(9.b) Changes of potential, kinetic and
total energies for the stable case
where fault was on line 21-22.
(Constant current load model)

(9.c) Variations of real power load
energy. energy stored in trans-
mission lines and energy due to
generator angle displacement for
stable case when the fault is on

line 21-22.
load model)

(Constant current

(9.d) Changes of reactive power load
energy for the stable case when
the fault is on line 21-22 ...

(10.a) Variation of angles of generator
6, 7 and 8 for unstable case where
fault on line 21-22 is cleared in
0.36 seconds and the faulted line

is cleared

(10.b) Variations of potential, kinetic
and total energy for unstable
case when fault is on line 21-22.
(Constant current load model)

ix

67

68

69

70

71

72

73

74

75

'~
6"

H .

Ink

RI. 4

Fig. (10.c) Variations of loads real power
energy: and energy stored in
transmission lines and energy due
to generator angle displacement
for unstable case when fault is
on line 21-22. (Constant current
load model) ............................... 76

Fig. (10.d) Variations of loads reactive power
energy for the unstable case when
fault is on line 21-22. (Constant
current load model) ....................... 77

Fig. (11.a) Variation of angle of generators
6, 7 and 8 for the stable case
when fault occurs on line 21-22.
The fault is cleared at 0.13
seconds and line 21-22 is removed.
(Constant impedance load model) ........... 78

Fig. (11.b) Variations of total, kinetic and
potential energies for the stable
case when the fault is on line
21-22. (Constant impedance load
model) .................................... 79

Fig. (11.c) Changes of energy stored in
transmission line and energy due
to generator angle displacement,
and energy stored in real power
loads for the stable case when
the fault is on line 21-22.
(Constant impedance load model) ........... 80

Fig. (ll.d) Variations of reactive load power
energy for the stable case where
fault is on line 21-22. (Constant
impedance load model) ..................... 81

Fig. (12.a) Variation of angle of generators
6, 7 and 8 for unstable case when
the fault occurs on line 21-22.
The fault is cleared at 0.14
seconds and line 21-22 is removed.
(Constant impedance load model) ........... 82

Fig. (12.b) The variation of total, potential
and kinetic energy for the unstable
case when the fault is on line 21-22.
(Constant impedance load model) ........... 83

his

\u.

Fig.

Fig.

Fig.

Fig.

Fig.

Fig.

Fig.

Fig.

(12.c)

(12.d)

(l3.a)

(13.b)

(13.c)

(13.d)

(l4.a)

(14.b)

Changes of real load energy and

energy due to generator angle

displacement and energy stored in

transmission lines for the unstable

case when the fault is on line

21-22. (Constant impedance load

model) .................................... 84

Changes of reactive load power

energy for the unstable case when

the fault is on line 21-22.

(Constant impedance load model) ........... 85

Variation of generator angles for

the stable case where the fault

occurs on line 28-29, the fault

is cleared at 0.2 seconds and line

28-29 is removed .......................... 118

Variation of potential, kinetic

and total energies for the stable

case where the fault occurs on line

28-29, the fault is cleared at 0.2

seconds and line 28-29 is removed ......... 119

Comparison of different components

of potential energy for the stable

case when fault occurs on line

28-29, the fault is cleared at 0.2

seconds and line 28-29 is removed ......... 120

Variation of flux linkage and

reactive load energies for the

stable case when the fault occurs

on line 28-29, the fault is

cleared at 0.2 seconds and line

28-29 is removed .......................... 121

Variations of generators angles

for unstable case when the fault

occurs on line 28-29, the fault

is cleared at 0.21 seconds and

line 28-29 is removed ..................... 122

Variation of potential, kinetic

and total energy for the unstable

case when fault occurs on line

28-29, the fault is cleared at

0.21 seconds and line 28-29 is

removed ................................... 123

xi

....
HI.

Fig.

Fig.

Fig.

Fig.

Fig.

Fig.

Fig.

Fig.

(14.c)

(14.d)

(15.a)

(15.b)

(15.c)

(15.d)

(16.a)

(16.b)

Comparison of different components

of potential energy for the unstable

case when fault occurs on line 28-29,

the fault is cleared at 0.21 seconds

and line 28-29 is removed .................

Variation of flux linkage and

reactive load energies for the

unstable case when the fault

occurs on line 28-29, the fault

is cleared at 0.21 seconds and

line 28-29 is removed .....................

Variation of generator angles for
stable case when the fault occurs
on line 14-33, the fault is cleared
at 0.23 seconds and line 14-33 is

removed

Variation of potential, kinetic and

total energies for the stable case

when the fault occurs on line l4-33,

the fault is cleared at 0.23 seconds

and line 14-33 is removed .................

Comparison of different components

of potential energy for the stable

case when fault occurs on line 14-33,

the fault is cleared at 0.23 seconds

and line l4-33 is removed .................

Variation of flux linkage and

reactive load energies for the

stable case when the fault occurs

on line 14-33, the fault is

cleared at 0.23 seconds and line

14-33 is removed ..........................

Variation of generators angles

for the unstable case when the

fault occurs on line 28-29, the

fault is cleared at .24 seconds

and line 14-33 is removed .................

Variation of potential, kinetic

and total energies for the unstable
case when the fault occurs on

Line l4-33, the fault is cleared

at 0.24 seconds and line 14-33 is

removed

xii

124

125

126

127

128

129

130

131

I
t

 

Fig.

Fig.

Fig.

Fig.

Fig.

Fig.

Fig.

Fig.

Fig.

(l6.c) Comparison of different components
of potential energy for the unstable
case when the fault occurs on line

(16.d)

(17)
(18)

(19)

(20)

(21.a)

(21.b)

(21.c)

14-33,

the fault is cleared at 0.24
seconds and line 14-33 is removed

Variation of flux linkage and
reactive load energies for unstable
case when the fault occurs on line

14-33 '

Popov criterion for branch ij

the fault is cleared at 0.24
seconds and line 14-33 is removed

Regions of stability and instability
in branCh angle space ....OOOOOOOOOOOOOO...

Comparison of oi.f(o..) for
stable and unstaBle

fault occurs on line 26-27.

C386 when the
The

fault has been cleared and line

26-27 has been removed.
energy function)

Values of V
axis general

(t
or

t

)
fdee

(Lyapunov

for the single
1 and the

constant current reactive load

model for the fault on line 26-27

Comparison of o..f(o..) for stable
and unstable caég whed the fault

occurs on line 14-33.

The fault

has been cleared and line l4-33

has been removed.

function)

(Lyapunov energy

Comparison of o..f(c..) for

different lines

ldn seable and

unstable case when the fault

occurs on line 14-33.

The fault

is cleared and line l4-33 is

removed.

(Lyapunov energy function)

Comparison of oi.f(o..) for
different lines 3n sfable and
unstable case when the fault

occurs on line 14-33.

The fault

is cleared and line 14-33 is

removed.

(Lyapunov energy function)

xiii

132

133
152

163

177

178

179

180

181

Fig. (21.d) Comparison of oi.f(o..) for
different lines 3n szble and
unstable case when the fault
occurs on line 14-33. l4-33 is
removed after the fault is cleared.
(Lyapunov energy function) ................ 182

Fig. (21.e) Comparison of o..f(o..) for stable
and unstable casé whed the fault
occurs on line 14-33.‘ The fault
has been cleared and.line 14-33
has been removed. (Lyapunov energy
function) ................................. 183

Fig. (22) Values of V (tf,t ) for the single
axis genera or moSel and the
constant current reactive load
model for the fault on line 14-33 ......... 184

Fig. (23) Comparison of a..f(c..) for stable
and unstable ca§3 whéfe the fault
occurs on line 26-27. The fault
is cleared and line 26-27 is removed.
(Constant current load model) ............. 185

Fig. (24) Values of V (tf,t ) for the constant
voltage behind trgnsient reactance
generator model and the constant
current reactive load model ............... 186

Fig. (25) Comparison of oi.f(ci.) for different
lines on stable 3nd unstable case
when fault occurs on line 26-27.
The fault is cleared and line 26-27
is removed. (Constant impedance
load model) ............................... 187

Fig. (26) Values of V (tf,t ) for the constant
d

voltage behin trgnsient reactance
and constant impedance load model ......... 188

xiv

CHAPTER 1

 

INTRODUCTION

The present method of determining retention or loss of
stability due to transients caused by electrical faults or
loss of generation contingencies is time-step integration of
the- differential equations. The transient stability
programs that perform the time step integration can take as
long as 15 CPU minutes on the largest and fastest computers
to determine retention or loss of stability for a single
fault, cleared at a specific time with a specific line
switching action, on a specific system* with a particular
network configuration; load level, type and distribution;
unit commitment, generation dispatch, and base case load
flow. To determine the transfer capabilities from the
Northwest (Oregon, Washington) and the East (Utah, Nevada,
Arizona) to southern California for each season (summer,
winter, spring, fall) requires solutions of hundreds of
fault cases that keep a VAX 785 running continuously for
three months. Although not all utilities must be as

concerned with limiting power transfers to maintain

 

*The computation time depends on the model and the
complexity and size of the system. The CPU time reported
above is obtained from results on the WSCC system.

R.»

5‘4

 

2
transient stability over a wide range of operating
conditions and fault cases; all utilities perform transient
stability simulations and set transfer limits on generating
units and across interfaces that are based on transient
stability simulations for the operating conditions expected
over each season or year.

Transient stability simulations are performed daily on
the Ontario Hydro System in order to adjust the seasonal
transfers based on changes in network configuration, unit
commitment, load level and distribution, voltage profile,
and load flow conditions expected over the next day.

Operators at the control center desire the ability to
perform transient stability simulations to assess the
stability margin of the system for a current operating
condition that was not anticipated by operation planners in
their seasonal or daily assessments. The capability of
performing transient stability simulations on line is
completely impractical at present based on the very large
computation time required even for the fastest computers
usimg the simplest of power system models. Thus, the
operator cannot adjust transfer limits on-line based on
transient stability simulations for faults given the
operation condition presently experienced or anticipated to
occur over the next few hours.

Adaptive protection schemes that would only trip
generators when a loss of transient stability would be immi-

nent for a particular operating condition, given occurrences

3

(of a particular fault, would be feasible if transient
stability simulation couLd be performed with computational
requirements equal 1x) load flow (steady state solution of
the network equations). The tripping of large radially
connected generators can cost a utility $20,000 an hour in
increased fuel costs. The generator tripping may require
the generator boiler to be shut down and restarted, which
can take up to a week in some cases. A utility could save
fuel costs by not tripping a particular generator for every
possible contingency but only for those that would cause a
possible loss of transient stability.

The time step integration of the transient stability
model, includes

1) 2N nonlinear algebraic equations that repre-
sent the N bus transmission network;

2) the kM cdifferential equations that. describe
the synchronous generator, exciter, power
system stabilizer, and turbine energy system
for the M generating units in the model. The
order of the model k can be as low as 2 for a
classical machine model and as high ten for a
more sophisticated model that requires tremen-
dous CPU, I/O and memory.

A fast transient stability method has long been desired
that could determine retention or loss of stability with the
computation requirements of a load flow. Two developments
are necessary to make such a development a reality:

1) the determination of a Lyapunov function and a
characterization of the region of stability
for the particular system and fault case;

2) a method for predicting whether the system
fault trajectory will remain within the region
of stability and remain stable or else enter
the region of instability and thus lose
stability.

ad

.6
‘4-

4

Significant research effort has been expended on
developing the energy function and characterizing the region
of stability. Up to now, these efforts have been only
marginally successful at best. The research reported in
this thesis is directed at developing Lyapunov functions or
energy functions for much more detailed transient stability
models than have been previously developed. A charac-
terization of the region of stability and region of
instability is then established that greatly extends
previous results. This characterization of the region of
stability is then used to justify a stability criterion that
is related to the potential energy boundary surface (PEBS)
methods developed for the single machine energy function.

It should be noted that no effort will be made to
develop a method for predicting whether a fault trajectory
will remain within the region of stability or enter the
region of instability given a particular fault, fault
clearing time and clearing action, and a particular system
and its operating condition. Such a method was developed in
a previous thesis for the classical power system model and
could be extended to the more complex models. [21]

The ultimate development of a fast transient stability
assessment method can

1) greatly reduce the computational requirements

for determining transfer limits in seasonal
transmission performance assessments;

2) permit daily transfer limit adjustments based

on hundreds of transient stability simulations

at modest computational requirements and
costs. At present, only a few transient

5

stability simulations can be performed to
adjust transfer limits based (n1 the expected
operating condition over the next day;

3) permit on-line transient stability simulation
and transfer limit adjustment by operators.
At present, the capability to simulate faults
on-line is impossible;

4) permit adaptive generator tripping, that would
only trip a generator when the particular
contingency that occurred on a system with a

given operating condition would cause loss of
transient stability.

Power System Stability Models
Consideredgin This Thesis

.A review of the literature on the development of Lya-
punov functions and energy functions is given in Chapter 2
along with discussion of the equal area and PEBS methods for
determining retention or loss of stability. In Chapter 3, a
topological energy function is derived from first principles
which retains the network without aggregation back to
generator internal buses and which includes a description of
real and reactive power load models as a function of
voltage. Similar models have been hypothesized [9, 16] but
never derived from first principles.

A Lyapunov function is derived for a model that

l) retains the network and does not aggregate the
network back to internal buses;

2) models the real power load as constant power
and the reactive power load as constant
current;

3) models the flux linkage decay in the synchro-
nous machine as well as the electromechanical
component of the machine model.

6

This Lyapunov function is derived using the Popov
Criterion of the Moore Anderson [10] and the construction
method of Willems [3].

In Chapter 5, a discussion of the region of stability
based on the Popov stability criterion and the definition of
the loss of transient stability is presented. Three
different theorems stating necessary and sufficient condi-
tions for loss of transient stability are established. A
stability criteria; based on the region of stability and the
region of instability is proposed for determining the
critical cutset and critical clearing time. AA second
stability criterion is proposed that based on a performance
measure determines the critical clearing time.

The stability criteria developed in Chapter 5 are
reviewed and then applied to four fault cases. The results
indicate that both stability criteria can determine the
critical clearing time with no detectable error. These
results confirm the theory developed in Chapter 5, upon
which they were proposed. Chapter 7 reviews the research
performed in the thesis and its contribution. A discussion

of extensions of the research is also given.

transient stability model is introduced to indicate the
general form of the multimachine transient stability model
which will be presented and discussed later.

criterion is presented

CHAPTER 2

REVIEW OF LYAPUNOV STABILITY THEORY
APPLIED TO POWER SYSTEMS

The objectives of this chapter are

l)

2)

3)

4)

5)

EgalyArea Criterion for the Single Machine Infinite

to introduce and review the simplified single
machine infinite bus power system transient
stability model;

to discuss the equal area and potential
energy boundary surface methods for deter-
mining retention or loss of stability for
this single machine infinite bus model. This
discussion provides understanding and motiva-
tion for the analysis that derives the region
of stability and the region of instability
for multimachine transient stability models
in Chapter 5;

to review Lyapunov stability theory:

to review the literature in [3, 9, 17] and
developing Lyapunov based and integral based
energy functions for various power system
transient stability models;

to present the form of the energy functions
derived for different power system transient
stability models.

Bus Model

The simple single machine

7

infinite bus power system

The equal area

in order to indicate the energy

 

N

\I‘

on

8
transfers durimg a fault and during the post fault period
after the fault is cleared. The fault is cleared by an
appropriate line switching action that isolates the faulted
branch.

The single» machine infinite bus transient stability

model has the form:

ME: Pm-C SinG (W) (2.1)

where

M: generator inertia constant,(J.s)

6: generator angle, (Radians)

Pm: mechanical input, (W)

8.3:;

C: x where E is generator's internal voltage, Va
is infinite bus voltage and x is the equivalent
transmission line reactance connecting the generator
internal bus and the infinite bus.

The prefault model and the prefault stable equilibrium

point 681 are defined by

$1 (2.2)

M6 = o = Pm-Co S1n 6
where Pm' the prefault power angle curve CO Sins , and 631
are shown in Figure l.a.

The faulted system model is defined as:

M6 = Pm - C1 Sin 6 (2.3)

where the faulted power angle curve C1 Sin 6:1 is shown in

Figure l.a. The post fault system is defined as:

where the post fault power angle C2 Sin 61 is shown in
Figure l.a. The post fault stable equilibrium point 652 is

defined as:

9

Fig. (l.a) Poweraangle curves showing the critical clearing
angles 6c and areas Al and A3.

 

p
4
Pre fault
postfault
/A
.. //

 

 

 

 

 

 

 

In

10
. $2 _
Pm - C2 S1n 6 - 0 (2.5)
and the unstable equilibrium point 6u is defined as:
. u _
The system is assumed to remain at the prefault stable

31 until the fault occurs at t=0. .At

equilibrium point 6
t=0+ and until the fault is cleared at t-tc, the mechanical
power Pm is larger than the electrical power Cl Sin6 (t)
causing the machine to accelerate. The energy

31) (2.7)

31
Vpe(t) = Pm (5(t)-6 ) + Cl (Cos6(t)-Cos6
increases until the fault is cleared at t=tc. The accelera-
tion energy, which is the kinetic energy increase on the
inertia of machine, is:

c 31 c 31
A1(tc) = Vpe(tc) = Pm(6 -6 ) + 01 (C056 -Cos6 ) (2.8)

where §a=51(tc). At to, the fault is cleared and kinetic
energy A1(tc) must be absorbed by the post fault network for
the accelerated generator to reverse direction and return to
the stable equilibrimm point 532. The acceleration energy
absorbed by the post fault network at time t is:

A3(t,tc) = v em, to) = pm(6(t)-<S°) + C2(Cos<5(t)-Cos<3c)

p (2.9)

This value of A3(t,tc) is negative and its absolute
value reaches a maximum at t=tB(tc) when the system is
stable where:

A1(tc) + A3 (tB(tC),tc) = 0 (2.10)

At t=tB(tc), the angle 6(tB(tc)) reaches its maximum excur-

sion as shown in Figure l.a and the accelerating kinetic

11

Fig. (l.b) Equal area criterion net acceleration energy
E(t) versus time.

 
    

  
  

\

’I’unstab] e

-------

 

F
i

12

energy is totally absorbed by the network and the area
A1(tc) and A3(tB(tc)) in Figure l.a are equal reflecting the
condition given by (2.10). If the fault is cleared at the
critical clearing time tcc condition (2.10) still holds but
5(tB(tc))=5u. If tc>tcc, the area A3(t,tc) is less than
A1(tc) for all t>tc and thus the velocity of the machine
angle 6(t) never reaches zero and 6(t) passes the unstable
equilibrium point 5“. For 5(t)>6u, A3(t,tc) increases and
A1(tc) + A3(t,tc) increases.
The net acceleration energy, defined as:

{A1(t) °£F<tc
E(t,tc) a A1(tc) + A3(t,tc) t>tc (2.11)

increases during the fault Ogtitc and decreases immediately
after tc when the fault is cleared. If the system is
stable, the net acceleration energy satisfies

E(tB(tc).tc) = 0 (2.12)

when the angle 6(t) reaches its peak excursion 6(tB(tc))
where both tB(tc) and 6(tB(tc)) depend on tc. If tc>tcc,
E(t,tc) never reaches zero for t>tc, as shown in Figure l.b.
Thus if E(t,tc) determines the minimum as a function of t
and defines the value as:

E*(tc) 8 min E(t,tc) (2.13)

then E*(tc) should equal zero for all tc<tCC and should be

an 1ncreas1ng function of tC for tc>tcc.

13
2.2 Potential Energy Boundary Surface, (PEBS) Method

The potential energy boundary surface method is now
introduced for single machine infinite bus system model.
The PEBS method assumes that there is a maximum potential
energy absorption capability of the post fault network.
This maximum potential energy absorption capability is never
approaches for tc<tcc but is utilized in an attempt to
decelerate the accelerated machines for tc>tcc. Figure 2.a
shows the change of generator's angle for different clearing
times. The system remains at 681 until t=0 and reaches a
maximum angle 6(tB(tc)) at tB=tB(tc) for a clearing time tc.
A measure of the energy absorbed by the post fault network
is

v*(tc) = Pm(5(tB)-531) + C2[Cos5(tB)-C05531)] (2.14)

as shown by the area in Figure 2.b. The value of V*(tc)
increases for tc<tcc and remains constant for tc>tCC for a
single machine infinite bus model. The energy V*(tc)
actually decreases for tc>tcc in the multimachine model
making the function V*(tc) reach a sharp maximum at tcc’ In
the multimachine system the function V*(tc) is evaluated
based on the potential energy component of the single
machine energy function given in equation (2.24). This
potential energy boundary surface method for multimachine

power systems will be shown to be related to the integral

criterion method developed in Chapter 5 to determine when

14

Fig. (2.a) Maximum values of generator angles for
- different clearing times.

a...
I
I
I
l

 

 

I-finm

15

Fig. (2.b) The potential energy stored in the system for
different clearing times.

 

post fault
I’-

 

 

 

 

 

16
loss of transient stability based on the region of stability

defined in that chapter.

2.3 Review of Literature on Energy and Lyapunov Functions

The two methods discussed in Sections 2.1 and 2.2 of
this chapter present a physical understanding of how the
single machine infinite bus system retains or loses
stability. The practical power system model includes
hundreds of generators. In order to apply the stability
assessment methods developed in sections (2.1) and (2.2) to
a multimachine power system, Lyapunov's direct method has
been used. The initial research focused on development of
energy functions that could be used in this method.

The first energy functions for classical transient
stability model were derived by Magnuson and Aylett [1, 2]
via the energy integral method. Later on, the systematic
construction of Lure type Lyapunov energy functions for
power system was proposed by Willems [3]. This Lyapunov
function did not contain transfer conductances and it was
realized that these terms were not small and could not be
neglected since they accounted for the load impedances.

There has been an attempt to prove that the energy
function which includes the transfer conductances is a
Lyapunov energy function [4], and another attempt to
approximate this energy function with a Lyapunov type energy
function [5], but none of these efforts have been

successful. Moreover, this energy function has several

l7
drawbacks in addition to the fact that it has not been
proven to be a Lyapunov function. The first of these
drawbacks is that this energy function is not topological,
which means that the transmission network has been
aggregated back to internal generator buses. Thus all load
buses have been eliminated and all the elements in the
resulting transmission network do not represent the actual
components. The aggregation masks the effects of operating
conditions (load level and distribution, unit commitment,
generation dispatch, network configuration, and load flow)
on stability because aggregation does not leave the network
and its operating conditions intact and reflects the actual
network operating conditions throughout the aggregated
network. It is thus difficult to access the specific
operating conditions that directly contribute to loss of
stability for a specific fault.

The second drawback is that the load at each bus is
assumed to be constant impedance since this assumption must
be made to eliminate the load buses. It would be desirable
to allow real power and reactive power load models to be
selected separately and independently at each load bus. The
real and reactive load should be either constant impedance,
constant current, constant power or a combination of these
or some polynomial function of voltage and frequency.

The third drawback is that in most cases the generator

model does not include flux decay effects.

18
Kakimoto [6], developed a Lure type function for a

nontopological transient stability' model with flux decay
using the generalized Popov criterion developed by Moore and
Anderson and a construction method developed by Willems.
The energy function presented in [6] is global and has all
the drawbacks previously mentioned.

Rastgoufard [17], developed. a :single machine energy
function for the nontopological model of multimachine power
systems. This energy function describes the kinetic and
potential energies for a single critical generator, but
since it is derived for an aggregated model it masks the
effects of operating conditions in the transmission network.
The region of stability determined based on the single
machine energy function is more accurate than can be
obtained based (n: the global energy function since the
single machine energy function contains all the information
required to predict retention or loss of stability if the
critical generator for which the energy function is written
can be properly identified [17].

Bergen and Hill [7], developed a topological Lyapunov
energy function by assuming that the loads were small
generators with small inertias which disappear after the
energy function is constructed. The resulting energy
function when load bus generator inertias are set to zero is
a constant real power load model (with a term that depends
on frequency). vo1tage at both generator and load buses is

assumed to be constant.

l9
Athay [8], constructed a topological energy function by

attempting to account for kinetic energy in the generators
and the potential and magnetic energies of the transmission
grid and loads. Real load power was permitted to be a
function of voltage and thus the load bus voltage was
allowed to change. The reactive load energy was ignored and
thus the resulting energy function was not complete.

Mussavi [9], developed an energy function by hypo-
thesis. fha did not attempt to derive kinetic energy,
position energy and magnetic energy of the transmission
grid, from the first principles and thus the kinetic energy
term given is with respect to a synchronous reference and
position and magnetic energies are based on an inertial
center reference. A very important insight of this paper is
the inclusion of the reactive energy of the load which has
never before been included. This term is necessary as
indicated by the analysis in [9] that shows the resultant
energy function is conservative. The drawback in this
energy function is that the kinetic energy is based on a
different reference than the other energy terms.

Sastry [16], gives a hypothesized topological energy
function for a model where flux linkage and saliency effects
of generators have been considered. The load modeling in
[11] is very practical. There has again been no effort to
derive a Lyapunov topological energy function for the

recommended modeling.

l

20

This thesis extends the energy function of Mussavi by
deriving the energy function from first principles rather
than hypothesizing that form. The kinetic energy term is
thus based on the same inertial center reference as the
potential energy terms. The load model at every bus allows
real power to depend on voltage which was not true in [7].
Finally the real and reactive loads are derived for the case
where both real and reactive loads at each bus can be unique
percentage of constant power, constant current, and constant
impedance load models. The energy function for the special
case where reactive load is constant current and real power
load at every bus is constant, is derived as a Lure type
Lyapunov function using the Popov criterion of Moore-

Anderson [12] and the construction method of Willems [3].

2.4 Lyapunov's Direct Method and Present Energy Functions

Since the function of this chapter is to facilitate the
understanding of the development in the following chapters,
the Lyapunov's direct method is reviewed:

The equilibrium point 0 of a system with model

3 =- yy (2.15)
is asymptotically stable if there exists a continuously
differentiable function V(x) which is locally positive
definite and for some real s>0 satisfies:

v(3_)<o (2.16)
for all

||x||<s (2.17)

21
This method could be applied to power systems, but the

problem is to identify the equations of system as (2.15) and
then be able to construct an energy function with properties
(i) and (ii). Based (“1 specific models used, several
Lyapunov energy functions have been derived for power
systems as indicated in the previous section. There are two
kinds of power system stability models for which Lyapunov
functions have been derived:
1) an aggregated network model

2) topological model.

2.4.1. Lyapunov Energy Function for the Aggregated
Network Model

In the classical power system model, all of the network
and load buses are aggregated and only the internal gene-
rator buses are retained as shown in Figure 3. The buses
and branches in the dotted box in this figure are aggregated
based on the assumption of constant impedance load at these
buses. A constant voltage behind transient reactance
generator model is also assumed in the classical model.
Willems [3] derived the Lyapunov energy function for the
classical model. Kakimoto's Lyapunov energy function [6] is
an extension of this work and assumes that the generators
are modeled by a single axis generator model. The network
is aggregated back to internal generator buses as in the
classical model. The power system transient stability model

used in [6] is

22

 

  

Transient
reactance

 

 

 

 

 

 

 

 
 

....)

Load./4'

(— ————— — — *1
generator I generator terminal bus I
Internal bus ' Load bus 1
_.__ I
I l
C»! I TRANSMISSION I
I NETWORK |
I . ‘ I

—— 0 .
'1‘ I . ' I
I t 3 I
l vflgi ‘ \iui. ‘
‘l
I
I

 

 

Fig. (3) Power system model for: the classical model
all the elements inside the doted box are

aggregated.

23

N
and
dEi o m o o
Tdoi'_3t = -(Ei-Ei )-(xdi-x di) jil Bij(E] Cos 6ij
-830 C05 613-) for i=1'oooo’N . (2019)

where for each generator i:

Pmi: mechanical power input

Yij igij: Post fault admittance between the ith and jth
internal generator buses of the aggregated
network

eij: compl1ment of ¢ij

Eilgi: internal voltage

(3”: 6.-6.
1] 1 J

E L§.: voltage magnitude and angle at the internal

q1 1 generator bus, 61 is referenced to
For the modeling given by (2.19), Kakimoto derived the
following Lyapunov energy function using the Moore-Anderson

theorem and the construction method of Willems [3].

2
V = ( ) Z 2 M-M. (w.-m.) +
2 “‘1 i=1 j=l 1 3 1 3
N N o o o o
C O O O Q .- 0 U - 6. 0-6 O I I I
i=1 §=IBIJ[E1E3(C036 1] C0351J ( 13 lJ)E 1E J

. o o o o
S1n 6 ij (Ei E i)(Ej E j) Cos 6 ij] +

O ;

24
The superscript "0" denotes the stable equilibrium point of
the post fault system.

In order to obtain his Lyapunov function Kakimoto
assumed that each internal voltage lags behind the q axis of
each generator by a constant angle ¢, and the transfer con-
ductances in the reduced admittance matrix are negligible.
Despite the significant contribution of this research, it
has the following deficiences:

a) assumes that the real power load can be
modeled as constant impedance in order to
aggregate the network and then assumes the
load is constant real power which is incor-
porated in P i' Aggregating the load back to
external buses and assuming it to be a
constant power load model, significantly
modifies the model and the critical clearing
times obtained:

b) aggregates the network and thus loses the
ability to determine or understand the funda-
mental cause of the loss of stability and the
operating condition changes that would improve
the stability margin for that particular
fault;

c) ignores the reactive load and the magnetic
energy stored in the loads that can influence
the stability of the system as voltage changes
occur:

d) neglects the transfer conductances of the
equivalent branches in the network.

2.4.2 Individual Machine Energy Function

Another type of energy function derived based on the
classical model of power system is the single machine energy
function. The single machine energy function has been
derived in Michigan State University by Rastgoufard [17] for

a synchronous reference frame model and by Fouad [18] using

25
a center of angle reference transient stability model. The
model used for this single energy function is identical to
that used by Willems [3] and neglects the flux decay in
(2.19). The voltages at generator internal buses are thus

constant and the model has the form

and
61 a mi (2.22)
where:
N
iii
a _ 2
Pi Pmi Bi G11
P = mechanical input power

mi

E. = constant voltage behind transient reactance

0')
ll

rotor angle

w. a rotor speed

M. = moment of inertia

The energy function derived in [17] has the form

V=LI§ M(w-w)2-1—-D§]I [PM-PM]
1 MT 3:1 3 i ' MT j=l i j j i
6i=6j+6o
sl

[Sij-Gij] - CijICos6ij-Cos6ij] + DijCos6ij d(5i+5j-5o)

5 Sl+5 $1_ 5 51
i j o
where

1 N N
5 = -— 2 Mo5- and MT = Z M

° "T i=1 1 1 i=1

An equal area method and a potential energy boundary
surface method [21] were derived based on this energy
function. These methods are similar to those described in
Section 2.1 for the single machine infinite bus model.
These methods have no detectable error in determining the
critical clearing time once the proper critical generator is
identified in order to write the energy function and then
apply these methods. The theoretical development of a
region of stability and the region of instability for this
model, based on the research performed in Chapter 5, is
underway, which would hopefully provide a theoretical
justification for the equal area and potential energy
boundary surface methods developed for this individual

machine energy function.

2.4.3 A Topological Lyapunov Energy Function
A Lyapunov energy function for a topological power

system model, where the network is not aggregated back to

27
internal generator buses, was derived by Bergen and Hill

[7]. The power system model has the form given below:

.. N
o . - o
M16i + D16i + §=lb13 S1n (6i 6j) - Pmi
jfi for i=1,..., Ng (2.25)
and for each load bus i:
. N o
jfi for i=Ng+1’....’N (2026)
= o '
where PDi PDi + D16i (2.27)

and

Mi: generator inertia constant

61: generator internal angle or angle of voltage of
bus i

Di: damping ratio

bij: Bijvivj

Pmi: mechanical input of generator i

PDi: electrical real power of bus 1

Bergen and Hill [7], obtained the following Lyapunov

function based on their topological model of power system:

N L 0
V = l 29 M w2 + 2 b I k (Sin u-Sin6o) du
2 k k k k
k=l k=1 of

(2.28)
where L is the number of the lines and Ng is the number of

generators.

28
Bergen and Hill [7], derived their Lyapunov energy

function using this model. This work is a significant
contribution especially in giving the opportunity to under-
stand the reason of lack of synchronism in the power system
in the case of instability. However, it has some deficien-
cies:

a) Real load power is not a function of voltage.
The real power for the sake of derivation has
been assumed to be a function of frequency but
it doesn‘t have any dependence on voltages.
This assumption significantly affects the
critical clearing time predicted.

b) The flux linkage and field effects of genera-
tor model have been completely ignored. This
will cause an error in kinetic and potential
energy terms.

c) Load bus voltages are kept constant. This
assumption is inconsistent with any transient
stability simulation model and would cause a
conservative critical clearing time since load
would not decrease during the fault.

d) The reactive load power has been ignored in
this model so the magnetic energy due to

generators and reactive power load is missing
in the potential energy term.

2.5 Summary
The literature on the application of Lyapunov's direct
method to power systems has been reviewed. The focus of

this literature review has been on the development of

Lyapunov energy functions in order to make apparent

l) the contribution made in the derivation of
energy functions for a topological network
model that permits a general real and reactive
load model in Chapter 3;

2)

3)

29

the contribution of deriving the Lyapunov
energy function for a topological network
model , constant real power load model ,
constant current reactive load model , and a
single axis generator model in Chapter 4;

the contribution made in Chapter 5 in deriving
the region of stability and region of insta-
bility for a power system based on the Popov
stability criterion assumed to be satisfied in
the derivation of the Lyapunov energy function
for the topological model in Chapter 4.

CHAPTER 3

 

TOPOLOGICAL ENERGY FUNCTION

A topological energy function is derived in this

section that
l) retains the actual network and thus does not
aggregate the network back to internal gene-
rator buses;
2) allows a general real power load model that
can include constant impedance, constant
current and constant power components;

3) allows a general reactive load model that can

include constant current, constant impedance
and constant power components.

The energy function is constructed using the integral
method that associates with each term in the differential
equations that describe the system. A term is added to this
energy function to account for the reactive load at each
bus. The energy function derived based on the integral
method for the special case of constant current reactive
load and constant power real load can be shown to be a
Lyapunov function based on the results of Chapter 4.

Although the form of the energy function has been
hypothesized in both [9] and [16], it has never been
analytically derived from the integral method, Lyapunov
construction methods, or any other procedure. The energy
function hypothesized in [9] thus utilizes a kinetic energy

term based on the synchronous reference frame but utilizes a

30

31

potential energy function based on a center of angle
reference frame. This inconsistency could not occur if the
energy function were analytically derived. The energy
function of Mussavi [9] was the first to successfully model
the reactive load energy component and this model is also
used in our derivation. A recent unpublished Ph.D. thesis
[11] has hypothesized an energy function that is very
similar to the one analytically derived in this section.
The energy function derived in this thesis and in Sastry
[16] permits use of a general reactive load model and
specifically itemizes the magnetic energy associated with
the transient reactance of the synchronous generator model.
A constant impedance load model was used in [9] and no
accounting for the magnetic energy associated with the
synchronous generator transient was made. The real power
load is assumed to be a general load model in this thesis
but was restricted to be constant real power in [16].

The energy function is derived in Section 3.2 and then
is proved to be conservative in Section 3.3 for the case of
constant real power load. The results of a simulation of a
39 bus system for two different fault cases each using a
constant current and constant impedance real and reactive
load model are presented in Section 3.4. The kinetic
potential, and total energy is plotted for the critically
stable and unstable simulation runs. The components of the
potential energy are also plotted for both stable and

unstable simulation runs for each fault case with both the

32

constant impedance and constant current load models.

3.1 Derivation of the Energy Function

System modeling:

a) Generator model:

It has been assumed that the voltage behind transient
reactance in each generator is constant. The swing equa-

tions are of the form

M161 = Pmi-Pei (3.1)
where
S 0
31 = Bi (3.2) for 1 = 1,2,..., Ng
EiVi
P . = ——7— Sin (ai-ei) (3.1-a)

e1 xd i
b) Load model:

The loads include both real and reactive power compo-
nents which are nonlinear functions of the applied voltage
and frequency. However, the constant real power case is
considered also. It has been assumed that there are
ficticious generators connected to the load buses which are

governed by

0 ’ Qli-Qi
6 = - .- .
631 1 P11 P1 where i=1,....,N (3.3)

c) Transmission system equations:
For the assumed transmission system the power flowing

through the i-th bus can be written for load bus 1.

33

 

 

N
Pi - Z_ Bijvivj Sln (61-6j)
3-1
(3.4)
N
Qi = - §=1 BijViVj Cos (5.1-5.3) (3.5)
for generator bus 1:
P. = Sin (e.-5.) + V.V.B.. Sin (e.-e )
1 "7Fi;;' 1 1 ggl 1 J 13 1 3 (3.6)
V2. -E.V. Cos (6-‘6-) N
Q. = 1 1 1 1 1 'XB .v. v. Cos (ei -ej )
1 , ij i j (3.7)

x di j= -1
Derivation of the energy function:
The swing equation with respect to a center of inertia
is formulated as
“116 a PcoA
Equations (3.1) and (3.2) could be written in terms of

(3.8)

the center of angle reference 11‘11'10‘35

. M.
2 = _ _ _1_ .=
Mini Pmi PCi PcoA for 1 1,..., Ng (3.9)
"T
1 Egi. ' 3 10
T

Multiplying both sides of (3.9) and (3.10) by (Di and
then summing both equations and then taking the integral of

this sum from ts to t the new energy function will be

NG N Ng
_ l - _ 2 1 2 5
VB — 2 2- 111 (mi mo) + 5 §_ emi(wi-mo) “i Pmi (51-51)
1-1 1- =1
g (6(t) g It
+ Pl. d6. + 5.P.dt
1:1 6(t ) 1 1 i=1 t 1 1
S S
N
29 t w P dt + Q
+ . .
1:1 1; 1 1 L (3.11)
S

The term QL is added to include the reactive energy in
the load that must be included to make this energy function

conservative one must utilize the fact that

{1 Em. a3. + 3:19 34.53.30

i=1 1' 1 i=1 1 1
in deriving (3.11).

The first term in (3.11) represents the kinetic energy
of generators and the second term represents the kinetic
energy at the load buses. By setting e‘—- 0 the second term
will vanish. Setting eu—e 0 eliminates the ficticious
generators at load buses that have been included here solely
because they must be included to derive the Lyapunov energy
functions in Chapter 4. The third term in (3.11) represents
the potential energy component associated with linear
displacement of rotor angles and the fourth term is
representing the potential energy associated. with linear
displacement of load bus angles. These two terms are called

the position components of potential energy. The fifth term

35

stands for the magnetic energy stored in the network. This
term will be discussed in the next section. The last term
represents the change in magnetic energy stored in the load.

Thus whene:—o 0 the energy function becomes:

 

N
l g 2 g s
= - - 6 -6
vs 2 {_ Mi(wi mo). 2_ Pmi( 1 i) +
1-1 1—1
E ei(t) g t
Pl.de. + 6.P.dt + Q
131 I 1 1 i=1 I 1 1 L (3.12)
61(ts) ts
N t
The term 2 ‘Lt BiPidt, represents the magnetic energy
i=1 8
stored in the network. The magnetic energy“' of each line
is:
1 2 ...—1 2 a 1 ‘W'ij’k‘vfl"
1 (I I2 I I I I ~

considering (3.13). (3.4) and (3.6)

N t l N N s s
X I mi? dt = -2 E 1 vivjsij Cos 913 + V1 vj
t
1.31 3 i=1 j=1
s “g s s
' ' .. 2 2 _ .. .-
Bij Cos e 13 + gal N1 + Bi 2 EiVi Cos (Si 51)
52 51 _ S s s s L
(81 + V1 ZEivi C05 (61 '51 )1 2x6. (3.14)
l

 

*The energy stored in an inductor is viewed physically as
kinetic energy based on the electron motion, but when
interpreted in this model is viewed as potential energy.

36

where B1 = - B.. and Bii does not include generator tran-

11 13

erraz

sient reactances for i=l,2,..., Ng or load impedance for
i=l,2...,N. The first summation in (3.14) represents the
magnetic energy stored in the branches that are connected to
nongenerator buses, and the second summation gives the
magnetic energy stored in the generator transient
reactances.

The last term in (3.12) represents the magnetic energy
of the load buses. Letting the reactive power at load bus 1
to be a nonlinear function of voltage Vi' 11(V11' 0 could

L

be written as:

N V- f.(x.) dx.
0L =2 fs1 _1_1_ 1 (3.15)
i=1 Vi xi

The right side of (3.15) is justified as a measure for load

energy since represents current andeidVi is a

1
measure of energy.

3.2 Proof of Conservativeness of the Energy Function
By substituting the results of (3.14) and (3.15) in

(3.12) the new form of energy function will be:

N N
VE = 1 lg M ((1) -(1) )2 - g P (O '58)
2 1:1 1 i 0 i=1 mi 1 i
N V.
i=1 t 1 dt t 1-1 s x. x
S . Vi .1

-% Y E Bi.(ViV. Cos ei.-v: v? Cos 9:.)
i=1 j_1 3 J J J J
N
g 2 2 2
s 2 s s s
+ i=1 ((131 + vi -2EiVi Cos (6i 91)) - (Ei + vi
-2ESV3 Cos (67-65)) 1
i i 1 i 2x5i

(3.16)

The first term in (3.16) is representing the kinetic
energy and the rest of the terms represent the potential

energy. Thus the potential energy in (2-9) PE, is:

 

N9 5 N t dei
1=l 1=l ts
V - N N
+ 1 I l f1(Xi)dxi - l E Z B..(V.V. Cos e..~
1:1 3 x 2 1=1 '=1 13 1 3 13
v. i 3
1
N9 2
s s s s 2 s
Vivj Cos eij) + i=1 [((Ei +Vi -2EiVi Cos (6i-ei))
s2 s2 s s s s l
- (Ei +vi - 2EiVi Cos (Si—61))1 3;,—

di

(3.17)

38

The energy function presented in (3.16) is a conser-

d VE
vative energy i.e., at = 0. This has been shown in the

following discussion:
Differentiating PE in (3.17) with respect to bus

voltage vi:

8
a)“ —v.—' 3?": ' 7— ‘ Bijvj C05 13-
1 1 d1 di j=1 (3.18)

for 1.11000! NC

a - B..V. Cos 6..
av v E 11 u
i 1 3:1 for i=N 'ooo’N (3.19)
+1
Thus
dPE ._ .
av;.= (Q1i + Qih/Vi for 1—1,2,....,N (3.20)

since by definition Qli= +fi(vi) is the reactive load of bus
i, and Qisatisfies (3.5, 3.7). It is clear based on (3.3)
that

i i=1,2,...,N (3.21)

Now differentiating (3.17) with aspect to 6i:

d N Eivi Sin (5i-ei)
PE _ . _ —7—
36—; - + §=1 Bijvivj S1n eij xdi

for i=l,....,Nb (3.22)

39

N
g3§ Z Bi.ViV. Sin ei.
91 3:1 3 J J

for i=Ng+l'°"°'N (3.23)

Considering (3.4) and (3.6), and applying their results

to (3.22) and (3.23):
do 1 i=1,2,3,....,N (3.24)

It can also be shown that

d8

N
dPE = i
“HE i=1 P11 HE‘

i=1,2,3,....,N (3.25)

Combining (3.24) and (3.25), one can show that

E dPE , dei + dPE = § (P +91 ) dei = 0
1:1 dei at dt 1:1 at
i=1,2,....,N (3.26)
since
Pi + pei = o i=1,2,3,....,N

from (3.3).

Now differentiating with respect to 6i:

dgE = -Pm. + Eivi S1n (6i-ei)
a . 1 xi.
1 d1

for i=1'oooo'Ng (3.27)

Recalling that the kinetic energy is:

'g 2

N
KE ‘”2 1:1 Mimi-(1)i )

40

and then differentiating it with respect to mi produces

BKE o
—&’—i-=Mi((lli-Qli)
for i-l,2,....,Ng (3,23)
Combining equations (3.27) and (3.28), (N) i can show
that
N9 3K ‘1‘1’1 an “1 $19 .
~ ' + ' = (M a - Pm.+P )w = 0
§=1 5&1 at 36i at 1:1 1 1 e1
(3.29)
since
ani a Pmi-Pei
(3.30)
Considering (3.21). (3.26) and (3.30)
(IVE- (3.31)
’3? 0

and thus the energy function has been shown to be conserva-

tive.

3.3 Energy Function with Different Load Models
3.3.1 A General Real Power Load Model_
In the case that real load is represented by nonlinear

function of applied voltage, the second term in (3.16) could

be written as:

41

N N
Z (6(t) Plidei = % { (6(t)-8(ts))(Pli(Vfts)'
1=Ng+1 6(tS) 1=Ng+1

+ Pli(Vi(t)) (3.32)

using the trapezoidal integration rule. Substituting the
right hand side of (45) the following generalized function

is obtained:

N _
N V(t ) V(t)
_ 1 ~ 2 l _ s
9 ' N V.(t)
- 1 Pm.(6.-6§) + , 1 fi(xi) dxi
- 1 1 1 ______
1=l 1=l V (t ) x
i s 1
1 N N s s s
-.. B i V.Vo C0560 u - V.Vo C ea c
2 i=1 1:1 13( 1 J 13 1 3 OS 13)
N
2 2 2

9 s
+ z [(Ei-Wi -2EiVi Cos (61-61))-(Ei +Vi
i=1

X .

3 5- 3 . 1 3.33)
-2EiVi C08 (51 Gi)fl i—El (

3.3.2 The Energy Function With Different Reactive Load
Models

 

The reactive power at load buses is assumed to be a
general nonlinear function of voltage. Initially the case

of constant current reactive load is considered.

42

a) Constant Current Reactive Load Model
In the case that reactive load is constant current

reactive power of the load will be:
-0:

f11xi) = aixi where ai = V:

 

and Qi8 is defined in (3.5) and (3.7). The energy

associated with the reactive load 0L will thus have the form

V. f.(x.) N V
.._1_l.._ dx. = —1 Q.
X- 1. =1 1

(3.34)
Substituting (3.34) in (3.33) the following energy function
will be obtained:

N V(t ) V(t)

N
3 1 9 _ 2 1 2 (e.-e?) (Pl$ s + p1. )
VB 2 E31 “1(“1 ”0’ + 2 i=1 1 1 . 1 1
N
9 N V N
- X Pmi(Gi-6:) _ Z (1-—§) f VEVSBi. Cos 8:.
i=1 i=1 vi j=1 3 3 3
N9 2
i s s s s s l
+ 1:1 (1-53) (visi Cos (ai—ei) - vi ) 131
i
N
g 2
2 2 2 s
+ 1:1 [(Ei +Vi -2EiVi Cos (6i_Bi))- (Ei +Vi
2 N N
s s s l 1
-2E.V. Cos (6.-6.)) -—7 - - 2 Z B..(V.V. Cos (6..)
1 1 1 1 (Zxdi 2 j=l i=1 1] 1 j 1]
- VFVF Cos (8.3)) (3.35)

1 3 13

43
b) Constant Impedance Load Model.

In the case of constant impedance reactive load model,

the function 11(V1’ is a quaradic function of voltage:

Hwn

f - 2 h Q
i(xi) - aixi w ere ai - - 2

V

 

rmm

and the energy associated with the reactive load QL has the

form
2
N V. f.(x.) N V. N
1 1 1 l s 1 s
Q: ———dx.=- g Q.—-——+ Q.
1 i=1 IVs x1 1 7 i=1 1 vis i=1 1
i

(3.36)

Substituting (3.36) in (3.33) the following energy

function will be obtained:

N
9 N \I(t ) V(t)
_ 1 _ 2 l _ s s S
N 2
9 N V.
- Z Pmi (61-6?) + Z (1- 1 ) Z viv331. Cos 913
V.
1
N 2
9 V1 3 s s s 2 1
+ I (l- _—_2)(ViEi Cos (6i-Bi) - Vi ) -;7-
1=1 Vs d1
1
N
g 2 2 2 32
+ z [(E.+V. - 23.v. Cos (5.-e.)) - (E.+V.
. 1 1 1 1 1 1 1 1
1=l
2 N N
s s s l 1
- 2E.V. Cos (6.-6.)d' - Z Z B..(V.V. Cos 6..
1 1 1 1 2x5i 2 1:1 j=1 13 1 j 13
s s s
- Vivj Cos eij) (3.37)

44

c) Constant Reactive Load Model.

When the reactive power of loads are constant then:

s

i

and the energy associated with the reactive load 0L has the

form

N . ) N
- l 1 i _ _ s s

1
(3.38)

Consequently the form of energy function will be:
N
1
VB = 2 1
i=1

9 V V

o 2 l s s
“1‘“1'11’ + 2 (91'61) (P11 + P11)

Fonz

'=l

N
9 N y
_ _ s s s s 3
Z pmiwi 61) + Z Ln(Vi/Vi) . vivjs.. Cos sij

i=1 i=1 . i=1 1]
N
g 2 1
s s s s s s ) ° —-
1 i=1 1n(V1/V1’ (ViEi C°S (51'611'V1 xdi
N9 2
2 2 2 s s s
+ g: B51 +vi zsivi Cos (51'91”“Ei +vi cos (61-61))
.%;,_ - % § § Bi.(ViV. Cos ei.-v§v§ Cos 6:.)

(3.39)

45
d) The Generalized Load Model.
In the case that reactive load is a combination of
constant current, constant impedance and constant VAR Model

the function fi1xi) is of the following form:

2
£i1xi) a aixi + bix1 + C1 (3.40)
where

s s
Q Q-

._ i _ _ J. _ _ 5

a1 - a1 2 b1 - 81 s and C1 - 7iQi
8 V.
Vi 1

This model assumes that a percentage of the load is

constant impedance (01), constant current (Bi) and constant

power (Y1) such that:
(11 + Bi +‘yi . 1 (3.41)
The total energy function in this case has the

following form:

N
1 9 o 2 1 " s V V
v: - Z N (m -m ) + (e -e ) (913 + pl
2 1-1 1 1 1 2L1 1 1 1 1’
"9 N
- p (6 -6’) + Ln(v v”) v’vss Cos 63
1.1 “'1 1 1 *1 (.1 i/ 1 .1 1111 1)
N9 2
+ r Ln(V v3) (v32 Cos (53-93) - vs ) ° 1,—
1_1 1/ 1 1 1 1 1 1 xd1
N v N v
1 s a s g i s s 32
+8[2 (1- —) Z VVB Cone +12 (1-——)(VECosé-8 )
1 1.1 v: j_1 1 j 1] ij 1_1 v: 1 1 1 1
2
N v
s l i s s s
-V)‘—r-+u (Z (1- )2 VVB C030)
1 ”d1 1 1-1 vs! j-l 1 5 13 15
1
N
9 vi a a s s 32 l
+ z (1--—I (v 3 Con (5 -e )-v ) - -4—

c (5 -e ))-(z 2+v32-21: v52 c (53-93)
°‘ 1 1 1 1 11 °9 1 1)

.1.

9 2 2
I (E +V -ZE V
1_1[ 1 1 11
’1 N N
. l s s s
(3.42)

46

3.4 Simulation Results

The components and total energy for the energy func-
tions presented in (3.35) and (3.37) are plotted as a
function of time based on transient stability simulations of
particular fault on a 39 bus New England system shown in
Figure 4. The purpose of computing and plotting the energy
functions for both a constant current and a constant
impedance reactive load model for two different faults is

1) to indicate that the total energy remains
constant after the fault is cleared when the
real power load model is constant power
confirming that the energy function is conser-
vative when the real power load is constant;

2) to show the relative magnitude of the kinetic
and potential energy components and show that
since the total energy is conservative that
an increase in either kinetic or potential
energy results in a decrease in the other
energy component;

3) to display the relative magnitude of the
position, network magnetic, real power load,
and reactive power load components of the
potential energy as a function of time. The
differences in the potential energy components
as well as the swing curves for constant
current compared with constant impedance
reactive load models will also be shown;

4) to show how a loss of stability causes a
sudden sharp drop and oscillations in the
various potential energy components as well as
a sharp increase in the kinetic energy in the
system. The total system energy no longer
remains constant once the loss of synchronism
occurs and may be due to numerical problems
caused by the very large and rapid changes in
the system after the loss of synchronism.

The two fault cases considered are a fault on line
26—27 where line 26-27 is removed and a fault on line 21-22

where line 21-22 is removed. It will be observed that the

47

Fig. (4). New England Study System.

 

 

 

 

23

 

 

14

%
j
1

 

 

 

 

 

 

 

1 n [:—___l I
n
N
1.-
d)

 

 

 

 

 

21
15

ll
L

 

 

 

 

 

 

 

 

 

 

 

n
_ 12
i u I

I!

 

 

3%

 

32

 

 

 

 

 

 

 

 

 

33
_34

 

 

 

 

 

i.
g

39

 

 

 

48

critical clearing time for the fault on line 21-22 is much
shorter than for the fault on line 26-27. The results for
the simulation of the fault on line 26-27 with that line
removed is discussed first for both the constant current and
constant impedance load model. The results of the second
fault~case for both the constant current and constant
impedance load models are then discussed. It should be
noted that the MVA base used is 10 MVA base due to the fact
that the Detroit Edison transient stability program used
this MVA base.

1) Fault on Line 26-27 and Constant Current Reactive

Load Model.

In this case the fault occurs on line 26-27 and the
critical clearing time by simulation has been determined as
.36 seconds.

Figure (5.a) shows the variations of angles of gene-
rators close to the fault on the stable case. Figure (6.a)
shows the variations of the same generator angles for the
usntable case. In the stable case the angle of generator
No. 9 increases until it reaches 180 degrees and then starts
decreasing and continues oscillating until it remains
constant. For the unstable case the generator No. 9's angle
continues to increase and is the first generator that loses
syncronism. Figure (5.b) shows the total, potential and
kinetic energies for the stable case; the same kind of
energies are shown in Figure (6.b) for the unstable case.

In both cases whenever the potential energy is changing the

49

kinetic energy experiences a negative change of the same
magnitude such that the total energy remains constant. The
small variations in total energy are due to the computa-
tional errors and also the assumption that real power
delivered to the load buses is constant. Whenever the
potential energy is at its maximum the kinetic energy has
its minimum value. This observation supports the P888.
method discussed in Chapter 2. In the unstable case when
the system loses synchronism the potential energy
experiences a large drastic sudden drop and the kinetic
energy starts increasing. Furthermore, there is a sudden
drop in total energy which may be due to numerical problems
in attempting to accurately simulate the large rapid changes
after the point in time when the loss of synchronism occurs.

Figure (5.c) shows three components of potential energy
for the stable case. These three components; energy due to
real power load, energy stored in the transmission grid, and
energy due to the generator's angle displacement, are shown
for the unstable case in Figure (6.c). In the stable case
the energy due to the generator's angle displacement
decreases until minimum is reached and then starts
increasing again. In the unstable case, since the
generator's angle continues increasing, the position energy
component is a monotone decreasing function of time. The
other two forms of energy in the stable case increase until
a maximum is reached and then decrease smoothly but in the

unstable case these energies start oscillating. The sum of

50

the real load power energy and the magnetic energy stored in
the transmission grid is larger as a function of time in the
unstable case than in the stable because the system
trajectory's peak angle excursion is larger in the unstable
case.

Figure (5.d) shows the energy stored in the reactive
power loads for the stable case. Figure (5.d) shows the
same energy for the unstable case. The amount of this
energy is much smaller than the other energies components.
Since this energy depends on variation of voltage of the
load buses in the stable case, the energy approaches zero as
time increases and the system trajectory approaches the 88?.
In the unstable case the oscillations of voltages in the
load buses cause this energy to experience large oscil-
lations.

ii) Fault on Line 26-27 and Constant Impedance Reactive

Load Model

For this case the fault occurs on line 26-27 and the
critical clearing time by simulation is determined as .31
seconds. The .05 seconds reduction of critical clearing
time due to the change of reactive load modeling is expected
because there is a greater reduction in the magnetic energy
stored by the reactive load clearing the fault when voltage
drops. This reduction in magnetic energy results in
increased kinetic energy because the energy function is

conservative.

51

Figure (7.a) shows the changes of angles of generators
7, 8 and 9 which are closer to the fault for the stable
case. The variations of the angles of the same generators
are shown in Figure (8.a) for the unstable case. In the
stable case the angles oscillations die out, but in the
unstable case the angle of generator 7 as well as generators
8 and 9 continue increasing. It is important to note that
for this kind of load modeling, not only the critical
clearing time changes but also the first generator that
loses syncronism is different from the constant current load
model.

Figure (7.b) shows the changes of total, kinetic and
potential energies for the stable case; the same energies
for the unstable case are shown in Figure (8.b). In both
cases whenever the potential energy has its maximum value
the kinetic energy has its minimum value. This observation
again supports the P388 method. In the stable case the
total energy is almost constant. The small variations in
the total energy over time shown in Figure (7.b) are due to
computational error and the assumption that real load at the
load buses is constant. For the unstable case as shown in
Figure (8.b) the total energy has dropped after the
generators have lost synchronism. The potential energy has
a sudden drop and the kinetic energy has a sudden increase
after the loss of synchronism. This sudden drop occurs

because the simulation program is unable to give accurate

52
solution with the large and sudden changes of voltage
magnitudes and the changes in angles.

Figure (7.c) shows the energy due to real loads and
energy stored in the transmission grid and energy due to
generators angle displacement for the stable case. These
three energies are components of potential energy. Figure
(8.c) shows the same potential energy components for the
unstable case. The energy due to generators angle displace-
ments for the stable case reaches a minimum and then starts
increasing but in the unstable case this energy is always
decreasing. This continual decrease is because in the
unstable case the angles continue increasing. The real load
energy and energy stored in the transmission lines for the
stable case are smaller and smoother than in the unstable
case. In the unstable case these energies have oscillations
which are due to large and fast oscillations of voltages and
angles.

Figures (7.d) and (8.d) show the energy stored in
reactive loads for stable and unstable cases. The reactive
power of the load for each bus is a function of the square
of the voltage of the corresponding bus. By recalling the
equation (3.36) it could be seen that the larger is Vi2 ,
the more negative will be the energy due to reactive loads
and this explains why this energy in this case is negative
compared to constant current load model. In the unstable
case, there is an oscillation in this reactive load energy

as can be seen in Figure (8.d).

53
From both cases that fault was on line 26-27, it could
be concluded that the total energy stored in the system for

the stable case is less than the unstable case.

iii) and iv) For the case that fault has occurred on
line 21-22 the same tests of the last two sections were
performed, the results have been presented in Figures 9, 10,
11 and 12. ‘These experiments were performed to show that
the cases discussed in the preceeding two sections are not
special cases. The results of these tests support the

conclusions from the tests for the faults on line 26-27.

54

Fig. (5.a) Variations of angles of generators 7, 8 and 9 for
stable case where the fault on line 26-27 is

cleared in 0.35 seconds by removing the faulted
line.

 

RNGLE DEGREE

i ' w ‘3
‘ TIME IN SECONDS

 

55

Fig. (5.b) Variations of potential, kinetic and total energy
for the stable case where the fault IS on line
26-27.

total‘energy

‘\\\tpotential

energy

 

 
 
 

   
 

   

¢._kinetic
energy

 
  

EIEIOV P.U.

1!

 

 

1

”J

VINE IIiOECONOI

Fig.

=2 '1?

I.

ENERGY

56

(5.c) Variations of real load energy and transmission

115.

RS.

-IISJ

-3”;

 

-01§J

line energy and energy due to angle displacement
of generators rotors for the stable case when the
fault is on line 26-27.

real load energy
energy stored
in transmission lines

i ;
TIN! ll SECOND.

.--. ..
G“
“4
q
no

energy to generator angle displacement

57

Fig. (5.d) Variation of reactive load energy for the stable
case where fault is on line 21-22.

- ENERGY P.U.

.Js' f i V ' 12
VINE [I SECONDS

 

58

Fig. (6.a) Variation of angles of generators 7, 8 and 9 for
unstable case where the fault on line 21-22 is
cleared in 0.36 seconds and line 26-27 is
removed.

 

RNGLE DEGREE

 

 

. _. i , 4?
‘.‘ TIME IN secouos V\\‘\~l_§x<:::;)/ '

59

Fig. (6.b) Changes of potential, kinetic and potential
energy for the unstable case when the fault is on
line 26-27.

331 total energy

13‘

  
 
    

 

.ur )
"”1 kinetic
energy

ENERGY '0“.

potential energy

-ns¢ o—”

 

-us‘

60

Fn1.(6.c) Variation of real load and transmission line

energies and energy due to angle displacement of
generator rotors for the unstable case when the
fault is on line 26-27.

   
 
    

 
 

 

 
   
 

5‘
real load energy
159.!
1 energy\ stored in
. . transmission lines
1 f i i
i ' z W
-15.) )
a: 4
e. -§l01
é .
3 energy due to generator angle
- 19L displacement
-mo‘

 

It"! [I SECOND!

61

Fig. (6.d) Changes of reactive load energy for the unstable
case where fault is on line 26-27.

C
L

MOV P .0.
O

I
l
.3" ‘ -
"H! [I C

WW,

 

Fig.

62

(7.a) Variation of generators angles for stable case
that fault occurs on line 26-27. The fault is
cleared at 0.30 seconds and line 26-27 is
removed. (Constant impedance load model)

MOLE DEGREE

 

    

 

TIME IN SECONDS

63

Fig. (7.b) Changes of kinetic, potential and total energies
for the stable case when the fault is on line
26-27. (Constant impedance load model)

 

 

 

I",
‘ r“‘_‘ total energy
l
.I-( .2.
z I kinetic energy
a I
I
2 u. I
l
y ' potential energy
I I ' . . -

 

 

TIRE IN SECONDS

64

Fig. (7.c) Changes of energy due to generator angle

displacement and loads real power energy and
energy stored in transmission grid for the stable
case where fault on line 26-27 is cleared in 0.3
seconds. (Constant impedance load model)

real load energy

   
 
    
 

 

energy stored in tra ission lines

 

 

      
   

generator angle displacement

energy \\\‘

TIME "I SECOND.

65

Fig. (7.d) Changes of loads reactive power energy for the
stable case when fault is on line 26-27.
(Constant impedance load model) '

“' i i ' fl
FINE IN SECOND.

ENERGY 9.0.

l
-

A

 

 

I
Q
L;

66

Fig. (8.a) Variation of generators angles for unstable case
when the fault occurs on line 26-27. The fault
is cleared at 0.31 seconds and line 26-27 is
removed. (Constant impedance load model)

RNGLE DEGREE

 

 

 

TIME IN SECONDS

67

Fig. (8.b) Variations of kinetic, potential and total
energies for the unstable case when the fault is
on line 26-27. (Constant impedance load model)

total energ
‘ kinetic energy 1 /
2». Fr-/

>f<><=<)

 

: T v i 3.
.. nut III are
’- " 1'1 potential
§ ‘ energy
-m‘

 

—‘

Fig.

ENERGY P.“.

68

(8.c) Variations of real load energy. energy due to
generator angle displacement and energy stored in
transmission line for unstable case when fault is
(Constant impedance load model)

on line 26-27.

energy stored in
transmission lines

real load energ

  
  

 
 

   

  

 

 

i

TIME IN RE

 

, w
/

generator angle displacement
energy

 

 

69

Fig. (8.d) Variations of reactive load energies for the
unstable case when the fault is on line 26-27.
(Constant impedance load model)

    

r

TIME IN SECONDS

hr-
4>

    

 

70

Fig. (9.a) Variation of angles and generators 7, 8 and 6 for
the stable case when the fault on line 21-22 is
cleared at 0.15 seconds and line 21-22 is
removed. (Constant current load model)

115
d
95.

851

RNGLE DEGREE
to
k-

 

IS fi
' i I 3

YIHE IN SECENDD

71

Fig. (9.b) Changes of potential, kinetic and total energies
for the stable case where fault was on line
21-22. (Constant current load model)

 

 
 
 
      

  

 

 

9|
1
total energy /////’.—‘\
..4) —\
=3
e I
5 1 potential
I ener
I kinetic
1 energy
I
, \
I. ' v r . . .
TIME IN SECOND.

72

Fig. (9.c) Variations of real power load energy, energy

stored in transmission lines and energy due to
generator angle displacement for stable case when
the fault is on line 21-22. (Constant current
load model)

real power load
energy

  
       

 
 

 

HQ
31 energy stored
in transmission lines
I _. . , I.
T' I 25
-3‘
J
'28 4

ENERGY P .U.

'3" 1

 

-515
4

      
   

generator angle displace-
ment energy

\

TIME IN SECOND!

73

Fig. (9.d) Changes of reactive power load energy for the
stable case when the fault is on line 21-22.

 

 
   

TIME IN DEW

 

Fig.

(10.a)

 

LEE__._L-

MOLE DEGREE
1 -

 

 

74

Variation of angles of generator 6, 7 and 8 for
unstable case where fault on line 21-22 is
cleared in 0.36 seconds and the faulted line is
cleared.

 

 

I I
TIME IN SECONDS

75

Fig. (10.b) Variations of potential, kinetic and total
energy for unstable case when fault is on line
21-22. (Constant current load model)

total energy [y

 

 

 

 

1! ' ‘7 ITIME IN arconos i ‘3

5 potential
4 1 energy kinetic energy
S-nfl
I
3 .
U

-m‘

-53.

 

76

Fig. (10.c) Variations of loads real power energy: and

 

energy stored in transmission lines and energy
due to generator angle displacement for unstable
case when fault is on line 21-22. (Constant
current load model)

real load energy

“\a

 
   
 
  

 
 

\ .
energy stored in
transmission lines

 

  
  

generator angle
displacement energy
4r—-"""

  

TIME IN SECONDS

77

Fig. (10.d) Variations of loads reactive power energy for
the unstable case when fault is on line 21-22.
(Constant current load model)

       

Vi

I
TIM IN SECOND

'4

..J

ENERGY PM.

 

78

Fig. (11.a) Variation of angle of generators 6, 7 and 8 for
the stable case when fault occurs on line 21-22.
The fault is cleared at 0.13 seconds and line

21-22 is removed. (Constant impedance load
model)

#4"’
In ,/
. Afl" ‘\\\\
IN. /

 

E\‘
E.

\Oll.

MOLE “CREE

1 ' 1 ' 3
r1»: IN seconos

79

Fig. (11.b) Variations of total, kinetic and potential
energies for the stable case when the fault is
on line 21-22. (Constant impedance load model)

kinetic energy

 

 

a: ”I J
: :/\ total \e|nergy
§ * :-
I
I:
I /\
. I
,1! I
I ‘ TIN! In as s i 544]
4'4 potential energy

 

 

Fig.

80

(11.c) Changes of energy stored in transmission line
and energy due to generator angle displacement,
and energy stored in real power loads for the
stable case when the fault is on line 21-22.
(Constant impedance load model)

real load energy

__~\\‘\_

  
 
    
      

 

energy stored in transmission
lines

 

energy due to generator angle
displacement

 

 

TIME IN SECONDS

81

Fig. (ll.d) Variations of reactive load power energy for the
stable case where fault is on line 21-22.
(Constant impedance load model)

       

TIM IN SECOND

 

41

Fig.

82

(12.a) Variation of angle of generators 6, 7 and 8 for

——

HNGLE DEGREE

 

unstable case when the fault occurs on line
21-22. The fault is cleared at 0.14 seconds and
line 21-22 is removed. (Constant impedance load

model)

 

 

I

TIME IN SECONDS

83

Fig. (12.b) The variation of total, potential and kinetic
energy for the unstable case when the fault is
on line 21-22. (Constant impedance load model)

 

 

 

W z/\
32.. .' .
4 I
° I
a I
IE IS. .
,_ I
O 4 I
I: .
‘5‘ o '
"‘ 4T ' i ' i ' T
. . . TIME IN SECONDS
-164

Fig.

ENERGY P.D.

84

(12.c) Changes of real load energy and energy due to

’n I

 

generator angle displacement and energy stored
in transmission lines for the unstable case when
the fault is on line 21-22. (Constant impedance
load model)

  
    
    
  
 

real load energy

 
 

ener;;\§tored in
transmission Iines

 

8‘

energy due to generator angle
displacement

TIME IN SECONDS

85

Fig. (12.d) Changes of reactive load power energy for the
unstable case when the fault is on line 21-22.
(Constant impedance load model)

 
   

TIII IN SECOTDS

 

86

APPENDIX 3-A

 

Notation and system modeling:

number of buses

number of generators

real load at load bus i

mechanical power input at generator i
inertia constant at generator i
voltage of generator bus i

voltage of nongenerator bus i

imaginary part of admitance of line i-j

real power leaving bus i, through the transmission
system .

reactive power leaving bus i, through the transmis-
sion system

fictitious inertia for load buses where Mf=emn

1
N9 N
MT 31:1 Mi + sign mk
g+1
where e is pertubation factor
center of inertia angle 5 = l_ g9 M.6 +e:§ 5
0 MT 1 i k=N Mk k

i=1 g+1
angle of generator bus i

angle of nongenerator bus i

87

N
1 9 N
mo center of inertia frequency wo=fi_ 2_ M.m.+ 62 Mkmk
T 1—1 =N
g+1
~ ~ N
m1 mi = w1-m where 2_ M1 1 - 0
i—l
Di damping power coefficient
di fictitious damping power coefficient for load buses,
D =6d.
i i
Tdoi d-axis transient open circuit time constant.
xdi’ xéi d-axis synchronous and transient reactances

The symbol ""‘ when appearing above matrices or vectors is
standing for transposed.

The symbol ”3" when appearing above characters stands for
stable equilibrium point.

CHAPTER 4

A LYAPUNOV ENERGY FUNCTION

A topological Lyapunov function is derived in this
chapter for

l) a nonaggregated network model;

2) a constant real power load model;

3) a constant current reactive load model;

4) a single axis generator model that includes
flux decay.

The Lyapunov function is derived based on assuming that
a fictitious single axis machine exists at every load bus
that is eliminated by allowing the inertia to approach zero
and appropriate impedances to become infinite. The resul-
tant Lyapunov function is constructed based on showing that
the power system transient stability model is a Lure type
model and can satisfy the Popov stability criterion, and the
conditions imposed by the Moore-Anderson theorem [9] in some
domain of the state space. The construction method of
. Willems [11] is then applied to determine the resultant
Lyapunov function. The fictitious generators at load buses
are then eliminated and the resultant Lyapunov function is
shown to be identical to the energy function derived in the
previous chapter if the real power load is constant power,

the reactive power load is constant current, and the single

88

89

axis generator model is assumed to be a constant voltage
behind transient reactance model.

The resultant Lyapunov function and its region of
stability S provide the basis for determining necessary and
sufficient conditions at some time t* that determine a
region of instability; These region of instability and the
region of stability derived and discussed in Chapter 5 are
both based on the network conditions imposed by the Popov
stability criterion and provide an accuracy in charac-
terizing retention and loss of stability that was previously
impossible.

The energy components of the Lyapunov energy function
associated with different components of the transient
stability model are plotted for different fault cases in

Section 4.

4.1 Theoretical Background
The theoretical background required to derive the

Lyapunov function is now presented.

a) Popov Criteria

For the system:

it'ljgg-g yg) (4.1)
2* 2’; (4.2)
where gig) is a nonlinear function of g satisfying the

following conditions:

90
N
i) £(g) is continuous and maps RN into R

ii) for some real constant matrix Q:

_f_’(g) 59 > o for all _e_ ERN (4.3)

and _f_ (2) = .0.

for g =‘Q

1 N

iii) there exists a function V1€C mapping R into R
such that

V1(g_) 2. o for all _€_ERN

and V1(g) =‘Q for 0= 0

and for some constant real matrix Q

AV1(9_) = g”; (g) for all 5an (4.4)

b) It can be verified [6] that the transfer function
for the linear part of (4.1) is an NxN matrix:
[(3) -- g’Isyy'lg

and‘W(w) ==0 (4.5)

The Moore-Anderson theorem that establishes sufficient
conditions for a system to be asymptotically stable is now

presented.

Moore-Anderson Theorem 1 [9]

If there exists real matrices g and 9 such that:

Ԥ(s) = (EIQF) fl(s) (4.6)
is positive real, then the system defined by (4.1) and (4.2)
satisfying the Popov criterion asymptotically stable in the
large providing (§’+ Qs)does not cause any pole zero cancel-

lation with fl(s).

91

The conditions for fl(s) to be positive real are:

1) 2(5) has elements which are analytic for Re(s)>0
(4.6.a)

2) 2*(s) = Z(s*) for Re S>0 (4.6.b)

3) ZT(s*) + Z(s) is positive semidefinite for ReS>0
(4.6.c)

The notation (*) stands for conjugate of a complex

variable.

Theorem 2 [6]

 

Given fl and Q satisfying the condition of the Moore-
Anderson theorem, then a Lyapunov energy function of the

form:
V(x) = .ng 35 + vl(g) (4.7)

exists. 3 is obtained as a positive definite matrix
satisfying the following set of nonlinear algebraic

equations:

é’£+£§.='££’ “-3)
£§’£§’+AEQ_"£1VO (4-9)
EOEO*Q.EE+§._C_9. (4.10)

where L and W are auxiliary matrices,the derivative of this
Lyapunov funEEion:
. 1 ’ I ’
WK) 3 -Z[§ C - game] [E g-floggH-gqmg (4.11)
should be negative definite.

After recalling these theorems the next step is to test

and apply the conditions embodied in these theorems to

construct a Lyapunov function for a power system model that

92

1) retains the network and does not aggregate the
network back to internal generator model:

2) utilizes a one axis synchronous machine model
that incorporates flux decay for each
generator in the system:

3) a fictitious one axis synchronous machine
model at each load bus that is ultimately
eliminated from the model after the Lyapunov
function is constructed. This one axis
synchronous machine model is eliminated by
letting its inertia approach zero, transient
reactance approach infinity and the diffe-
rences between the synchronous and transient
reactance approach infinity. These three
assumptions eliminate the electromechanical
model of the machine, and eliminate the
armature reaction effects of stator current on
the magnitude of the induced voltage on the
stator:

4) a constant current reactive load model:

5) a constant power real power load model.

4.2 pLyapunov Energy Function Derivation

The notation used in this chapter is identical to that
used in Chapter 3 except that Bi is used to represent the
voltage magnitude at internal generator buses and network
buses. Furthermore, since there is no need to distinguish
the energy associated with generator transient reactance,
x81 from the energy stored in any network branch the
generator internal buses will be numbered l,2,....,Ng the
general terminal buses will be numbered Ng+l,...2Ng, and the

network buses will be numbered as 2Ng+1,...N. The generator

transient reactance will be included as

-— '3- =—_)_‘_
311 ‘ B1+Ng 1 31, i+Ng xd
Ng+N
Bi+N i+N = 7 X Bi+N ' ' x1
9' 9 j=Ng+1 9'3 d
sij = o ifj; ij# j + Ng, j and i, Ng+i

i'j 1'2'eeeepNg

Since the internal generator admittance of the
fictitious generators is zero and ultimately will disappear,
the transient admittance of these machines is omitted and
the internal generator bus of the fictitious generator at
load bus i is assumed to be the load bus i itself. These

assumptions allow derivation of the Lyapunov function.

4.2.1 System Modeling

a) Generators -

The flux linkage is considered in generator modeling
hence the generator i is governed by the following equation

[6]

.. N
. s s 0 s. e
= 0 e s - s e 5s. .

M16i + Didi §=13ij(3133 S1n 51] 8183 Sin 1]) (4 12)
° 3 N s s

= - e a- I - e as e 600- e 60. .
Ei 01(81 El) 81 §=1 313(8) Cos 1] EJCos 13) (4 13)

is“

Where i=1’eeeepNg

a. = [l-(x .-x’.)B..]/ . and B. = (x .-x’.)/ ,
1 d1 d1 11 Tdoi 1 d1 d1 Tdoi

(4.14)

94
b) Load Buses -

Loads are modeled by fictitious generators connected to

the load buses hence:

 

 

S 3 o S .
1 1 1 1 i=1 Bij(EiEj Sin Gij-EiEj Sin aij) (4.15)
and
s = —a (E -Es)- I B (E8 Cos 63 -E Cos 6 )
1 1 1 1 %.j=1 ij j ij j ij
j#i
for 1=Ng+l,...., N (4.16)
where
a = (l-(xdi-Xdi) B11’/Tdo1
1 E
= (l-(xdi-xdi)Bii)/Tdoi (4.17a)
B = (xdi-xéi)/Téo1 = (xdi'xdi)/T'
i 8 doi (4.17b)

It is also assumed that the real power delivered to each

load bus is constant and

N
P1i = i=1 BijEiEj Sln fiij = Constant (4.18)

and that the mechanical input power to each generator is
constant and given by:

_ s . s .=
Pmi — g E.EjB.. Sin Gij for 1 Ng+l'°°"

(4.19)

95
All the equations presented in sections a and b are true for

the post fault network.
4.2.2 Satisfaction of Moore-Anderson Theorem Conditions

Having defined the power system model we must

1) express it in Lure form (4.1)(4.2);

2) show that Popov statutory criterion is satisfied

3) establish that £(s) is positive real.

In order to show that it satisfies the conditions in
the Moore-Anderson theorem, having shown that the above
conditions can be satisfied, then Theorem 2 can be applied
to actually construct the Lyapunov function in the next
subsection.

The first task in derivation of the Lyapunov function
is to express (4.12), (4.13), (4.15) and (4.16) has Lure
form given in (4.1) and (4.2). The following definitions

are necessary:

x =(§;_ gfg_f) (4.20)
where

5r =51(1+1)'5171+1) for i=l,2,....,N-l

mi = 61 for i=l,2,....,N

AE = Ei-EE for i=1,2,....,N

and also:

5 (o) = §i(o) §§(0) (4-21)

96
where £1(g) is an m-vector defined by:

_ . s _ s s . s
fl(o) — Bij(EiEj Sin (0k+61j) EiEj Sin Gij) (4.22)
for i=l,2,....,N-l and j=i+l,....N
and k=l,2,....m where m = Elglll

and f2(0) is a N-vector defined by:

= s s _
f2i(o) (Ej Cos 6.. Ei Cos 51.)

13 3
for i=l,2,....,N (4.23)

U- M2

isij

The members of vector 0 used in function f1(0) are de-

fined as:

ck a 5..-6.. for k=l,2,....,m (4.24)

The members of vector 0 used in function f2i(0) are:

S -
0k 31.81.. for k-III+1,....,M+N (4025)

Considering the above definitions, the state equation matrices

 

 

 

 

are:
F , 0' F 0 ‘
° §N(N-l) °
-1 _ -l
0 -92 BNN 0 ' 9. " 1.“. ZNN 0
B
5 = o o -3NN o —NN
and
F 1
9(N-1)m °[
3 = 0 0
0 ENNJ (4.26)

 

 

97

 

where:
ENN = diag (Mi) (4.27.a)
QNN a diag (Di) (4.27.b)
ENN = diag (ai) ' (4.27.c)
ENN = diag (Bi) (4.27.d)
and
§N(N-1) = 1 '
-1(N-1) (4.28.a)
-£(N-l)(N-1)J
9(N-1)m =

I
[—(N—l)(N-l) -T(N-l)(m-N+1) ] (4.28.b)

The matrix ENm is of the following form:

INm = [$1 32 33,...., gNm] (4.29)

2i is an Nx(Nfl) matrix of the form:

P 0 -1
-(i-l)(N-i)
1 1 ............... 1
.2 -1 0 ............... 0
-1
0 -1 ............... 0
b o .............. -1 (4,30)

 

 

98

 

and:
11(N_1) 01(m-N+1)
INm ’ 'I(N-1)(N-1) 5(N-l)(m-N+1)
‘ (4.31)
and:
INm ‘ 5NIN-1) 9-(N-1)m ”'32)

The system equations are now in the form of equations
(4.1) and (4.2) and hence the conditions (i) and (ii) and

(iii) of Popov criteria could be Checked at this point.
Let:

1m 52

9NN (4.33)

lo )QIH

H:

With this choice of H the inequality of (4.3) is in the

following simpler form:

f1(ok) Gk 3’0 for all k=l,2,....,m (4.44)

defines the region of stability 5.

The above inequality is satisfied for the range of:

Ominkf‘ok: (0 or ask) and(.0 or'os):0 :0!“an

where

omink = -"-(5:j+61j) , ask = agj—aij

and

qmaxk=: - 5§j+éij and aij = Sin -1 (3:3? Sin sij/Eigj)

which will be discussed further in Chapter 5 when the region
of stability is more thoroughly discussed.
Considering the transfer function for the linear part

of the system:

99

sees-14‘ 2 >' (.4 '_r_ 0
fl(s) = g’(§1-§) g = o 5(sx+a)’l e
= 31(3) 0
0 fl2(s) (4.34)

The matrices land _Q_should be found satisfying the
conditions Moore-Anderson theorem including the constraint
(4.3) of the Popov criterion. A possible choice of 5 could
have been the zero matrix but this would have caused 3T9? to
have a pole zero cancellation with 3(3) which has a pole at
s-O, this justifies the choice of §_as in (4.33). The next
step is the Choice of V1(O) in (iii) of the Popov criteria,

and consequently matrix Q.

V1(o) is chosen as:

m 0k l N N s
v o — f f 0 do = - B.. 3.8. C 6..-C 6..
1( ) i=1 0 1( ) 2 i=1 §=1 13( 1 3( OS 13 OS 13)
- (6 -GS )ESES Sin 63 )- 1 § ) (E -Es)(E -ES)B Cos 6?.)
ij ij i j ij 2 i=1 i=1 i i j j ij 1]
3541
(4.35)

It is obvious that V1(0) is positive just for a range
of 3 around gfg. This condition is very practical in the
stable case because the angle separation of the lines
remains almost the same and in the stable case there is no

voltage collapse or large increase.

100

After the choice of V1(0) matric_Q_could be chosen to

satisfy (4.4), that is:

V1(g_) = Q’_f(_q) for all gen” (4.4)

The components of §V1(O) are:

321—: B (E E Sin (o +53 )-ESESS 53 .) 4
dok ij i j k ij 1 J in lj ( .36)
for k=l,2,....,m

and
dvi 1;

8 - B. E. -E. C 5. .+ - .. =
33; j= 1 ij(J 8) os ij B. lej (C036S lj Cosélj)

#1
If s s
B .(E. Cos 6..-E. Cos 6..) for k=m+l,....,m+N

j=1 1] 3 13 J 13 (4.37)
#1

It could be concluded that [6):
V ..

and the value of‘Q by comparing (4.4) and (4.38) is:

9- ’ E-Im+N)(m+N) . (4-39)

Having found the matrices g and Q for the Moore-Anderson
theorem, the next step is verifying the positive realness of
matrix Z(s), defined in (4.6) as:

5(8) = (flst)fl(s) (4.6)

Substituting for g and Q and fl(s) in (4.6):

101

l -1 -1

<<§+s>rtsifm .4 2) 0

0 s[slfg]-l 5

1(8)

EIIS) 0
0 Ԥ2(s) (4.40)

Since 5(3) is block diagonal with diagonal blocks 11(5)
and E2(s), each one of them can be looked at separately.

It is obvious that conditions (4.6.a) and (4.6.b) hold
for 51 and 32. Condition (4.6.c) holds under the condition

q>mi/d and Bi>o for i=1,2,....,N.[6]
i

4.3 Construction of the Lyapunov Function

It could be shown that it is possible to find matrices
H and Q to satisfy the Moore-Anderson theorem condition such
that a Lyapunov function of the form:»

V(X) = -;-x' g x + vl(g) (4.7)

can be constructed, where P is a positive definite symetric

matrix satisfying 4.8 and 4.9 and 4.10 and has the form:

31 0

g a 0 32 (4.41.a)

L and we in 4.8 and 4.9 and 4.10 are auxiliary
matrices.

Considering the fact that Z(s) is positive real and
recalling Lemmas 4-A-3 and 4-A-4 and results of 4-A-1 from

appendix 4-A, 2(3) + Z(-s) can be factorized as:

Z(s)+Z'(-s) = Y'(-s)Y(s) (4.41.b)

 

102

where Y(s) has the minimum realization (A,B,L) i.e.

1

Y(s)-Y(m) = L'(SI-A)' B (4.41.c)

and

W0 = we) (4.41.d)

The matrix fll(s) in (4.34) can also be written as:

_ ’ -1
where:
0 K' o g-
A a o -M‘1D B . M-lT d = o
_1 .. _. I ..1 _ an E.

(4.43)

It is trivial that 912180 and 2412130 and so from (4.8)
and (4.9) and (4.10) the following matrix equations hold for

21:

931.121 + 2151 = -_1_I_._'1 _<_ o (4.44)

9—121 3 9.1931 + 5'1919'1 (4.45)

Expanding (4.44):

 

 

 

P11 312 ° 35' ° 0 P11 312
-1 + -1 =
2.21 2.22 0 ‘5. E E. '31 9. P21 322
ro P K'-P M-lD 1 'L ‘
—11— -12— - -11
= - [L' L' ]
_ -1 _ -1 _ -1 —-11 _ 12
K311 5. 2321 £21§'+§£12 322$ P. 5 £22 £112
a _ E11 E'11 £11 9312
£12 11'11 £12 9'12 ”'45)

Considering (4.46), the following equations are found

to hold:

L.

L 11

—11

and:

2.11 E."

2123' +

103

0 => 211 = o (4.47)
MID = L' a o
-12- - —11— 12
-1 -1 _ _
E 312’2225 P. " L“. 2 E22 ’ E12E'12 (4-43)

Evaluating the component matrices in equation (4.45)

"-1
E11 E12 .. .T.
P a a P “'1 M 1 T (4 49)
—l-l —12 -22 —22 — - '
§_(N-l)m l/q _I_mm 0 a l/q-g
gag} 0 0 0 0 (4.50)
and using the property of E §_= 3
0 0 g_ I 0 0 0
-1 8 3
4321931 5. :4 2 0 0 .I. 5.9 2 WM)
hence the equation (4.45) becomes:
-1
2.125 1 E/q
-1 a
and so:
-1 _ l
P M'1 T = T (4 54)
—22- — — '

While P22 is
expressed as:

-1
(322.!

NxN and symetric and T is Nxm (4.54) can be

)_'r_=o

-Nxm (4.55)

-NxN

The solution £22 is obtained by multiplying (4.55) by gpto

obtain

104

Z(El‘lxm = 9Nxm (4'56)
and then solving it for Y =

1 - g [3221fl (4.57)

2qu]
Y is an unknown symmetric matrix (NxN). P22 is then
obtained knowing the solution for Y. The procedure for
determining Y is now discussed.

Considering the structure of matrix 3, each column
contains only two nonzero elements, 1 and -1. This shows
that all elements on the same row of Y are equal, and since
Y is symmetric it follows that all elements of Y are equal.
Hence a necessary and sufficient condition for a symmetric
matrix Y to be a solution of (4.56) is that it has the form

YauU, where )1 is'a scalar constant and U is an NxN matrix

with all its elements equal to 1.

Now 5 gezfl-lfg, is symmetric so from the above discus-
sion
3 -1- =
X 53225 5 no (4.58)
and:
-1
E2222! §+ 112 (4.59)
-1
£222 11!! 2+1. (4.50)
-1
322 3&3“! Q! (4.61)
Since 312 is not symmetric consider (4.48):
-1 a
£115. ‘ 2.125 E 9. (4.62)
and
-l
5515' 12.13.1214 2=9. (4.63)

Since 5 £11 13' is symmetric, then 5 _P_12 5'12 is symmetric

from (4.63). From (4.53), it is apparent that

-1 _ l _ 1
hence:
-l l -1
Matrix 5 glzg‘lg-égis symmetric since 5 3125-12 is symme-
tric. Multiplying both sides by 251:
-l -1 l -1 , a
and following the same argument as before
_¥_ 2 (I_(_ 2125 Q qD)D 02 (4.67)
and thus:
1
531235-‘+ 9.9.! (4.68)

-l

q

The solution for all could therefore be concluded as:
..1.

£125 2 q . 2+ 9232 (4.69)

5 £113.. ' K

In the above equations u and D are nonnegative scalars
and g is an NxN matrix with all elements equal to 1 and
satisfying the following inequality:

2(D-M/q) + (u-p) (D U M + M U D) 3_0 (4.70)
The inequality (4.70) is originally obtained from an assures
satisfaction of the following inequality:

£151 + 5'31 _<_ o (4.71)
Considering the right hand side of (4.46) and recalling that

L - 0, (4.71) would reduce to:

11

-L 0 (4.72)

12L'12 5-
Now by comparing the left hand side of (4.72) and right hand

side of (4.46) it could be concluded that the following

inequality should hold:

-1 -1
321‘“ +5312 "2.22! 9.7" 02225.0 (4.73)

106

Since the element 32110 has not been defined yet at this
point it is necessary to find the solution for £215"

From (4.46) recall this equation:

_ -l
52.11 ‘5 9.2.21

multiplying (4.74) from the right by K' and taking the

3 o (4.74)

second term to the right hand side:
-1

52115.. = fl 2.2.213... (4.75)
Substituting the result of (4.69) in (4.75):

l -l

a-2+DQ!Q'E 9.3215. (4°76)

multiplying (4.76) from left by 2‘15:

5.61- + g g Q =- 221 5' (4.77)
now substituting (4.77) and (4.68) and (4.54) in (4.73):

Z(Q-Q/Q) + (Ll—'0) (2214 + g g Q) g 0 (4.70)

The following discussion gives the condition for (4.70)
to hold when it is equal to zero.

Let u* - u-p, and consider:

2(2-§+u* (29§+§92) =0 (4.78)
then

.1... _M -1 =_

2n ‘2 a) (2.2!!1‘EEE) in“, (4.79)

Since 11* is a scalar and the maximum rank of Q g 11 +
91 U Q is 2, consider all the first and second order
equations in terms of 11*. The matrix i— [p_ - %]-1[ ]_.')_ g E +

g g Q] has the form:

107

 

 

 

 

 

 

 

F
201141 D2M1+M2Dl ...... MnDl+Dan
D.- . - -
1 El) 2101 M16 ) 2(Dl Ml/ )
1 -1 q
7(0-5) (DUM+MUD)= q
q
M1D2+D1M2 2132M2 ..... . M D2+M2Dn
2(1)2 M2/ ) Z(DZ-Ei) 2(02-53)
f q q : q
M ............. 2.4:. Lu
2(Dngg) (ON-EH)
. q q
(4.80)

Since the rank of equation (4.80) is 2 there are only
two sets of equations in terms of 11* that should be con-
sidered in (4.79). The first set are obtained by setting
each diagonal element of the r.h.s. matrix in (4.80) to 1.
These equations are

*
Dl"1" + 1 = 0

ll
0

DNMNu* + 1

D -M
N Nél

The other set are the second order equations that are
the determinant of block diagonal 2x2 matrices in the r.h.s.

matrix of (4.80). These N-l equations are of the form:

 

108

 

 

 

2 2
2 (D M -D M) *
.. 2 2 2 2 1 2 2 l )1
u - - - =
(4D102M1M2 ZMleDlDZ 02 M1 M2 D1 )_1 4(D -M )(D _M
4(1) -M)(D-M) l—l 2"2
1 1 2 2 q q
q q
2 2
-(D M -D M ) u*
2 3 3 2 -1=o

q- . Ei-
(N-l) terms

Adding these equations the following quadralic equation in

i* is obtained

 

 

N‘l N 2 N D.M.
_u*2 Z Z (DiMi DiMi) "’ (N’1)+u* Z ———Df_:l.+N=0
i=1 j=i+l 1(Di-lg.) (Dj-Mj/ ) 1:1 1 _i
E_ q q
(4.83)
N-l N (D.M.-D.Mi)2 * N DiMi
u*2 ITi—l——%T—— -u 5 D.dM -130
_— q
q q (4.34)

If u* lies between two roots of the solution of the

quadratic equation, the above inequality holds.

When the damping torques are uniform u* reduces to

“o where:

q) (4.85)

-1:

0

109

Having found 21, the following is a description of

determining g2 of matrix 2.

Recalling the transfer function W2(s) from (4.34)

-1
912(3) 3 Q'2(S_I_"A2) 5.2 (4.86)
where:
A2 a ..QNN’ B2 ’ éNN' Q2 ’ lNN (4°36°a)

Since §2(s) is positive real, using Lemma 4-A-3, 22(8) +
2'2(-s) can be factorized as

22(5) + Z'2(-s) = 1'2(-s),¥_2(5) =

s/IB. s/IB.

 

 

s(s£fa)-l §|+ S‘Sl11)-1 §'= diag (s+ail)' diag (s-ail
and so: '
3/78. (4.87)
Y2(S) a diag (s+ail )

Considering Lemma 4-A-4[14]. Y2(s) has the minimal reali-

zation as following:

-1
solving (4.88) for £2:
L2 = -diag (ff—ai/ ) (4.89)
B.
1
But L2 and P2 are related as:
3252 + £222 = ‘EZE'Z (4.90)

By solving (4.90) with (4.86.a) and (4.89), P2 is determined
as:
-1
P2 = a8 (4.91)
The Lyapunov energy function can now be found by

substituting £1 and 2 into (4.7):

2

llO

 

 

 

 

311 E12 0 Er
V(X)=%[§'r'9-"AE']' 15’21 P22 0 Q + V1(9_) =
0 0 P A
. ‘Zj _ j
1 . 1 ,
76 (9/ +7 292): +57 (M/q+PBQ_l~_1)_Q
+ 1 ' 1 I ‘1
2 9’. (E + u! 2 EM». + 24.5.1 _ ii éfi + V1(o) (4.92)

The term 1/q was introduced initially to avoid pole and
zero cancellation in (§+Q_s)§_(s). Now that the energy
function is obtained, setting q -—-°° will cause all the terms
that are factors of 1/q to vanish. Substituting for V1(0),
from (4.35) and using the results of Appendix 4.8 in (4.92)

will result:

VE=T——l 1% £1 MM (w-w)2+l(p*-u)(§ Mm)2
gmi i=1j=113 13 7 "i=11L1
i=1
N N

l s 2 1 s 2
+ p 2 (D.(5.-6.) + M.w.) + 2 (E.-E.) _
7 [i=1 J. i 1 1. 1:} 71:1 1 1 /(x3 x
N N
l s s s
z z - . -
+ 7 i=1 j=l Bi]. EiEj (Cos Sij Cos 613” .2- (sij 51].) EisEj Sin 6.1].

1 ,
-§- (Ei -E. Sj)(E ~33.) Cos 65d] . (4.93)

111
l

where “o = - XE:

If damping torques are uniform, )1* equals “0 and the
second term in (4.93) vanishes. ()could be chosen zero to
simplify the form of energy function. Also if damping
torques are un1form or zero and with; mo = ZMiad/ the

energy function is: 2 i

N N N
_. 1 2 1
v1: — 51) M. (cu. m) + 2 g: 32 Bij((EiEj(Cos<5:j-Coséi.)

1 =1 3

_ _S SS . S_ _s _s
(aij sij) EiEjBij S1n aij (Ei Ei)(Ej Ej) Cos dij)

(Ei “E 32)

l//(xdiw

d)i (4.94)

The only task left in this derivation is to show that 3;:

for VB defined in (4.94) is negative.

Considering equations (4.12) and (4.13) and (4.15) and

(4.16):
N N
dVE s .
= Z 2' B. E. ES 6 - .E. s 5.. w.-w.

HE— i=1 j=1 ij(l j Sin ij El 3 1n 1]) ( 1 J)

N N SE S S
- Z 2 B E. ' 6.. - E.E. S' 5.. w.-w.

1-1 j_ _1 1j(E i 3 Sin 13 1 J 1“ lJ)( 1 3)

N dEi 63 5
+ Z Z - E. ..

i=1(a——) j= l Bij(EjC°S ij 3 Cos 13)

N dB.

2

- z ( 5/dt) _

1-1 /(Xdi x'di) (4.95)

112

Equation (4.95) has been proven to be true in Appendix 4.C.

Hence:
dVE Iq 2
..-.--X T’ ~(dE. ) \-
a—— . d _ .
t, 1:1 01 b/dt /(xdi x di) (4.96)

Since the right hand side of (4.96) is always negative, VB
is derived Lyapunov function for the specific power system.
Now that the Lyapunov energy function is derived, it is
time to drop the ficticious generators connected to the load
buses. This could be done by setting €—-- 0. The

differential equations become algebraic equations

N

s s . s .
2 3.3.3.. 5 5.. - 3.3.3.. S1n 5.. = o 4.97)
j=1 1 3 13 in 13 1 3 13 13 (
N s s

B.. E. Cos 6.. - E. Cos 6.. = 0
§=1 13( 3 13 3 13) (4.98)
and

s
31:31

where B? can be a function of time. Since

N s s . s . 4 9
Pmi = €31 EiEjBij Sin 6ij 1-l,....Ng ( .9 )
N 5 S s ' 4 100
N (4 101)
I . = 2 B..E. Cos 6.. .
d]. j=l 1] J 1]

These algebraic equations specify a constant power real
power load power from (4.97) and (4.100) and a constant

current reactive load model since from (4.98) and (4.101)

8-
Idi - Idi

113

Since the transient reactance of the fictitious generators
was omitted from the model entirely, these terms do not
appear in the energy function, letting E -- 0 the inertias

are

M1 = 0 for 1=Ng+1: Ng+2,....,N

and the susceptances

xdi xdi for i=Ng+1. Ng+2,....,N

N N N
1 9 2 l 3
VB = - Z M.(m.-m ) + — 2 Z B..(E.E.(Cos 6 .-Cos 6..)
2 1:1 1 1 o 2 1:1 j=l 13) 1 3 13 1]
s s s . s _ s _ s s
- (513 Gij)EiE Sln 5ij) (Ei Ei)(Ej Ej)Cos 5ijfl
s 2
1 29 (Bi-Bi)
+ 7 ——+——T—

(4.102)

The energy function matches the energy function
developed in (3.35) if the flux linkage term (last term is
in 4.96) and the real power load term is assumed constant in

(3.35). This energy function can be written as

114

 

8
VB = Z W
i=1 1
where
N
w = 1 1g M (m -m )2 (4 103)
i 7 . i i 0 °
1=l
N
g 2 2 32 s2 s s s 3
W2 - iglkEi +Vi -2EiVi Cos (61-81))-(Ei + Vi -ZEiVi Cos (61-61))
_1_
2x di - (4.104)
N N ‘
1 s s s s
W = - — Z Z B..(V.V. Cos (8.-6.)-V.V. Cos (8.-8.))
3 2 i=N +1 j=N +1 13 1 3 1 3 1 3 1 3
9 9 (4.105)
N5 N9 Riv:
g _ _ s S _ _ _ S —T—" o S_ S .
W4 2. (6 61) Pmi E. (<51 61) xzdi Sin (61 61) (4.106)
1—1 —1
N N N
. s s
w = -z (e -e?) 913 = -z (e.-e?) ((2 v?v?s.. Sin (9 -e.))
5 1=Ng+l 1 1-1 1 1 i=1 1 3 13 1 3
s s
E V.
i 1 . s s
+ i-rd—T- Sln (Bi-61)) (4.107)
1
PB 32 32 s s s 3
N6 = +1_ -((Ei +vi -2EiVi Cos (61-81))
1—1
s s s s s s l
N N V
W7 = + 2 z Viv; Bij Cos (GE-8:) (—§ - 1)
1=Ng+1 3=Ng+1 Vi (4.109)
s 2
N (E.—E )
w _ 1 £9 1 _11 (4.110)
8 - 5 1:1 ’51 xd1

115

The terms Wi can be associated with the following
1) W represents the kinetic energy and W repre-
3 nts the magnetic energy stored in the
transient reactances of generators.

2) W3 represents the magnetic energy stored in
transmission lines and transformers.

3) W represents the position energy due to angle
displacement of generators.

4) W represents the position energy stored in
tRe real power of the loads.

5) W is the energy due to the reactive power
pébduced by the generator and the energy
stored in the transient reactance.

6) W is the magnetic energy stored in the
rZactive loads.

7) W is the form of kinetic energy and is due to
consideration of flux linkages.

The terms W2, W3, W4, W5 and W6 and W7 are all the
different components of the potential energy of the
generator. 1

The procedure followed in the development of this
Lyapunov energy function could be utilized to produce a
topological Lyapunov function for the case where

l) the synchronous generator model is more
complex and includes amortisseur effects:

2) excitation system models are included;
3) power system stabilizer models are included:

4) governor turbine energy system models are
included

as long as a constant current reactive and constant power
real power load models are assumed. The development of
Lyapunov functions for the case where general real and

reactive load models with constant impedance, constant

116
current, and constant power components would require the
development of a different construction procedure than were

used in this Chapter.

4.4 Simulation Results

The Lyapunov energy function derived as (4.97) is
tested on the 39-bus, 10 generator system of New England.
The schematic of this network is given in Figure 4. The
tests have been performed on two fault cases. The fault
cases considered is a fault on line 28-29 and a fault on
line 14-23 where in each case the line is removed to clear

the fault.

1) Fault on Line 28-29 -

The fault occurs on line 28-29 at t=0 and the critical
clearing time has been determined by iterative simulation
runs to be .21 seconds. Figures (l3.a) and (14.a) show the
variations of angles of generators which are close to
faulted line. As can be observed in (14.a) generator 9 is
the first generator that loses synchronism and thus would be
called the critical generator for this fault.

Figures (13.b) and (14.b) show the total and potential
and kinetic energies for both stable and unstable cases. In
both cases the kinetic energy at the time that fault is
cleared is very large and the potential energy is very
small. The kinetic energy gradually decreases and the
potential energy increases, but in the unstable case when

the generators lose synchronism the kinetic energy has a

117
sharp jump and continues increasing due to loss of synchro-
nism. The total energy in both the stable and unstable
cases is nonconstant because of the damping caused by
generator voltage variation.

Figures (13.c) and (14.c) give the three different
components of potential energy. It shows that in the stable
case energy due to generator angle displacement decreases
until a minimum is reached and then it starts increasing
again, but in the unstable case it is always decreasing.
The fact that generator angles continue increasing with
respect to the stable equilibrium point justifies this
behavior. The energies due to real load power and
transmission grid in the stable case have a maximum and then
decrease. This is caused because the load bus angles and
voltages oscillate around stable equilibrium point. In the
unstable case, the continuous increase in angles cause these
energies to increase for a time and later drop when the loss
of synchronism occurs (when bus angle differences across
particular branches exceed 360°). The flux linkage and
reactive load energies are shown in Figures (13.c) and
(l4.c) for both stable and unstable cases. These energies

are relatively small compared to other forms of energies.

ii) Fault on Line 14-23 -
The results of this test are the same as case (i). The

angles of generators and different components of energy for

stable and unstable case are shown in Figures (15) and (16).

118

Fig. (13.a) Variation of generator angles for the stable
case where the fault occurs on line 28-29, the
fault is cleared at 0.2 seconds and line 28-29
is removed.

“1.
9t

“4

CL

.4

0001.0

BNMJEOEGHI

 

 

r f fl

h " 10 u
TIME IN SECONDS

119

Fig. (13.b) Variation of potential, kinetic and total

ENERGY P.--U .

1S0.

 

—

energies for the stable case where the fault
occurs on line 28-29, the fault is cleared at
0.2 seconds and line 28-29 is removed.

 

 

TIME IN SECONDS

Fig.

ENEROY (P M)

120

(13.c) Comparison of different components of potential

_3

energy for the stable case when fault occurs on
line 28-29, the fault is cleared at 0.2 seconds
and line 28-29 is removed.

“WIN”

NW WWI-diam...)

 

1'.

1'1

 

1 ' i ' ‘
TIME IN SECONDS

0!". M! 0“!“ tall"
r

121

Fig. (13.d) Variation of flux linkage and reactive load
energies for the stable case when the fault
occurs on line 28-29, the fault is cleared at
0.2 seconds and line 28-29 is removed.

euenoy P.U.

'0‘“ Who. m
4",—

 

_r

' j\\ \/ 57
_“ TIME IN SECOND!

 

122

Fig. (14.a) Variations of generators angles for unstable
case when the fault occurs on line 28-29, the
fault is cleared at 0.21 seconds and line 28-29
is removed.

(m

11!!

104+

RNOLE DEGREE

 

 

 

i .'
TIME IN SECONDS

123

Fig. (14.b) Variation of potential, kinetic and total energy
for the unstable case when fault occurs on line

28-29, the fault is cleared at 0.21 seconds and
line 28-29 is removed.

154

¢mm m x: an

L’ «*5'
V l
\ i '
WM "'0'
TIME IN saconos \///

 

 

”1

 

ENERGY 9.0.
E

124

Fig. (14.c) Comparison of different components of potential
energy for the unstable case when fault occurs
on line 28-29, the fault is cleared at 0.21
seconds and line 28-29 is removed.

)

 
       

Ill DAD m
\

 

WWII“:

..ir” ---- i ,__ 3-
TIME IN sccouos

 

 

      
 

MAM man MAC”! (Him

ENERGY (9.0)

 

125

Fig. (14.d) Variation of flux linkage and reactive load
energies for the unstable case when the fault
occurs on line 28-29, the fault is cleared at
0.21 seconds and line 28-29 is removed.

IL

 

ENERGY P.U.

   

RU! IN”! m

. ("0:30. , s
.21 ' I \‘2
TIME IN SECONDS
O

 

 

126

Fig. (15.a) Variation of generator angles for stable case
when the fault occurs on line 14-33, the fault

is cleared at 0.23 seconds and line 14-33 is
removed.

105
and

75.

45.

.... /

TIHE IN SECONDS

RNGLE DEGREE

 

 

 

-IS.

127

Fig. (15.b) Variation of potential, kinetic and total
energies for the stable case when the fault
occurs on line 14-33, the fault is cleared at
0.23 seconds and line 14-33 is removed.

 

=3 l“. mam-um
L |
\
» O
O I
m o
W I
“z" n
“4 I
I
: 'Onm‘l [Hm
I
o
H °
. thud'

 

1 IN
TIME IN SECONDS

128

Fig. (15.C) Comparison of different components of potential
energy for the stable case when fault occurs on

line 14-33, the fault is cleared at 0.23 seconds
and line 14-33 is removed.

 

     
   

. lul m m
100‘
u.___..~.1=f. A
:F i

"‘4 ' TIME IN SECONDS
5?
Q. ' ”.4
H annulus minim."
S x
m
g s.
w .

n:

 

129

Fig. (15.d) Variation of flux linkage and reactive load
energies for the stable case when the fault
occurs on line 14-33, the fault is cleared at
0.23 seconds and line 14-33 is removed.

 

4
w
I ' * /
1 .n . l [J
-41 ~\\\\\\\\\\\\“‘~_—_—_ TIME IN sscouos

IKYM W “'01

ENERGY (P .U)

130

Fig. (16.a) Variation of generators angles for the unstable
case when the fault occurs on line 28-29, the

fault is cleared at .24 seconds and line 14-33
is removed.

 

-5 i / i

3
1 INE IN SECONDS

 

Fig.

ENERGY P.U.

131

(16.b) Variaton of potential, kinetic and total

2‘.)

 

energies for the unstable case when the fault
occurs on Line 14-33, the fault is cleared at
0.24 seconds and line 14-33 is removed.

Iran. lfllldl
‘,.~—"JL _“fl"””/,//’zr~—~\\

'OYINIIN OW

/'__uuncy

J‘ r i 1, I ‘ . 4_fi
1..
TIME IN SECONDS

   
     

132

Fig. (16.c) Comparison of different components of potential
energy for the unstable case when the fault
occurs on line 14-33, the fault is cleared at
0.24 seconds and line 14-33 is removed.

 

 
 
   
 

 

m
. m mu mm
nu ‘1
3». mm mum um

 

} 4
4%
4 i 1

 

 
 

. "t TIME IN SECONDS

,3; -...i

2 I ‘ "In.” IINWENEBY
'i' J” -

U 1

 

--m.)

Fig.

ENERGY P .U .

.0‘

-II

 

‘28

133

(16.d) Variation of flux linkage and reactive load

energies for unstable case when the fault occurs
on line 14-33, the fault is cleared at 0.24
seconds and line 14-33 is removed.

 

TIDE IN SECONDS

 
    

llACTM DAD OHIO?

134

APPENDIX 4‘. A

4.A.l) Minimum realization of matrix Z(s) in (4-6) - expanding

Z(s) with respect to A,B, and C:

Z(S) = (H+Qs)W(S) =

l l

B+QC’((SI-A)+A)(SI-A)- B

= QéB+(Né+QéA)(SI-A)-IB

NcIsI-A)'

if Z+m Z(w) = QcB so:

Z(s)-2(a) = (Nc+QcA)(s1-A)'1B

thus the minimum realization of Z(s)-Z(m) is:

(AIB'CN’+ A’CQ’)

Lemma 4 .A. 2 [14]

Let Z(s) be a matrix of rational transfer functions
such that 26”) is finite and 2(3) has poles which lie
in ReS<0, or are simple on ReS=0. Let (A,B,C) be a
minimal realization of Z(s)-Z(m). Then Z(s) is posi-
tive real if and only if there exists a symetric
positive definite matrix P and matrices W0 and L such

that:
PA + A’P = -LL’
PB = C-LWO
wgwo = 2(a) + z’(w)

Lemma 4.A.3 [14]

Let the nxn matrix Z(s) be positive real and suppose
that 2(3) + Z’(-s) has rank r almost everywhere, then
there exists an rxn matrix Y(s) such that

2(3) + Z’(-s) = Y’(-s) Y(s)
Lemma 4.A.4 [14]

Let Z(s) have a minimal realization (A,B,C) and Z and
Y be related as in Lemma 4,A,4,Then there exists a
matrix I.such that (A,B,L) is a minimum realization
for Y.

135
APPENDIX 4 . B

In this Appendix an expansion for each term in (4.92)

is presented:

4.8.1 Term w’P w

22

The solution obtained for P22 in Chapter 4 is:

P22 = M + )JMUM (B-l-l)
Substituting for M and U

r1414101412 uMle ........ leMN -1

P22 8 uN}M2 .. ........... ~ ........... MMZMN
3 2
“M1MN ...... . ................. MN+uMN (B-l-Z)

 

 

Substituting for P22 in m’Pzzw the following form will be
obtained:

. _ 2 2

szzw-uml ml +M1M2mlw2+MlM3wlw3+MlM2wlw2

2 2

2 2
+112 “2 +M2M3“’2‘”3+M1M3‘”1“’3+M2M3‘12“’3+M3 “3 +

+ .......... +MN2wN2)+(Mlm12+M2m22+ ..... +MNwN2)

(B-1-3)
But u=p*+p, substituting this value for u in (B-l-3) and

rearranging it.

136

, y 2 2 2 2 .... 2 2

2 2 2
+ 2M1M3”1”3+ ‘ '+2MN-1MN(”N-1”N’+(M1”1 +M2”2 +°"'+MNwN ’

2 2 2 2 2 2
= * o oo o o o o.
u (Ml ml +M2 wz + +MN mN +2M1M2wlw2+2MlM3wlm3+

2 2
+2MN_1MNwN)+p(Ml ml + +2M _ M m )

2 2 2
+(M1wl +M2w2 + +MNwN ) (B-1-4)

Now considering the last term in (B-l-4) and calling it

 

 

 

F:
.. 2 2 ..... 2 _ _
F - Mlml +M2m2 + +MNmN (B 1 5)
N
Multiplying and dividing the r.h.s. of (B-l-S) by 2 Mi
i=1 :
N
Z M.
i=1 1
2 2 2 2 . 2 2 2 2
a) no. 0000‘
F: (“1 1 +M1M2‘”1 +' +M1MN‘”1 ’ + (Mlemz +142 ”2 + “2MN”N )
(M1+M2+ ' ' ' '+MN) (M1+M2+o o o o ‘+MN7
2 2 2
+ ..... +(MNleN + ..... +MN mN )
(M1+M2+----+MN) (B-l-6)
Rearranging (B-l-6) and considering the fact that uo=-N__—
Z M.
1:11
2 2 2 2 2 2 _ 2 2 2
F=‘"o‘M1 ”1 +M2 ”2 +’ """ +MN ”N ) (M1M2”1 +M1M2”2 +M1M3”1
2 2 2 w 2
+M1M3w3 +M1M2m1 +M1M2°’2 + """ +MN-1MN N )“o

(B-1-7)

137

 

N N
Now adding and subtracting u 2 2 M.M.w.m. to the r.h.s.
oi=1 j=l 1 3 1 3
of (B-1-7) will result: j¢i
N ZuoN N
F=-uo(2 Mim.)- Z Z M..M(O) -(O.2) (B-1-8)
i= 1 1 71 i=1 j= —1 1 3 1 3
Substituting (B-1-8) in (B-l-4):
N N u N N
l 1 2 l 2_ o 2
—wP w- *(2 M.m.) -—uo (Z M. m. ) —— Z Z M.M.(w.-m.) +F
2 22 -2“ i=1 1 1_1 i.i 4 i;l j=l 1 3 1 3 l
1 N N M M (m -m )2+1( *- )(g M w )2+F
=———4N z z 1313 2“ “o 1:111 1
2 M1 1‘1 3‘1 (3-1-9)
i=1
where F =p(M 26 2+---- +2 w ) (B-l-lO)
1 1 1 MN-1MN N
N
_ z M1”i
now let w = i=1
N
2 ”1
i=1
and consider
)1 N N - 2 2 .....
1§_£ 2 Mi Mj (mi wj )2 =7r(2M1M2ml +2M1M2w2 +
i=1 3= -11
+2MN-1MN”N 2‘”‘“‘111"12‘*’1 2 4M1M3”1” 3 '''' ‘4MN-1MN”N-1”N)
_2 N
l 2 2 2 2w — M.w.))
= (2M m +2M w + °°°°° +2 w + -4w(2 1 1
I 1 1 2 2 MN N N M_ i=1
1
i=1
_ 1 g M.(m.-6)2 (B—l-ll)
--§ 1 1
i=1

138
4-B-2) Term S’rpllar
The solution obtained for P11 in Chapter 4 is:
K Pll K’= D U D (B-2-l)
The matrix K is:
11(N-1)
'1(N-1)(N-1) (3-2-2)

The left inverse of matrix K is:

-1- -
KL ‘ [0(N-l)l I(N-1)(N-1)]

(B-2-3)
the right inverse of matrix K is:
-1 °1(N-1)
KR g -I
(N-l)(N-1) (B-2-4)

multiplying the right hand side of (B-Z-l) from the left by

Kgl and from the right by K-l, the first row and column of

matrix D U D will vanish and P11 will be obtained as:

2 .......
92 D203 D2DN

6 D ............ D 2 (3-2-5)

Now multiplying the both sides of P11 by 6’r and Sr:

, _ 2 _ s 2 2 _ s 2 .... 2 _ s 2
5: P115;:‘9‘02 (512 512) +D3 (513 513) + +DN (51N 51N)
+ 29 D (a -55 )(6 -5s )+ ----- +29 D (5 -a )-(5 -5S )
2 3 13 13 12 12 N-l N 1(N-l) 1(N-l) 1N 1N

(B-2-6)

139

Each term of the r.h.s. of (B-2-6) could be written as:

2 s

2
(512‘512

2

2

2- _ s 2 _ s 2_ _ s _ 3

But since 61 is the reference angle 51-6: could be assumed
to be zero and so by adding some zero terms to (B-2-6) the

form of 5’rPllér would be

’ _ 2 _ s 2 2 _ s 2 ..... s _ s 2
5rP226r— (D1 (5161) +132 (52 52) + +DN(6N SN)
5 s s s
+2D1D2(51-61) (52-52)+ ..... +ZD(N-1)DN(6N-6N) (6N-l-6N-l)

(B-2-8)

Now combining the F1 from (B-l-lO) and (B-2-8) and term

‘ G’rPlzw the third term of energy function will be:

N
%D(Z (Di(Gi-6:)+Miwi))2
i=1
The other terms of energy function of (4.93) could easily

be related to corresponding terms in (4.92).

140
APPENDIX 4 . C
In this appendix it is proved that the derivative of

total energy in (4.94) is negative definite.

If the first term is called Vk for kinetic energy:

aVk = 1 g MiM.(wi'w-) (c-1>
Bwi 2 N M .=1 3 3
z i 3
i=1
dm. N
. 1 _ 1_ s s . s _ . _
EE— - Mi gngij(EiEj Sin Gij EiEj S1n Gij) (C 2)
3V dw. N
k . 1 = l s s . s _ . _
‘5; EE‘ E'N' M §=1Mj31j‘313j 51“ 513 EjEi 51“ 51j)‘wi wj)
i=1 1 (c-3)
thus
N BV (1111. dV N N
k 1 k l s s . s . .
Z ——— - = -—— = z z B..(E.E. S1n 6..-S1n¢..)(w.-m.)
i=1 mi at dt 2 i=1 j=l 13 1 3 13 13 1 3

(C-4)

Now considering the potential energy terms in (4.94) and

calling them Vp:

3v -1 N s s. . s
ij jgl J J 3 3
d6.-
—_£1 = — C‘6
dt “1 “3 ( )

141

SO
av d6.. N
. 11 _ l X n(E Es Sin GS -E.E. Sin 6.-)(w.-m.)
53g. dt ‘ '7 j=1Bl j j l 3 13 l J
13
and
N .dv d6. N N
Z 332 ' dtil = -% Z 2 Bi ji(E HE Sin 65 .-E.E. Sin 61-)(w.
1=l ij i=1 j=11j1j l J 3
(C-7)
3V N
555 = §=lBij(Esi Cos 6. :j-Ei Cos dij) (C-8)

Equation (C-8) is obtained from second and fourth term in

(4.94) thus:
N N
g dEi . av2_ z dEi 2 Bi.(E: Cos 6:.-Ei Cos 51.)
i=1 EE‘ aEi j= -1 ‘E’ j=l 3 3 3
(C-9)

Now differentiating the last term in (4.94) with respect to

time and calling the last term V -

F'
s

dVF g g dEi . (Ei Bi)
dt 1:1 dt zxdi-xdi) (C-lO)
Substituting for (E.-E:) from (4.13)
dVF N dEi N dEi N

= _ z T’ ( —) _ , -Z 2 B. (E: Cos 68
at i -l Tdoi at /(xd Xdi) i=1 at j- -l ij1j
-Ei Cos éij) (C-ll)

Summing equations (C-4) and (C-7) and (C-9) and (C-11) will

give equation (4.95).

l

CHAPTER 5

 

STABILITY CRITERIA BASED ON
REGION OF STABILITY AND INSTABILITY

5.1 Introduction
The research presented in the previous two chapters has
for the first time provided the theoretical foundation
required to develop and theoretically justify stability
criteria for determining retention or loss of stability
given a particular fault, fault clearing time, and line
switching action required to clear the fault. Several
stability criteria and associated methods have been proposed
[18, 17, 7, 2] but could never truly be justified theore-
tically because Lyapunov and energy integral based energy
functions were never developed that
1) described the energy associated with both the
real and reactive current flows in every
branch, load bus and generator in the system;
2) determined the theoretical conditions imbedded
in the Popov stability criterion that permit
characterizing the region of stability and the
region of instability for a power system.
With this theoretical basis a discussion of the region
of stability based on the Popov stability criterion will be
given in Section 5.2. A definition of the loss of transient

stability will also be given. Theorems stating necessary

and sufficient conditions for the loss of transient

142

143

stability are established in Section 5.3. The conditions
stated in these theorems describe an estimate of the region
of instability. Although a precise boundary between the
region of stability and the region of instability is not
established and is a subject for future research, a
relationship is established between the region of stability
and the region of instability that allows a more precise
characterization of exactly what conditions are required for
retention and loss of stability.

In section 5.4, a cutset integral criterion is estab-
lished based on the theoretical description of the region of
stability and the region of instability. A brief discussion
of how this theoretical description of the region of
stability for this transient stability model could be
repeated for the classical transient stability model to
theoretically justify the P888 method based on the indivi-
dual machine energy function is given. .A theoretical
justification of the P838 method based on the outset energy
function developed frqm the topological energy function of
Bergen and Hill [7] is also given. This discussion
indicates that the methods of characterizing the region of
stability and instability either directly in terms of
constraints‘on the state variables or indirectly using the
integral criterion in this research relate to stability
criteria developed previously based on simpler transient

stability models.

144

5.2 Regions of Stability and Instability for Transient
Stability Models.

A loss of stability for the transient stability model
(4.12) can be observed by noting that the angles for one set
of internal generator and transmission network buses will
continue to increase with respect to the other buses (both
internal generator and network buses). ‘The angle
differences across the outset of branches that connect these
two mutually exclusive bus subsets will continue to increase
and ultimately exceed 360° indicating a loss of transient
stability has occurred. Although one could associate a pole
slippage and a loss of synchronism with an increase in
angles across all branches belonging to some cutset, one
might not always experience the continued increase in angles
beyond 360° normally associated with a loss of stability.
For the purpose of the theoretical developments in this
chapter a loss of transient stability will be defined as a
continuing increase in angles differences across branches in
a outset that ultimately lead at some point in time t** to

Gij(t**) > 360 + afij
for all branch pairs ij belonging to some cutset. The
cutset for which the angles exceed 360° is called the
critical outset and the bus pairs belonging to this cutset
is denoted ije J. It is assumed that this "loss of
synchronism”, that is called a loss of transient stability,
occurs on only one cutset for a particular fault. This
assumption may or may not be true and in no way effects the

theory to be presented.

145
.A more precise definition of loss of transient
stability will be given after defining what is meant by the
critical group and stationary group of buses. The critical
group and stationary group are characterized by

1) maximum bus angle differences of less than1A°
for all bus pairs ij in the stationary group;

2) maximum bus angle differences of less than A0
for all bus pairs in the critical group;

3) minimum bus angle differences for any bus ”i"
in the stationary group and any bgs ”j” in the
critical group must exceed 360°+6 j(t) where

1
"‘1’“ fiijm- 5ij(t)}_>_ 360°
E§E§ Sin 6:.

o _ . -1 1 J
513*“ - Sm [W]

J

A model relating the inertial center of the critical
group 5c(t) and stationary group 63a”) is now derived in
preparation for proving the theorems that are to be stated.
The voltage for each generator or load bus is represented by

Bi and the angle by 5i. There is a fictitious generator

connected to each of the load buses with inertia constant

Mi 3 Em. i=Ng+lpoooopN

1

and N9 generator buses with generator inertias M1 for

i=l,2,....,Ng. The N buses are reordered and then divided

into two different groups. A group of L buses is called the

 

critical group and N-L remaining (generator and load) buses

are called the stationary group.

146
Group 1 i=l,....,L critical group
Group 2 i=L+l,....,N stationary group

The following center of angles are defined for groups 1

and 2.
L
6 = Z M 6 L (5.1)
c k=l kk/z Mk
k=1
N
5 = 2 M 6 N (5.2)
53 k=L+1 k t//z Mk
k=L+1

and the following center of inertias are defined as

(5.3)

(5.4)

Since load buses have fictitious generators, a diffe-
rential equation can represent the variation in the bus
angle at every bus in the system.

N

M.6.=P.-Z E

1 1 m1 j=1 Sin (Si-Sj) (505)

iEjBij
where i=1,2,....,N

Considering (5.1) and (5.2):

 

N F2 Niki l N .
MsaGSa = i=L+lMi sfi=L+1 ' = Z Mi51
2 Mi i=L+l
i=L+l _

 

 

(5.6)

147

Substituting (5.5) in (5.6), one obtains:

"'6' N N N
= Z P '2 X EBOBooSi 6.
ea sa i=L+l mi i=L+1 j=l i 3 13 n 13
N L 88
+ z z s B B Sin 53.
i=L+1 5:1 1 j ii 13
N N 33
+ 2 z s B B 31 s
i=L+l j=L+1 1 j ij n Gij
N L
- Z 2 E E 8 Si 6
i=L+1 j=1 1 5 11 n 11
N N
-Z Z EBB Sin6.
1=L+1 ,=L.1 i i ii in <5-7)
N N
The terms that are 2 2 °

are zero because they are
i=L+1 j=L+1

the sum of real power in different directions.

Hence the equation (5.12) can be written as:

N L
s s . s
- 2 X E 8.8.. Sin 6..
sa sa i=L+1 j=1 i 3 13 1
N L
- Z 2 E EoBoo Sin6.. (508)
1=L+l j=1 i J 1] 13

Using the same procedure used to develop (5.8), it can be

shown that:

N
Z E.B.B.. Sin 6.. (5.9)
3

148

A formal definition of loss of stability for a
topological transient stability model is now given so that

l) a discussion of the region of stability can be
given, and;

2) the theorems stating necessary and sufficient
conditions for loss of stability can be stated
and proved.

DEF Given angles<5(t) for buses in the critical group
and angles 6}(t) for the stationary group, the
system described by (4.14, 4.15) is called
unstable if for t>t*

' - 6
5c(t) sa(t) ; 0
until at some time t**

61(t**) - 6j(t**) g 360 + 52j(t**) for all ijeJ

where

0 it = - ‘1 s s - 8 ** **
6ij(t ) Sin EiB' Sin Sij/Ei(t )Ej(t )

J

This definition of instability for a power system is

based on the Popov criterion condition (ii) that states

- - s (t)
s s . s _ s
[EiEjBij Sln aij] [6ij(t) aij] 3,0

Although this condition is assumed to hold for all Oij =
§j(t) - 6gj, it only holds for ranges

_ _ O S
n 6ij(t) g aij(t) g sij (5.10)

O O
6ij(t) g 6ij(t) g_n - 6ij(t)

149

if either Bi (t)Bj (t) < E SE: and Gij Z 0 or Ei(t)Ej(t) >

Ei “E and 5?. < 0

J 13
-1r -5°j(t) <6 jm§6°jm
s3
o. t < - 5.1].
61336 3-() 1r 6°j(t) ( )
if 31 (t)Ej (t) > E sEj and 5?j > 0 or Bi (t)Ej (t) < E iE? and
s
aijli 0
and
- fl - 6°. < j(t) < 1r- 6°. (5.12)
13 -1j
' . = 6.. . 6?. 6?. = 6? - 6?
1f El(t)Ej(t) E :Ej and 1J(t) 13 where 13 1 J
Gij the angle difference at buses i and j at
the post fault stable equilibrium point
ES voltage at bus i as the post fault stable

equilibrium point
Thus the system is not asymptotically stable in the large
but is only stable if the state 5TH) = (£T(t)£T(t)§_T(t))

satisfies 6(t)=gl§(t) eS(t) where

6 6 6 6 6

i = (512 613 IN 21 23°°°° 2N""" N-lN)

and S(t) is the angle subspace given by constraints (5.10-
*
5.12) for all tltc.

A more formal definition of the region of stability is

9(t*) = _C__:_c_ (t*) e S(t*) (5.13)
g(t*) = 0 (5.14)

E(t*) £5 (5.15)

150

for some times t*_>_tc where g? is the voltage magnitude at
the post fault equilibrium point. This is a very conser-
vative but mathematically precise manner of specifying the
region of stability because it constrains every element of
the state and places these constraints at a specific time.

The formal definition of the loss of transient
stability is consistent with the formal definition of the
region of stability because

1) the criterion that a loss of stability has

occurred requires that the angle exceed
360+6E t) for all ijeJ at some time t** where

the boundary of the region of stability is
51j(t):18°°'61j(t) at t* for all pairs ijeJ:

2) both are based on the Popov stability
criterion.

5.3 A Necessary and Sufficient Condition for Loss of

Transient Stability

The theorems will now be stated and proved to establish
necessary and sufficient conditions for loss of stability in
the power system. The conditions stated in the theorems can
be applied at any time tZtc+ and if satisfied indicate the
system trajectory for that fault, clearing time, and line
switching action will be unstable or if the system is

unstable indicate that these conditions had to hold at some

* ....
t _>_tC

Theorem 3
Given that

1) the maximum bus angle difference A between bus
pairs in the stationary and bus pairs in the
critical group is A=90°;

151

2) the angle differences across branches in the
critical cutset are all positive in the post
fault equilibrium state
6..(t) Z 6§11 o . all ijEJ

(5.16)

3) IE. (t) - E. s|<
then a sufficient condition for loss of transient stability

is

[Ei(t)Ej(t)Bij Sin 6 j(t)- -e. ”a a. Sin 6?.1161j<t)—6§j1 < 0

i j ij 1]
(5.17)
and
6 (t)-6 (t) > e > o i=1,2,....,L
c 3a ‘ 2 j=L+1,L+2,....,N
(5.18)
for some instant t*_>_tc where 5%. and E: are the bus angle

differences and voltage magnitude of the post fault stable
equilibrium point.

12:29.1

The differential equation describing the relative
motion of the equivalent machines representing the critical

and stationary group of buses is

6 6 1 .l 1' 1 6
- t=--+ Z B..E.tE.tS' --t
c(t) sa( ) (MC "daa;fl.j=L+1 13[ 1( ) J( ) in 13( )
- E. :83 Sin 6s. ] (5.19)
J 1j
Angle differences satisfy
6i.(t) - 63. 3 o .i=l,2,....,L
3 j j=L+1,L+2,....,N

by the definition of the stationary and critical group.

Thus, from condition (5.17)and (5.16), it is clear that

152

:04.

/////.~,.

.3

0:95 How cofiumufiuo >omom :3 .mflm

:

//

<_

 

153

. S
B 0 I u '-
lJ[81(t)BJ(t) Sln 6. .(t) B. faj Sin 6. :j] < o

i=1,2,....,L
j=L+1,L=2,....,N

. * -.. *
and 5i(t ) %(t ) Z 0.

Since 5 (t*)-5 (t*)>0 and 3 (t*)-5 (t*)>0 from 5u10
c sa - C sa ‘
6c(t)-6sa(t) will increase monotonically and exponentially

until for some t1 and some pair ijeJ

o
61j(t1) Z 360 + 6ij(t1) (5.20)
Although one branch ijEJ crosses the threshold (5.20)
not all branches must necessarily cross threshold (5.20) if

sufficient deceleration is experienced after time t1. A
branch causes accelerathma of the critical group based on

(5.19) over the range

180°-6§j(t) g 6. j(t) < 360° + 5°j(t t)

since

3 s s
B 1jEi (t)Ej (t) Sin 6i j(t) < BijBiBj Sin 6. 1j

and causes deceleration over

o+60. .. “6?.
360 i3(t) g 613(t) g 540 lJ(t)

since

. Es . s
BijEi(t)Ej(t) Sin 6i j(t) > BijagE j Sin 6ij
as can be observed from Figure 17.

This analysis is similar to that performed for the

equal area criterion based on Figure l.a in Chapter 2. Some

154

branches may cause deceleration on the critical group with
respect to the stationary group before a time t** when all
generators exceed the instability region threshold 360°+
62j(t). The basic question is whether the acceleration

energy due to each branch ij picked up as angle 6 (t)

. ij
increases from 180°-6Ci>j(t) to 360+6:j(t) will exceed the

deceleration energy caused by excursion of angle 6i (t)

beyond 360°+6:j(t). If the acceleration energy exceedsJ
the deceleration energy on all branches up to t** when all
branch ijEJ exceed y

6ij(t) 3 360° + 6§j(t) (5.21)
then a loss of stability will occur since all angles on all
branches ijEJ exceed (5.21) and (56(t)--3 “(112° for all t
satisfying t*$tgt**.

Since the bus angle differences in the stationary group

and the bus angle differences in the critical group do not

exceed A=90" then

max { 6.(t) }- min { §.(t)} 5’26 180°

ijeJ 13 ijeJ

Moreover, if all angles exceed 360+ 6?j(t) where branch iojo
satisfies
° ** - _ O **
m1n{ 6ij(t ) 360 6ij(t )}
1JeJ

= 6. . (t**) - 360 - 69 . (t**) = 0
030 1030

155

then the maximum bus angle difference satisfies

max 61.(t**) 5 6i . (t**) + 26 + 360°
ijEJ J 030

540° + 6i . (t**)
030

The energy A2ij defined by
540 + 69 . (t**)

1030
= . - S S .
360° + 6:j(t) (5°22)

is always less than

360° + 60j(t)

1
= B.. 3.3. Sin 6.. - E? 5 ' .. ..
Alij Jf 1]( 1 J 1] 1E] S1n 61])d613
180° — 6?.(t)
1]

(5.23)

assuming Bi is sufficently close to B? i=1,2,....,N from
observing Figure 17 because 6ij>0 ijeJ.

Thus, the net acceleration energy for each branch ij
will exceed the deceleration energy at t** so (5.21) is
satisfied. Thus, the conditions for the loss of transient
stability are satisfied and the theorem is proved.

The above theorem only concerns the case where

Gij Z 0
for all ijsJ on the critical cutset. Although this

restriction limits the generality of the result, a cutset

where some or all the angles 6:]. ijeJ are negative is

156

less vulnerable to a loss of synchronism than if Gij>0 since

the area A3ij defined by

800-st

A3ij = 6j:12:1j(ESE: Sin 6'1j - E. :3: Sin 6ij) d6ij
that represents the deceleration energy from the angle
across the branch at the post fault equilibrium point to the
boundary of the region of stability is larger for 61j<0 than
for'5:j>0. The cutset where loss of transient stability will
occur for a particular clearing time and line switching

action is thus much more likely to be a cutset where

51j>0 for all ijeJ than for a cutset where one or all

as
Gij < 0
for ijeJ. Although the case for 6?.<0 may be somewhat less

1]
important, the following theorem establishes a result for

some 6?.(0.
1]

Theorem 4
Given that

l) the maximum bus angle difference A for bus
pairs in the stationary and bus pairs in the
critical group is A =32 48°

2) some or all of the angle differences across
branches in the critical cutset can be
negative in the post fault equilibrium state

6Sj < 0 for some ijeJ

3)

[*1
II

5
Bi _ (5.24)

157

4) The angle difference 61 . (t**)=6: . satisfies
030 o o

5.. > 6. . = -A = -32.48° (5.25)

then a sufficient condition for loss of transient stability

is

. s s . s s
[Ei(t)Ej(t)Bi. S1n 6ij(t)-E.B.B Sin 6. lIGij(t)-Gi'] < 0

J 1 J 13' 13' J
(5.17)
and
6 (t)-6 (t) > o i=1,2,....,L
c 5’“ - j=L+1,L+1,....,N
s s (5.18)
for some instant t*2_tc where 61]. and B]. are the bus

angle differences and voltage magnitude of the post fault
stable equilibrium point.

212101

The proof is similar to that in Theorem 3. The
argument that 6c(t)- 63a(t)_>_0 and :sc(t)-:ssa(t)_>_0 for t"‘_<_.t_<_.t1
is identical to that in Theorem 3 where t1 is the time where
for a single pair ijEJ

o _ s
6ij(t) Z_360° + 6ij(t) — 360 + Gij

since 6?j(t)=6:j because 81(t)=B?.

The argument that the acceleration energy Alij on each
branch ijeJ is less than the deceleration energy is depen-
dent on A. In Theorem 3, A=90 but in this case A=32.48°.
All bus pairs ijEJ will satisfy (5.21) at t** when the bus

pair iojoeJ with the smallest angle difference satisfies.

6. . (t**) = 360 + a:

1.030 030

158

and the bus pair with the largest angle difference satisfies

6. . (t**) < 360 + 6? . + 2A
1030 — 1oJo
Given that 69 . (t**) = -A from (5.25)
1 J
o 0
then
360°+ZA-A=360+A
.. S S . _ . s
360°+6§.
1]
and
360°+A
S S . . s
Alij + AZij BijEiEj [Sin 6ij Sin 6ij] d6ij
180-6?
J

Since one de31res that Alij-l-AZiJ

(5.26)

(5.27)

.<0 for all ijeJ to assume

loss of transient stability. Since 6ij is negative and

6? . = -A < 6?.
1030 —’ 13
if
A . + A . . = 0
11030 21030
then
Alij + AZij g 0 for all ijeJ.

159

Bvaluat1ng A1i . + A2i j in terms of A,

030 o
360 + A
A . . + A . . s s 5 (Sin 6.. - Sin (-A))
11030 21030 - EiEjBij o 1]
180 + A
_ s s _ .
- EiEjBij ( 2 Cos A+nS1nA)
Ali ' + A
030 2i j is zero for all A=A where
o o o
_ 2 _ o . min _ s _ _ 0
tan AO — F" AO—32.48 . Thus 1f ijeJ -6iojo— 32.48 .

and if the bus angle differences within the stationary group
and within the critical group are within A=02.48° the
deceleration energy on every branch is less than the accele-
ration energy on every branch up to a time t** when
5 ** 53 --
ij(t )1; 360 + ij 13€J
Thus, a loss of transient stability is experienced.
Note that since
. O _ 0
min 51' = ‘A - -32.48
ijeJ 3'
and since the angle differences within the stationary group

and within the critical group is A,

S _
max 6ij —A

Although the minimum angle difference cross branches in
the critical cutset for stable equilibrium point might be
smaller than -32.48° and possibly the angle differences
within the stationary and within the critical group could be
larger than 32.48° and a proof of loss of transient

stability still might be obtained, the present result does

160

not appear overly restrictive and suggests that in order to

prove a loss of stability one must place a lower bound on

min 6?.
ijeJ J

and one must place an upper bound on the bus angle

differences within stationary and within the critical

group. A necessary condition for loss of transient

stability is now proved for the case where

1j 3 0 for all ijEJ.

Theorem 5
Given that

1) the maximum bus angle difference A for bus
pairs in the stationary and bus pairs in the
critical group is A=45°

2) that the angle differences across branches in
the critical cutset are all positive in the
post fault equilibrium state

6? 3 0 all ijeJ
S
3) Isi(t) - 81' 5, s

then a necessary condition for loss of transient stability

is that there exists a time t* such that

' S S . S S
[Ei(t)Ej(t)Bij Sln Gij(t)-EiEjBij Sln dij][6ij(t)-6ij] < o

(5.17)

and

6 (t)-6 (t) > o i=1,2,....,L
C 33 ‘ j=L+l,L+2,....,N
(5.13)

161

where 5:3. and B: are the bus angle differences and voltage

magnitude of the post fault stable equilibrium point.

Proof of Necessity
Given that for all t satisfying t*£t£t** and all ijeJ
6c(t)-6sa(t) ; o
and at t**

** 0 it --
61.j(t ) _>_ 360° + (Sij(t ) for all IJEJ

then from the definition of instability there is a t* for
which

° 6 -' *

6c(t ) 683“: ) > 0
and since the maximum spread of angle differences is 2A=90°

there is a time t* such that

o_ 3 t t - 0 t
180 61j(t ) S 6ij(t ) g_360 6ij(t )
for all ijeJ since

1) 6?.(t) must not vary by more than 90° over all
ij¥J in order for the load flow conditions to
converge:

2) 6..(t) must not vary by more than 2A=90° over
affl ijeJ by assumption.

Thus conditions (5.17) is satisfied at t* since all the
angles lie within the range where the inequality constraint
holds. Since both (5.17) and (5.18) have been shown to hold
at some t*, the theorem is proved.

It should be noted that the above theorems do not offer
help in identifying the critical group: stationary group: or
the critical cutset ijeJ of branches before one simulates a

particular fault. The theorem just states that if certain

162
constraints (5.17, 5.18) are satisfied, at some time t*
during the simulation and on some cutset, then a loss of
stability will occur. Note that the bus angle separation
with the stationary group and critical group can be 90° for
the case where the bus angles 6%(t), ijeJ are positive but
only A832° if 6%(6‘) are negative to ensure that the
acceleration experienced when

6ij(t) 3 180°—6‘i’j(t)
is not overcome before all bus pairs ijeJ satisfy

6ij(t) 3 360° + 6‘1’j(t)

The deceleration energy can be greater than the
acceleration energy only if 6%(t) is negative since when
623.01) is positive, the acceleration energy picked up as
angles pass through the band

180° - 6334);) g 6ij(t) $360 + 6‘i’j(t*)
will always exceed the deceleration energy as angles pass
through the band

360 +6‘i’j(t) g 61341:) <_54o - 6‘i’j(t)

If (595°, and 5%20 for all ijEJ, one can assure that
if a loss of transient stability occurs, then at some time
t* one has satisfied conditions (5.17, 5.18).

The constraints imposed in specifying the region of
instability constrain the entire state i(t),g(t) and §(t)
where the region of stability constrain only 5(t) and
specify both Q(t) and §(t). Thus the region of instability

better indicates conditions that cause instability than does

163

the conditions that 'specify the region of stability.
Furthermore, given that the region of stability is specified
by a sufficient condition for stability there is never any
way of assessing how conservative it may be but by exhaus-
tive simulation of many fault cases on a single system and
many such experiments on different systems. However, given
that for A532°, one has necessary and sufficient conditions
for loss of stability for both positive and negative values
on 6%.(t) ijeJ, then one can compare the region of

stability and the region of instability to determine how
conservative the estimate of the region of stability is.
Figure 18 shows the region of stability in branch angle
space and the region of instability for the case where J
contains two branches. Note that the region of stability

and instability touch at a point.

Region of
Instability

180-62 j (t*)
1 1 ‘1‘

 

Region of
Stability

 

 

 

‘\\.

Fig. 18 Regions of stability and instability
in branch angle space.

180-6. .(t*)
1232

164

One can say that if the trajectory enters the region of
instability at some t** a loss of stability will be experi-
enced and that if one experiences a loss of stability the
trajectory will enter the region of instability when A<45°
and 62j(t*)30. One can also state that if the trajectory is
within the region of stability for all time the trajectory
will remain stable. Nothing can be stated about stability
or loss of stability if the trajectory enters the area where
one or more of the angle differences satisfy

6ij(t) > 180 - 6fj(t)

but not all of the angle difference ijeJ satisfy (5.10-
5.12). Since one wishes to maintain a sufficient margin of
stability and one knows if all angles exceed (5.10-5.12) at
some time t* a loss of stability results. A prudent
constraint for assuring retention of stability is that no

branch angles 61 (t) even approach the boundary of the

3'
region of stability

6ij(t) 5 180 - 6§j(t) ijeJ
Thus the need to know whether the system is stable or
unstable in the areas where statements about retention or
loss of stability cannot be made may not be important in a
practical sense.

It should be noted that the results of this chapter
have

1) provided a theoretical basis for a region of

instability and its relationship to the region

of stability established from Popov's
stability criterion:

165

2) provided a theoretical indication of the
conservatism of the region of stability.
Although the‘ region of stability and insta-
bility do not cover the entire angle space,
they touch at a single point and this suggests
that the conservatimn of the region of
stability is modest. In a practical sense,
its conservatism may be considered small
compared to stability margins desired for
system security with respect to the particular
fault.

5.4 Stability Criteria Based on the Boundaries of the
Region of Instability and Stability

A performance measure

 

1:
v1 =12jea {o6ij(t)fij(6ij(t)) at /tf-to
tf
= fjea J[;(Ei(t)sj(t) Sin6ij(t)-E:E§ Sin dij)(5ij(t)-5?j)dt
tf_to

could be used to indicate if a loss of transient stability
occurs or does not occur as well as the time tc** when the
trajectory enters the region of instability. As long as the
trajectory remains in the region of stability the function
being integrated is positive and the integral will increase
monotonically with time. When the region of instability is
entered the function being integrated becomes negative and
remains negative.

The value of V1, which integrates components of the
Popov stability criterion on a cutset over the interval

t€[t0,t£], depends on the value of tC and the system

166
trajectory for this clearing time. The value of V1 for a

particular value of clearing time is denoted as V1(tc). If
tC increases for the case tc<tcc, the measure V1(tc) should
increase since the integrand will be larger in general over

telto,tf]. If tC exceeds the critical clearing time tcc'
the integrand for some portion of the interval [t°,tf) will
be negative. As tc increases for tc>tcc, the portion of the
interval [to'tf] when the integrand is negative becomes
larger and thus V1(tc) decreases. Computational results in

the next chapter confimm that the maximum value of Vl(tc)

should occur at the critical clearing time tc

c
.A second integral that represents potential energy on

the cutset

(t) fij (oij) dt

could also be used to indicate if a loss of transient
stability occurs or does not occur and the time to(tc) when
the trajectory enters the region of instability. If

Ei(t)=E? for i=1,2,....,N, then V2(t) is the outset energy

function
t
s s . . 52
= . . . .. - .. .. t
V2(t) I. 813E115):J (S1n 613(t) S1n 61]) d613( )
1jeJ
= - z ESESB [(Cos 6 (t) - Cos 651)-(6 (t) - 6§i) Sin 6??
i j ij ij ij ij 13 13

ijeJ

167

This cutset energy function increases monotonically

with increasing fiij

one clears the fault at tcgtcc the measure V2(t) increases

(t) in the region of stability. Thus if

until the trajectory reaches its maximum angular excursion
at tB(tc) and V2(t) reaches its maximum value at tB(tc). If

tC is increased but is less than tcc the value of the

maximum value of V(t)

v2 (t:B (to ))2---- mix {v2 (t)}
increases. The maximum value of V2(t) denoted V* would be
reached if at some point t if

6. . 180° - 6°. 5.20
13(1) j ( )

for all ijeJ. Experience has shown that the trajectory very

close approaches this point at tB(tcc) when tc=tcc. A

criterion for stability developed in [23] was that if
V2(tB(tc)) 5 V2*

where

2
V§=ijJBlJlJ{Cos6§+Cos6][6318111633.}

If Ei(t)=B? and in addition there are no load buses so that
N=Ng, then the model reduces to the classical transient
stability model. The critical cutset would correspond to
the branches connecting the internal generator bus of a
critical generator to the internal generator buses of other

generators. In this case, the criterion for stability is

V2(tB(tc)) 5 V2(tB(t*c))

168

where t*c is selected based on

max

V2(tB(t*C)) = to V2(tB(tc))

The results in [21] show that the critical clearing time can
be determined with no detected error, i.e. t*cgtcc' This
method is called the PEBS method in [21] and V2(t) is the
potential energy component of the individual machine energy
function.

It is clear that the boundary of the region of stabi-
lity and the region of instability touch at a point in angle
space satisfying (5.28) where the performance measure V2=V2*
and where the slope of S¥2(t)=0. Thus stability criteria can
be established based on observing measure V2(t) as a
function of time for a cutset J rather than checking the
angles 5ij(t) for all ij EJ and all tto determine loss of
stability» It is also clear tht the performance criterion
V2(t) is related to the cutset energy function stability
criterion if Ei(t)=B? for all i=1,....,N and is related to
the PBBS stability criterion for the individual machine
energy function if 82(t)=8? for i=l,2,....,N and N=Ng. Thus
the stability criteria developed based on simpler models is
justified based on the region of stability and region of

instability derived in this chapter.

CHAPTER 6

VERIFICATION OF
STABIILTY CRITERIA VIA SIMULATION

The objective of this chapter is to verify the accuracy
and validity of the stability criteria proposed in the
previous chapter.

The first stability criterion to be tested is based on
the Popov stability criterion. The theoretical results of
the previous chapter indicate

1) that the system will be stable if the Popov

stability criterion is satisfied for all
branches over the entire simulation interval

2) that the system will lose synchronism if one

can find a cutset of branches for which the
Popov stability criterion is violated at some
time t*.

The second stability criterion is based on measure V1
that will continue to increase if the state is in the region
of stability but will decrease sharply after the trajectory
enters the region of instability. The stability criterion
determines the value of the maximum value of V1 for the
interval [to,tf] for each clearing time for which the

transient stability simulation criterion is run. The value

of V which integrates components of the Popov stability

1'
criterion on a cutset over the interval telto,tf],depends<n1

the value of tc and the system trajectory for this clearing

169

 

170

time. The value of V1 for a particular value of clearing
time is denoted as V1(tc). If to increases for the case

tc<tcc, the measure V1(tc) should increase since the

integrand will be larger in general over te[t0,tf]. If tc

exceeds the critical clearing time, too' the integrand for
some portion of the interval [to,tf] will be negative. As
tC increases for tc>tcc, the portion of the interval [to'tf]
when the integrand is negative becomes larger and thus
V1(tc) decreases. Computational results in this chapter
confirm that the maximum value of V1(tc) should occur at the
critical clearing time tcc.

Although these stability criteria were developed for a
stability model that includes a topological network model, a
single axis generator model, a constant power real load
model, and a constant current reactive load model, the

stability criteria will be tested and verified for

a) a single axis generator model, constant power
load model, constant current reactive model;

b) a constant voltage behind transient reactance
generator model, a constant power real load
model, and a constant current reactive load
model;

c) a constant voltage behind transient reactance
generator model, a constant power real load
model, and a constant impedance reactive load
model.

The simulation of transient stability has been

performed for two fault cases on the 39 bus New England
system model. The two fault cases are a fault on line 26-27

and a fault on line l4-33 where in both cases the fault is

cleared by removing the line.

171

The first stability criterion is tested by plotting

°ij(t) fij(oij(t)) (6.1)

as a function of time for branches that are geographically
close to the fault location. It will be seen that this
function will remain positive for branches far from the
fault, will approach zero for branches near the fault when
the clearing time is close to the critical clearing time,
and will go negative for branches close to the fault when

t >t

c cc' and will reach very large negative values and

oscillate for branches in the critical cutset when tc>tcc.
Although the critical cutset is unknown before the simula-
tion runs are performed for a particular fault to determine
the critical clearing time, it is easily determined from the
simulation results.

The second stability criterion

 

tr
V(‘t )= E B..] [Eo(t)E-(t) Sin 5--(t)
C ijeJlJ to 1 J 13
S S . S S

is evaluated cn14a transient stability simulation of a
particular fault cleared at tc and the critical outset of
branches ijeJ determined based on the results of the tests
of the first stability criterion. The function V(tc) is

plotted as a function of tC and the value t*C for which

172

V(tc) has a maximum value should be too‘ Retention of
stability would be predicted to occur for tc<t*c.

Four different fault cases have been presented in the
results that follow:

Case 1) Fault on line 26-27 (flux linkage of generators
considered and constant current reactive load model): Using
the transient stability program the critical clearing time
for this case is determined to be .37 seconds.

By considering the variations of angles of generator
for the unstable case it was observed that the first
generator that loses stability is generator No. 9. This
conclusion could be anticipated before considering the
transient stability programs results because generator No. 9
is the closest generator to the location that fault
occurred.

The variations of Oij(t) f (‘fij(t) for the stable and
unstable case are plotted as functions of time for lines
28-29 and 9-29 and 26-28 and 26-25 and generator 9's tran-
sient reactance in Figure 19 for the first three seconds.

In the stable case the values of Oij f (oij(t)) are
positive and are often close to zero. This justifies the
claim that %j fij (Oij) should be greater than or equal to
zero in the region of stability and should approach and
remain close to zero when tcatcc.

The values of oij(t)f(Oij(t)) are also plotted for the
unstable case when tc=tcc* in Figure 19. The values

oij(t) fij (Oij(t)) are positive for all branches for the

173
first 1.2 seconds at which time the values on every branch‘
become negative suggesting the loss of stability is
imminent. The value of 0ij(t) f (Oij(t)) on the transient
reactance on generator drops to a large negative number,
remains constant for a short period and then begins
oscillating (the values below a certain value were not
plotted and when the curve returns to the range of the
coordinates, the curve is again plotted.)

The large negative and positive oscillations are not
observed on the plots for the other branches and thus the
transient reactance of generator 9 is the critical cutset.
It should be pointed out that after the first large drop in
generator 9's 0.

13'
transformer has a large increase. This appears to occur

f (oij), line 9-29 which represents a

because the transformer is trying to'absorb the excessive
energy of the generator and prevent it from loss of
synchronism.

The measure of V(tc) given by (6.1) was computed for
various clearing times when tf=3.0 and the transient
reactance of generator 9 is the sole branch belonging to the
critical cutset. The value of V(tc) increases as tc
increases. The maximum value of V(tc) occurred at the
critical clearing time determined from simulation results.
The function V(tc) decreases rapidly for increase in tc
above tcc’ These results confirm the accuracy and validity
of these stability criteria as well as the theoretical

foundations upon which the stability criteria were based.

174

Case 2) Fault on line 14-33 (single axis generator
model and constant current reactive load model): Using the
transient stability program the critical clearing time for
this case is determined to be .24 seconds.

In this case using the same procedure of case 1, the
pair of generators 2 and 3 were determined to be the first
generators that lose stability.

The variations of 0.. f-- (Oij) for the lines that are

13 1]
close to critical generators and 0;. f.. (Oij) for transient

1] 1]

reactance of these generators are plotted in Figures (21.a),
(21.b), (21.c), (21.d) and (21.a). With the same analogy as
in case 1, the transient reactances of the group of
generators l and 2 are identified as the critical cutset.
The large negative values and the large oscillations
observed on the branches in the critical cutset are not
observed (”1 the other branches. ‘The integral measure
defined in (6.1) is plotted in Figure 22 and again correctly
predicts the critical clearing time. The results of this
case are similar to the results of case 1 and both cases
support the theory presented in Chapter 5.

The values of Oij ij ij) remain positive and
close to zero for the critically stable case transient

(t) f (a
stability. The values of oij(t) f (oij) on most of the
branches go negative for a short period when the trajectory
enters the region of instability followed by a period of
large positive variation. The cij(t) f (Oij) values for the

two branches in the critical cutset go negative, attain a

175
constant value and then start oscillating. The magnitude of
the oscillation are large values.

Case 3) Fault on line 26-27 (constant current reactive
load model and constant voltage behind transient reactance
generator model): Using the transient stability program the
critical clearing time for this case is determined to be
.36 seconds. The results of this case are very similar to
case 1 because the models are identical except that the
single axis generator model is replaced by a constant
voltage behind transient reaCtance generator model. The most
significant difference between the cases is that generator 9
does not oscillate near the boundary of the region of
stability in this case but quickly enters the region of

instability. The values of o. f

13'
generator 9's transient reactance are plotted in Figure 23.

ij (oij) for lines and
The integral measure is plotted in Figure (24). These
results show that despite the fact that the model of the
power system in this case is a little different from the
model used in Chapter 5, the theory introduced in that
chapter could also be used in this case.

Case 4) Fault on line 26-27 (constant impedance
reactive load model and constant voltage behind transient
reactance generator model): By using the transient stability
simulation program the critical clearing time in this case
is determined to be .31 seconds.

Despite the fact that the model used in deriving the

theory in Chapter 5 is different from this case the

176

stability criteria based cum f ( ) is successful in

0.. 0..
13 1]
determining the critical clearing time and the cutset and
the integral measure can accurately determine the critical
clearing time.

The critical cutset, similar to previous cases, is
generator No. 9's transient reactance and the value of
integral measure for this cutset is positive and increasing

as tc gets close to tcc and drops to negative when tc>tcc.

The results of this case are shown in Figures 25 and 26.

 

177

Fig. (19) Comparison of o..f(oi.) for stable and unstable
case when the Eahlt obcurs on line 26-27. .The
fault has been cleared and line 26-27 has been
removed. (Lyapunov energy function)

' AA A /\ A___/\ [\JV
' v v v Y7 fi

- ' 'W
40 I a I

f l_ . F\j\fkmAAAl
I \l i V 3

f

LINE 28-20

 

r
I

{3:0
I
32:
x a
34“ f 3 ~~ .
l a Y —T

001 .
o f - fi/\_ A A ‘ AA;

7!”! [I SECONDS

LINE 23-29

 

:3 . NR -/\Ar/\V/\‘VAJ\./\t/‘

 

LINE 20-29
é
”)
<

i v , _B-fl/\ARAA
-fl ' %9%A/\H/UVUY

0i“. ”0. 9

178

Fig. (20) Values of V (t t ) for the single axis generator
model and he constant current reactive load
model for the fault on line 26-27.

12.

..-----

 

o.i oi of: 0.37 "c
-4. TIME IN SECONDS

 

~201

179

Fig. (21.a) Comparison of o..f(oi.) for stable and unstable
case when the fghlt otcurs on line 14-33. The
fault has been cleared and line 14-33 has been
removed. (Lyapunov energy function)

 

‘ . J 3A V/\ fi/\’\ /fi
2 V ,

LINE 12-13
<3

-1

 

 

«.7
2

LINE "-13
O C

O
_)

flflE IN SECONDS

Fig.

-OO

LINE 30-1!

—
b N

LINE 30-13

I-‘N
D.

0

LINE 33-34

S
.a
N
o

-2001

I
N
a N O
a e O O
A L L L A 4L_4 l A

I
O
A

(21. b) Comparison of <1.

180

f(0

for different lines on

stable and unstable else when the fault occurs
The fault is cleared and line
(Lyapunov energy function)

on line 14-33.
14-33 is removed.

 

 

-1

f

"1

 

 

 

VINE IN SECONDS

" I

 

“1
i)

7
F

181

Fig. (21.c) Comparison of Oi.f(0..) for different lines on
stable and unstable ddse when the fault occurs
on line 14-33. The fault is cleared and line
14-33 is removed. (Lyapunov energy function)

azi _

IS

0 . _,.__—-—-1———k I AVA A
is

4

LINE 33-32

-IS
78

FL
P

I I
9 N
cu m

1 v i v i VA V w}
-IZS

300.

LINE 38-36

.:>
f>
2

400:

2:1.
34 A /\
. . . . /\{
1 V 1

-24.

 

LINE 35-11
0
(a

2

3

I“

3

ca

In

r:

3

O

a»
I
K

182

Fig. (21.d) Comparison of Oi.f(01.) for different lines on
stable and unstable c33e when the fault occurs
on line 14-33. 14-33 is removed after the fault

is cleared. (Lyapunov energy function)

280

.— (3
O
0

LINE 35-3:

 

LINE 31-36
N a
O G
22

LINE 31-3.
I
O
O
1
(

 

LINE 37-3.

“£001 I INE IN SECONDS /

183

Fig. (21. e) Comparison Of 0 ) for stable and unstable
case when the fELmlt io’ccurs on line 14-33. The
fault has been cleared and line 14- 33 has been
removed. (Lyapunov energy function)

H Y I - 2 Y ‘7
‘200
4
-400.
1
‘500.

' v a

'75: YIHE IN secouos

 

GEN. N0 2

 

 

184

Fig. (22) Values of V1(t ,t ) for the single axis gene-
rator model and the constant current reactive
load model for the fault on line 14-33.

  

 

 

I
Nd_‘-—o - -- o. q - - _ c. - -- ...

0.1 0.2 ‘.
TIME IN SECONDS I
)

t
-4. \

 

185

Fig. (23) Comparison of o..f(o..) for stable and unstable
case where the faultldccurs on line 26-27. The
fault is cleared and line 26-27 is removed.
(Constant current load model)

 

 

 

 

 

 

 

 

 

 

HS,
gnzfi
$‘B.
“.26
z ~V/Y\AT/\* kM/\ /\
3‘25; I ' r - ' *1
:5}
m
7 m‘ /«\\a;~[>\=ZB\JCLJfi\
3 4——-r‘—~5J . . ,/NL1 f\\
“-15. I i \J j
5 1
_J-4S
03‘51
N J
S m. /\\\\ /\\ Z 5 [§ [5 /\
N (\1 ‘WV /\
“-15. V ' . I V' I
5 I
4‘45
7D.
:3 so}
ré 30.
.u ID) \
z . \/ v ‘\V. 1 ‘~
3 -HJ. r - , e 4—7
TIHF IN SFrnnD:

03

q: 6

Q

N 2

Eg ,/’I—‘=z::—77Z::ExzjCI7Z€§¥Z7Qfo\éZCBbz(:kl[fl\_\/[\\

3-2) ' K] I ' 5 -.

TlflE IN SECONDS

500 . - __H- h-_-_.
c”250] A A A
.o . .n ,. /\/\/\/~.
”mi \'/n)nH\/’\/\/\/‘ *3
Z
8'500 NSCDDS \

186

Fig. (24) 'Values of Vl(t ,t ) for the constant voltage
behind tran31en€ rgactance generator model and
the constant current reactive load model.

 

1
4. .
' Ii 0.2 1'3 '
_4 TIME IN SECONDS
-12.
-2O.l

 

187

Fig. (25) Comparison of Oi.f(oi.) for different lines on
stable and unstable ckse when fault occurs on
line 26-27. The fault is cleared and line 26-27
is removed. (Constant impedance load model)

60.
304
+
0

4 V v fi— 5 V
VINE IN SECONDS \V/

...J

LINE 20-29

 

LINE 09-29
a:

I
N
-

 

f 2 r w j
35‘ VINE IN SECONDS

LINE 26-20

25.

15:

5‘ </”—_‘§§E§__ ,//A\f:;

-5 ' ' r I r Y",14_3

\
A
v A . I

. ' v . ' Y I

-60.
400) T
[x , _F

I \/\ I
-60. TIN! In sscnnos ‘/

-IDO) ‘J

° ' i . .lifi
-I:O ' \ ’
-zco

LINE 26-23

 

LINE 26-25

 

GENN09
.
O

188

Fig. (26) Values of V1(t ,tc) for the constant voltage
behind transi nt reactance and constant
impedance load model.

   

 

4'1
. e E :
" .I .2 .3 .36 '
_4 TIME IN SECONDS
-12.

 

-201

CHAPTER 7

OVERVIEW AND FUTURE RESEARCH

7.1 Overview

 

An overview of the research performed in this thesis
and its contributions are now provided. A review of the
literature on the development of Lyapunov functions and
energy functions is given in Chapter 2, along with discus-
sion of the equal area and P888 methods for determining
retention or loss of stability.

In Chapter 3, a topological energy function has been
derived from the first principles that retain the network
without aggregation back to generators internal buses
and includes a description of real and reactive power load
models as a function of voltage. This energy function has
been tested on a 39 bus system and the results are presented
at the end of Chapter 3. This is the first attempt at
deriving a topological energy function using the energy
integral method. This topological energy function also
permits both a general real and reactive load model descrip-
tion for the first time.

In Chapter 4, a Lyapunov energy function is derived for
a model that

l) retains the network and does not aggregate the
network back to internal buses;

189

190

2) models the real power load as constant power

and the reactive power load as constant
current;

3) models the flux linkage decay in the synchro-

nous machine as well as the electromechanical
component of the machine model.

This Lyapunov function has been derived using the Popov
criterion, the Moore-Anderson theorem [10] and the construc-
tion method of Willems [3]. The last section of Chapter 4
contains the computational test of this Lyapunov energy
function. Although a topological Lyapunov has been derived
previously, it did include flux decay effects in the
generator model and completely ignores both the reactive
load power energy, and the reactive energy produced by the
generators. The Lyapunov energy function can be shown to be
identical to the energy functions derived in Chapter 3 when'
the real load power model is constant power, the reactive
load model is constant current and the generator model is
the constant voltage behind transient reactance generator
model.

In Chapter 5, a discussion of the region of stability
based on the Popov stability criterion and a definition of
loss of transient stability has been presented. Three
different theorems stating necessary and sufficient condi-
tions for the loss of transient stability have been
established. The conditions stated in those theorems
describe an estimate of the region of instability. Although

a precise boundary between the region of stability and the

region of instability is not established, a relationship

191

between the regions of stability and instability has been
established. At the end of Chapter 5, a cutset integral
criterion based on the theoretical description of the
regions of stability and instability has been introduced and
justified. The cutset energy function criterion [7] for the
topological function and the PEBS method [21] for the
individual machine energy function were also justified based
on this theory.

A description of the region of stability and insta-
bility determined has a common boundary point which
establishes for the first time the level of conservatism of
the region of stability developed based on the Popov stabi-
lity criterion. The necessary condition for instability
provides a condition that must have occurred at some
instance before a loss of stability oCcurs. These necessary
and sufficient conditions provide testable conditions on the
entire state of the system that can be clearly identified
with the loss of synchronism inherent in the loss of
transient stability.

Two stability criteria based on the Popov stability
criterion and based on the outset integral measure are
tested in Chapter 6. The criterion that established the
critical outset and the critical clearing time based on
a) the satisfaction of the Popov stability criterion for all
t_>_tc or b) satisfaction of the necessary and sufficient
condition for instability at some t*>tC appears very

accurate and reliable based on the four fault cases studied.

192
The stability criterion based on the outset integral measure
V1(tc) is shown to be a very accurate method of determining
the critical clearing time. Since the integrand will be
positive and generally increase for time te[to,tf] as tc
increases for tc<tcc, and since the portion of the interval
[to,tf] where the integrand is negative increases with tc
for tc>tcc, the value of clearing time tc at which V1(tc) is
maximum clearly and accurately determines tcc' The computa-
tional results on all four fault cases verify the accuracy
and validity of this stability criterion. These results for
the first time provide theoretical justification of
stability criteria and establish why these criteria have no
apparent error in determining retention or loss of transient

stability.

7.2 Future Research

The Lyapunov energy function introduced in Chapter 4
has been derived for one axis generator model for the case
that reactive loads are a first order function of voltage
(constant current). An important extension would be the
derivation of a Lyapunov energy function with a) a more
complex generator model including both the exciter and power
system stabilizer and b) a general real and reactive load
model.

Another useful and important investigation would be to
develop a precise boundary between the region of stability

and the region of instability introduced in Chapter 5.

193

The sensitivity of the defined regions of instability
or stability could be found for operating conditions such
as:

Network configuration
Energy transfers

Generating dispatch

Unit commitment

Voltage profile and support
HVDC.

mmbUNH
vvvvvv

The sensitivity of the region of stability could be used to
» derive security constraints that assure sufficient stability
margins for retention of stability for these various changes
in operating conditions.

Another investigation could be the development of a
fast computationally efficient program that determines
whether the system trajectory remains within or crossed the
boundary of the stability region. An initial task of such a
fast program would be the identification of the critical
outset and the critical and stationary groups defined in
Chapter 5. This fast computationally efficient method could
be developed based on either a Taylor series or RMS measure

methods introduced by other researchers.

10.

LIST OF REFERENCES

Magnusson, P.C., "Transient energy method of calculating
stability,” AIEE Trans., Vol. 66, 1947, pp. 747-755.

Aylett, P.D., "The energy integral criterion of tran-
sient stability limits of power systems,” Proceedings
of IEE (London), July 1958, pp. 527-536.

Willems, J.L., "Optimum Lyapunov Function and stability
regions for multi-machine power systems,” Proceedings
of IEE (London), Vol. 117, No. 3, March 1970, pp.
573-578.

Guduru, V., “A general Lyapunov function for multi-
machine power system with transfer conductances,"
Int. J. Control, Vbl. 21, No. 2, 1975, pp. 333—343.

Kakimoto, N., Y. Ohsawa, and M. Hayashi, ”Transient
stability analysis of electric power system via Lure
type Lyapunov function," Parts I and II, Trans. IEE
of Japan, v01. 98, No. 5/6, May/June 1978.

Kakimoto, N., Y. Ohsawa, and M. Hayashi, ”Transient
stability analysis of multimachine power systems with
field flux decays via Lyapunov's direct method," IEEE
trans. on power apparatus and systems, V01. PAS-99,
No. 5, September/October 1980.

Hill, D.J. and A.R. Bergen, "Stability analysis of
multimachine power networks with linear frequency
dependent loads," IEEE trans. on circuits and
systems, Vol. CAS-29, No. 12, December 1982.

Athay, Thomas M., and'David J. Sun, ”An Improved Energy
Function for Transient Stability Analysis," IEE
International Symp. on Circuits and Systems, 1981.

Musavi, M.T. and N. Narasimlanurthis, ”A general energy
function for transient stability analysis of power
systems," IEEE trans. on CAS, No. 7, July 1984.

Moore, J.B. and 3.0.0. Anderson, 1968, J. Franklin Int.
488-492. "

194

11.

12.

13.

14.

15.

16.

17.

18.

19.

20.

21.

22.

195

Willems, J.L., ”Improved Lyapunov function for transient
power system stability,” Proceedings of IEE, Vol.
115, No. 9, September 1968.

Desoer, C.A. and M.Y. Wu, "Stability of a nonlinear
time-imminent feedback system under almost constant
input," Automatica, Vbl. 5, pp. 231-233, 1969.

Pai, M.A. and P.G. Murty, “New Lyapunov functions for
power systems based on minimal realization, Int. J.
Control 1974, Vol. 19, No. 2, pp. 401-415.

Anderson, B.D.0., ”A system theory criterion for posi-
tive real matrices," SIAM J. Control. Vol. 5, No. 2,
pp. 171-182, 1967.

Kakimoto, N., Y. Ohsawa, and M. Hayashi, ”Transient
stability analysis of large-scale power system by
Lyapunov direct method,” IEEE transaction on power
apparatus and systems, July 1983.

Sastry, H.S.Y., ”Application of topological energy
functions for the direct stability evaluation of
power systems, a Ph.D. thesis submitted to .the
DeparUment of Electrical Engineering, University of
California at Berkley, March 1984.

Rastgoufard, P., ”Local energy function methods for
power system transient stability,“ a Ph.D. thesis
submitted to the Department of EE/SS, Michigan State
University, December 1982.

Michel, Fouad, Vitall, "Power system transient stability
using individual machine energy function," IEEE
trans. on circuits and systems, Vol. 5, May 1983.

Vitall, Michel, ”Stability and security assessment of a
class systems caverned by Legerange's with applica-
tion multi-cavern system,“ IEEE trans. on CAS,
November 1983.

Michel, Miller, “Qualitative .Analysis of Large Scale
Dynamical System," Academy Press, 1977.

Yazdankhah, A.S., ”Power systems security assessment for
faults using direct methods," a Ph.D. thesis
submitted to the Department of EE/SS, Michigan State
University, July 1984.

Chandroshekar, K.S., D.J. Hill, "Dynamic security
dispatch: basic formulation,” IEEE trans. on PAS,
July 1983.

196
GENERAL REFERENCES

23. Kimbark, E.W., "Power System Stability,” John Wiley and
Sons, 1948.

24. Crary, S.B., ”Power System Stability,” John Wiley and
Sons, 1945.

25. Anderson, P.M. and A.A. Fouad, "Power System Control and
Stability," Iowa State University Press, Ames, Iowa,
1977.