Citation
An investigation of the power system stability domains for a two-finite-machine system

Material Information

Title:
An investigation of the power system stability domains for a two-finite-machine system
Creator:
Dghayes, Hussein Abdullatif
Publication Date:
Language:
English
Physical Description:
94 leaves : illustrations, charts ; 28 cm

Subjects

Subjects / Keywords:
Electric power system stability ( lcsh )
Electric power system stability ( fast )
Genre:
bibliography ( marcgt )
theses ( marcgt )
non-fiction ( marcgt )

Notes

Bibliography:
Includes bibliographical references (leaves 66-67).
General Note:
Submitted in partial fulfillment of the requirements for the degree, Master of Science, Department of Electrical Engineering and Computer Science.
Statement of Responsibility:
by Hussein Abdullatif Dghayes.

Record Information

Source Institution:
University of Florida
Rights Management:
All applicable rights reserved by the source institution and holding location.
Resource Identifier:
19719705 ( OCLC )
ocm19719705

Downloads

This item has the following downloads:


Full Text
AN INVESTIGATION OF THE POWER SYSTEM STABILITY

DOMAINS FOR A TWO-FINITE-MACHINE SYSTEM
by
Hussein Abdullatif Dghayes
B.S. University of Colorado, 1985
A thesis submitted to the
Faculty of the Graduate School of the
University of Colorado in partial fulfillment
of the requirements for the degree of
Master of Science .
Department of Electrical Engineering and Computer Science
1988


This thesis for the Master of Science degree by
Hussein Abdullatif Dghayes
has been approved for the
Department of
Electrical Engineering and Computer Science
by
Edward T. Wall
William R. Roemish
Date / ? *> 5 6~"


Dghayes Abdullatif, Hussein.(M.S., Electrical Engineering)
An Investigation of the Power System Stability
Domains for a. Two-Finite-Machine System
Thesis directed by professor Edward T. Wall
Abstract- Analysis of the stability of power system
following a transient disturbance involves the study of a large
set of non-linear differential equations. This thesis presents
a systematic procedure for investigation of stability domains of
a two-finite machine, interconnected by passive loads, power
system model in which the classical representation is used for
each machine and resistance is neglected. The direct method of
Liapunov is used to determine the stability boundary of the power
system. The method is first compared with the equal-area
criterion and the phase plane technique, with all three approaches
giving identical results for the simplified model With constant
damping introduced into the system model, a larger, more accurate,
region of stability is obtained by use of Mansours modification
of the Energy-Metric-Algorithm.


IV
ACKNOWLEDGEMENTS
This thesis is submitted to the University of Colorado as
the final requirement for the degree of Master of Science in
Electrical Engineering studies at the Denver campus.
This study would not have been possible without the
cooperation and support of many people to whom I will always be
grateful.
I especially wish to express my sincere appreciation to
Dr. Edward T. Wall, for his cheerful support and encouragement
during the preparation of this thesis. Also, the enthusiasm of
Dr. William R. Roemish for his suggestions and encouragement
throughout this study.
I would further like to thank the faculty of the
University of Colorado at Denver for their efforts in providing
an excellent education.
Finally, I would like to thank my family, for their moral
support and inspiration.


CONTENTS
CHAPTER
I. INTRODUCTION ........................................ 1
II. STABILITY ........................................... 4
2.1 Introduction ...................................... 4
2.2 Steady-State Stability ............................ 4
2.3 Dynamic Stability ................................. 5
2.4 Transient Stability .............................. 6
2.5 Methods of Simulation ............................. 7
III. CLASSICAL MODEL
3.1 Introduction ...................................... 9
3.2 Classical Stability Study Assumption .............. 9
3.3 Power System Representation ...................... 10
3.4 System Equation of Motion ........................ 10
3.5 Pre-Fault Power System Model ..................... 14
3.6 During Fault Power System Model .................. 18
3.7 Post-Fault Power System Model .................... 20
IV. NUMERICAL SOLUTION OF THE SWING EQUATION
4.1 Introduction ..................................... 23
4.2 Point-by-Point Solution of the Swing
Equation .......................................... 23
4.3 Numerical Methods for Solution of Differential
Equation
27


VI
V. STABILITY DOMAIN
5.1 Introduction ...................................... 30
5.2 Power-Angle Characteristics ....................... 30
VI. APPLICATION OF THE PHASE PLANE TRAJECTORIES
AND THE STABILITY BOUNDARY
6.1 Introduction ...................................... 35
6.2 Phase Plane Trajectory ............................ 36
6.3 Stability Boundary ................................ 42
VII. ASYMPTOTIC STABILITY
7.1 Introduction ...................................... 46
7.2 Liapunov Stability ................................ 46
7.3 Energy Metric Algorithm(EMA) ...................... 47
7.4 A Suitable Liapunov Function and
an Estimation of X ................................ 53
VIII. MANSOURS MODIFICATION OF THE EMA
8.1 Introduction ...................................... 57
8.2 Power System Representation ....................... 57
8.3 System Equation of Motion ......................... 58
8.4 The Unmodified EMA ................................ 59
8.5 Mansours Modified EMA ............................ 60
IX. CONCLUSIONS ............................................ 64


vii
BIBLIOGRAPHY ............................................... 66
APPENDIX
A. LISTING OF COMPUTER PROGRAM FOR THE SOLUTION
OF THE SWING EQUATION...................................... 68
B. LISTING OF COMPUTED PHASE-PLANE TRAJECTORY
COORDINATES ............................................... 76
C. THE ENERGY METRIC ALGORITHM ............................... 89
D. NOTATION
91


Vlll
TABLES
Table
3.1 Generator Data (in p.u to Generator MVA Base) ........ 11
3.2 Transmission Line Data On 100 MVA Base
(Resistance Neglected) .............................. 11
3.3 Generator Data (in p.u to 100 MVA Base) .............. 12
3.4 Load-Flow Data ........................................ 12
3.5 Load-Flow Data (on 100 MVA Base) ..................... 12


FIGURES
Figure
3.1 Single-Line Diagram Of Two-Machine Power System .. 11
3.2 Pre-Fault System Reactance Diagram .................. 15
3.3 Pre-Fault Admittance Diagram ........................ 15
3.4 Pre-Fault Power-Angle Curve ......................... 17
3.5 During Fault System Reactance Diagram .............. 19
3.6 During Fault Admittance Diagram ..................... 19
3.7 Post-Fault System Reactance Diagram ................. 21
3.8 Post-Fault Admittance Diagram ....................... 21
4.1 Actual and Assumed Values of Pa, u, and S
as Function of Time ................................. 26
4.2 Plot of 61 and 62 vs. Time ........................ 28
4.3 Plot of Delta Differences (61-62) vs. Time .......... 29
5.1 Classical Power-Angle Curve ......................... 31
6.1 Phase Plane Plot Of X2 Vs. Xl ....................... 37
6.2 Power-Angle Curves And Phase Plane
Trajectories ........................................ 40
6.3 Plot of Stable and Unstable Phase Plane
Trajectories ........................................ 45
7.1 Determination of the Stability Boundary
V-Function Value, K, form the Postfault System
Power-Angle Curve ................................... 51


X
7.2 Relationship Between % and the Mechanical
. Shaft Power Pml2 ................................... 52
7.3 Effect of the Equivalent Mechanical Power (Pml2)
on the Phase Plane Postfault Stability Boundary .. 57


CHAPTER I
INTRODUCTION
The dependence of modern society on electrical energy
requires that major power failures be avoided. The loss of
synchronism between generators and concentrated loads, caused by
power surges or system faults, presents a potential cause of power
failures. In order to evaluate the hazards and to take steps to
prevent loss of synchronism, accurate methods of analysis of
on-line system stability must be developed and put into use.
Current practice usually involves analytical studies which result
in system design or operating procedure modifications. Even with
the use of modern digital computer modeling, this approach does
not fully protect the modern large scale system during emergency
fault conditions.
This thesis considers the simulation of a proven system,
which is modelled by a two-machine equivalent system, during
possible fault and erratic operating conditions. A goal is to
seek, by analytical investigation, the range of operating limita-
tions for continual safe and reliable operation of such a power
system. The basic approach used involves a computer stability
simulation of the composite generator groups and equivalent loads.


2
A comprehensive presentation of the class of problems involved and
the analytical models currently in use, together with selected
models involving a stability analytical tool known as the second
method of Liapunov, is given. This method together with the
complex problem of developing Liapunov functions is developed and
compared by example with the classical semi-graphical methods in
current use, for relative accuracy and convenience in application.
Transient stability involves the study of the systems
capacity to recover, after a disturbance, to synchronous steady-
state operation. Such a study should consider the action of
voltage and speed regulators which would normally require the use
of Parks transformation[14] and also a consideration of vari-
ations in mechanical power. In addition, damping parameters which
complicate the analysis, should be included. However in this
thesis a first approach to the problem which gives conservative
results is used. This is an extension of the standard procedure
used in transient stability studies known as the step-by-step
method.
The problem may be formulated as follows: given a system
initially in a steady operation and assume a disturbance at time
tg. Then the question: is there a stable equilibrium position for
the system after the disturbance is cleared, and if so, what is
the critical clearing time, that is, the maximum time that the
perdisturbances may remain, before the system loses its capability


3
to return to steady-state. The advantage of the use of Liapunovs
second method as developed in this thesis, is that the extensive
computing time required in the last phase, is replaced by a
stability criterion (which is dependent on the construction of
a Liapunov function) which allows the determination, in the state
space, of the stability domain which surrounds the stable asymp-
totic solution of the post-fault system. This method is discussed
in this study. The success of the approach depends upon the
development of the Liapunov function, which is a main concern of
this thesis. It is essential to develop the best possible
Liapunov function for the model under study. This function is
developed to produce the largest possible asymptotic stability
domain. The method selected for the generation of the Liapunov
function for the two-machine equivalent system is the Energy
Metric Algorithm as extended by Mansour. This algorithm permits
the Liapunov function to be generated directly from the system
equations. This approach is very valuable in power system
stability studies.


CHAPTER II
STABILITY
2.1 Introduction
Power system stability concerns the operating condition
in which the various synchronous machines in a power system remain
in synchronism, or in equilibrium with one another. A stable
electric power system returns to normal operation after a dis-
turbance or change in the operating condition, the most severe
type of disturbance to which a power system is subjected is a
three phase short circuit.
Stability studies will be divided into three different
categories depending on the extent of the disturbance on the
system :
2.2 Steady-State Stability
Steady-state stability is the ability of the synchronous
machines of a system to remain in synchronism with each other
subsequent to the occurrence of small disturbances. These small
disturbances include very small gradual load changes which can
occur on a continuous basis during the systems normal operations.
This type of stability is realized by the use of manual control
devices, since the majority of the excitation systems in use have


5
sufficient deadband and long enough time delays so that the
voltage regulator action ( A E^= 0 ) can be neglected for small
disturbances. The steady-state stability limit refers to the
maximum flow of power possible through a particular operating
point without the loss of synchronism and when the power is
increased very gradually.
2.3 Dynamic Stability
Dynamic stability concerns the ability of the generators
of a large interconnected system to remain in synchronism over
long periods of time following small changes or normal random
disturbances in the dynamic system. The stability limit is
greater than the steady state limit. This type of stability can
be controlled by the use of automatic control devices and voltage
regulators. Since the excitation systems are used without
deadband, the voltage regulator action (A 0) cannot be
neglected for either small or severe disturbances.
In dynamic stability analysis the system equations may be
linearized about the operating point. Stability of the linear
time-invariant systems will be determined by use of the following
methods: Routh- Hurwitz criterion, root-locus, and Nyquist
criterion.
Using the state space form :
X = AX + Bu , (2.1)


6
the stability characteristics of a system may be determined by
examining the eigenvalues of the A matrix. In the state space
equation above, X is an n vector denoting the states of the system
and A is a coefficient matrix. The system inputs are represented
by the rth vector u, and these inputs will be related to the
differential equations by an n x r matrix B. This description has
the advantage that A may be time-varying and that u may be used to
represent several inputs including system load responses.
2.4 Transient Stability
Transient stability is the ability of the machines to
remain in synchronism subsequent to the occurrence of large dis-
turbances in the system. The disturbance may be a sudden increase
in load, or switching operations which result in the loss of a
major system component, such as a large generator or a long trans-
mission line. However, the most severe type of disturbance is the
three phase fault. The limits of stability in these cases are
defined to be the maximum flow of power possible through a given
operating point without the loss of synchronism at the time of the
disturbance.
The system equations for a transient stability study are
nonlinear in nature and the study of the first swing cycle is most
important. The system is described by a set of coupled non-linear
differential equations of the form :
X = f( x u, t) (2.2)


7
Where f is an n vector of nonlinear functions. The determination
of the dynamic behavior of the system described by Eq. (2.2), is
more difficult than that of a linearized system such as for Eq.
(2.1). Digital computers are extensively used to analyze the
stability off line for reliability studies but difficult to use
on line for emergency decisions during faults.
2.5 Methods of Simulation
The first step in a stability study is to make a mathemat-
ical model of the system during the fault. The elements included
in the model are those affecting the dynamics of the machines
(acceleration or deceleration of the machine rotors). Generally,
the elements of the power system that influence the electrical and
mechanical torques of the machines should be included in the
model. These elements are :
1. The network before, during, and after disturbance.
2. The parameters of the synchronous machines.
3. The loads and their characteristics.
4. The excitation systems.
5. The mechanical turbines and speed governors.
6. System components such as transformers, and capacitors.
The complexity of the model depends upon the type of
transient and system being investigated. Using the techniques of
modern control theory, the system stability of a large nonlinear
system can be determined without obtaining the solution of the


8
system of differential equations, that is, the stability limit
can be calculated directly from the system equations and used
during system operation.
In the subsequent chapters, the close relationship of
three different methods of analyzing stability will be discussed.
These are the classical equal-area criterion, the phase plane
trajectories technique, and the second method of Liapunov.
These will be demonstrated for the case of two connected finite
machines. Of the three methods, it will be shown that the second
method of Liapunov is exact for a two machine system but at
present conservative for three or more machines. These will be
carried out for an example involving the transient stability
region for the system using the three different methods.


CHAPTER III
CLASSICAL MODEL
3.1 Introduction
Transient stability studies are performed to determine if
a system will remain in synchronism subsequent to the occurrence
of a major disturbance. That is, will the machine rotors remain
in synchronism and will they return to a constant speed of
operation following the disturbance.
3.2 Classical Stability Study Assumption
The basic assumptions usually made in a classical stabil-
ity study are :
1. Mechanical power input to each machine is held constant.
2. Damping is neglected.
3. Loads are represented by passive impedance.
4. Each machine has a constant-voltage behind transient reactance
5. The mechanical rotor angle of each machine coincides with the
angle of the voltage behind transient reactance.
6. Machine saturation is neglected.
The main reason for making the previous assumptions, is
that the calculations are simplified and, thus, attention can be
focused on the main problem of determining the stability region


10
within a reasonable degree of accuracy. References [1] and [13]
consider the detailed reasoning involved.
3.3 Power System Representation
The equivalent two-machine power system of Figure 3.1
will be used as an example (problem 2.18 [1]) in the comparative
analysis of the stability boundaries for the three methods. The
computations are based on a two-finite machine system model in
which the classical representation is used for each machine.
In all cases resistances are neglected. The initial internal
voltages of the generators, the initial torque angles, and the
load equivalent passive admittance matrix are determined from
the load-flow data. In transient stability studies a load-flow
calculation is first made to obtain system conditions prior to
the disturbance. In this calculation all power system data is
converted to a 100 MVA Base for the two generators, transmission
lines, and system components. The load-flow data are given in
Tables 3.2, 3.3, and 3.5 respectively.
3.4 System Equation Of Motion
During transients in a two-finite machine system, the
system goes through the following stages:
(a) System before the disturbance or fault; prefault system .
(b) System during the disturbance or fault; faulted system .
(c) System after the removal of the disturbance or fault; post-
fault system.


11
BUS1
BUS2
Table 3.1 Generator Data (in p.u to Generator MVA Base)

GENERATOR NO. xd(p.u) X^P.U) H(MW.S/MVA) RATING(MVA)
1 0.28 0.080 5.0 50.0
2 0.250 0.070 4.0 120.0
Table 3.2 Transmission Line Data on 100 MVA Base
(Resistance Neglected)
Line Number 3 4 5 6
X p.u(lOOMVA) 0.08 0.06 0 .06 0.13


12
Table 3.3 Generator Data (in p.u to 100 MVA Base)
GENERATOR NO. X^(P.U) XT(P.U) H(MW.S/MVA) RATING(MVA)
1 0.56 0.16 2.5 0.5
2 0.21 0.05 4.8 1.2
Table 3.4 Load-Flow Data
BUS NO VOLTAGE MAG. ANGLE P.U DEGREE LOAD MW MVAR GENERATOR MW MVAR
1 1.03 0 0 0 30 23.1
2 1.02 -0.5 80 40 100 37.8
3 1.018 -1.0 50 20 0 0
Table 3 5 Load- -Flow Data (on 100 MVA Base)
VOLTAGE LOAD GENERATOR
BUS NO MAG. ANGLE MW MVAR MW MVAR
P.U DEGREE
1 1.03 0 0 .0 0.30 0.23
2 1.02 -0.5 0.80 0.40 1.0 0.37
3 1.018 -1.0 0.50 0.20 0 0


13
The form of the swing equation of the two-finite machine
the system during the previous three stages will be the same as
that given by Eq (3.1), except for the difference in the input
output power from one stage to another.
M
dh
dt2
=Pal2 = Pml2 Pei2
(3.1)
Thus, the swing equation for machine one is :
2
Ml = Pal = Pml Pel
dt
2
The swing equation for machine two is :
M2
d262
dt2
= Pa2 = Pm2 Pe2
(3.2)
(3.3)
The relative torque angle 6 = 61- 62 will be used, because
its value is significant in showing the stability or instability
of the two machine system.
Pal Pa2
d 2S
dt'
Mi
M2
(3.4)
Next, multiply each side of Eq.(3.4) by the equivalent
inertia constant.
M1M2
M =
M
d 2S
Ml + M2
Pal M2
dt'
d 2S
dt2
Ml + M2
Pa2Mi
Ml + M2
M2Pml MiPm2
Ml + M2
M2Pei MlPe2
Ml + M2
(3.5)
(3.6)
(3.7)
Which may be written more simply as in Eq. (3.1)
M
d 2S
dt2
= Pml2 Pel2
(3.8)


14
Where
Pal2= the equivalent acceleration output power
Pml2= the equivalent mechanical input power
Pel2= the equivalent electrical output power
M = the equivalent inertia constant
8 = the relative internal voltage angle between the
two machines(6l £2).
3.5 Prefault Power System Model
The single-line diagram of Figure 3.3 is a further simpli-
fication of the system of Figure 3.2. The elements of the bus
admittance matrix for the network are reduced by delta-wye trans-
formation to two nodes in addition to the reference node.
From the load-flow study of the pretransient network, the
following phasor values are determined:
II W 1 21 e*10 P.U
E2 = 1 15 jl2.65 P.U
II 1 789 + j 0.716 P.U
t-^5 CO II 1 039 + j 0.521 P.U
Pmi= 0 3 P.U
Pm2= 1 0 P.U
matrix for the reduced network is :
0.064 j 1.001 0.182 + j 0813
^matrix =
0.182 + j 0.813 0.528 j 1.358


15
Figure 3.2 Pre-Fault System Reactance Diagram
-0. 182-J0.813
Figure 3.3 Pre-Fault Admittance Diagram


16
The equation of electrical power output from a machine is given by
n
. = E? 6.. + r E. E. Y.. cos( 0.. - 6. + 6. )
l l li l j lj ij i J
P =
e
So that,
1
Pei=
Pe2=
Pe2=
Since Pe!2
Where
S2 G h
1 li
t G H
2 22
1=1
j*i
1 2 12
1 2 12
PelM2 Pe2Ml
Ml + M2
(3.9)
(3.10)
(3.11)
(3.12)
(3.12)
(3.13)
Ml=23XlO ^sec2./elec deg., M2=45X10 ^sec2./elec deg.
Then Pei2 = -0.177 + 1.135 sin ( 8\2 + 4.1 ) (3.14)
The maximum power that can be transmitted under rated
operating condition occurs at a torque angle of( 812 + 7 ) = -90,
where 7 = 4.1 and 612 = -94.1. Notice that the characteristic
is a sine function of 8 between the two rotors as plotted in Fig-
ure 3.4. This sine curve is displaced from the origin vertically
by an amount Pc, which represents the power dissipation in the
equivalent network, and horizontally by the angle 7, which is the
real component of the transfer admittance Y12.
The system is stable only if the torque angle 8 is in the
range from -90to +90, in which the slope dP/d8 is positive, so
that an increase in torque angle results in an increase in
transmitted power.


REAL POWER Pei2 (in per unit)
PRE-FAULT o Pc=-0.177 p.u.
-3
Figure 3.4 Pre-Fault Power-Angle Curve


18
So the maximum power resulting is
Pe, = -0.177 + 1.135 sin(-94.1) = -1.31 p.u
l2max. *
(3.15)
Pe = -1.31 p.u occurs at S12 = -94.1 ( -1.642 rad.)
l2max *
And the steady state synchronous speed, Pm12= Pe12, *s
^ PmlM2 Pm2Ml ,
Pm, = ---rr---rrr--- = -0.1452 p.u
12 Ml + M2 *
-0.1452 = -0.177 + 1.135 sin ( So +4.1)
arcsin ( So +4.1 ) =0.0318
(3.16)
So = -2.5 (-0.0436 rad.)
3.6 Power System Model During Fault
The single-line diagram for the faulted power system is
shown in Figure 3.5 and 3.6. The disturbance to be considered is
a three-phase short circuit which occurs near bus number 3 at the
end of line 5. The line-relaying schemes will sense the fault on
the line and the fault will be cleared in six cycles by simulta-
neously tripping of the line terminal breakers.
Using the same power output equation given in (3.9), and
the bus admittance matrix for the reduced network during the dis-
turbance the Y matrix for the reduced network is
0.0001- j 1.295 0.001 + j 0.044
0.001 + j 0.044 0.006 j 3.419
The new system equations are:
2
Pe = El Gil + E E Y cos (012 Si + S2)
1 1 a 1 Z
matrix =
(3.17)


19
Figure 3.5 During Fault System Reactance Diagram
-0.00077-J0.044
Figure 3.6 During Fault Admittance Diagram


I
20
Pei= 0.000015 + 0.0606 cos ( 89 0l + 82)
COS ( 012 -82+81 )
(3.18)
(3.19)
Pe = 0.0079 + 0.0606 cos ( 89- 82 + 81 )
L
(3.20)
Since
p PelM2 Pe2Ml
12 = Ml + M2
(3.21)
Then
Pei2 = _0-0026 + -061 sin < 512 + 1- )
(3.22)
3.7 Post Fault Power System Model
The post-fault system single-line diagram before and after
network reduction, is indicated in Figures 3.7 and 3.8. The mag-
nitudes of El and E2 computed for the prefault system operation
condition are used in the postfault, and fault model, to find the
electrical output power Pel2, according to the previous assump-
tions. 01 and 02 are time dependent variables which are changed
for each time interval.
line end breakers, the power-angle equation applies, since a net-
work change has occurred. This requires the use of Eq.(3.9) and
the bus admittance matrix for the reduced network post-disturb-
ance. The Y matrix for the reduced network is :
When the fault is cleared by simultaneous opening of the
0.0667- j 0.999 0.186 + j 0.807
Y
matrix =
0.186 + j 0.807 0.519 j 1.347


21
Figure 3.7 Post-Fault System Reactance Diagram
-0. 185-J0.810
Figure 3.8 Post-Fault Admittance Diagram


22
And the system equations are :
Pe = E2 G + E E Y cos ( 012 5l+ 62 ) (3.23)
1 1 XI 1 1
PeJ= 0.09760 + 1.15 cos ( 77 - 5l + 62) (3.24)
Pe2= E2 22+ E1E2Y12 CS ( 012 ~ Si + 6l ) (3.25)
Pe = 0.6860 + 1.15 cos ( 77- 62 + 61 ) (3.26)
it
Since
p PelM2 Pe2Ml
e!2 Ml + M2
Then
(3.27)
Pe = -0.167 + 1.120 sin ( <512 + 4.3 ) (3.28)
1 it
2 2
The steady-state values of Pe 6 and d 6/dt for the
1 it
stable postdisturbance system are :
Pm12 = Pei2(steady-state) = -0.1452 p.u
-0.1452 = -0.167 + 1.12 sin ( <$ss +4.3 )
arcsin ( 6ss + 4.3 ) = 0.0218
8ss = -3.2 (-0.0558 rad.), and d26ss/dt2 = 0
From the power-angle diagram shown in Figure 5.1.
Jrnax. = 7T <5ss
£max. = -176.8 (- 3.0857 rad.)


CHAPTER IV
NUMERICAL SOLUTION OF THE SWING EQUATION
4.1 Introduction
The analysis of first-swing classical transient stability
constitute the important tools for judging system performance.
The reason for its relative importance is that if the system is
stable on first swing, it will for most cases be stable on the
subsequent swings. The torque angle 8 is calculated as a function
of time over a period long enough to determine whether 6 will
increase without limit or reach a maximum and start to decrease.
4.2 Point-by-Point Solution of the Swing Equation
The swing equation governing the motion of each machine of
the power system, is
M = Pa (4.1)
dt2
Where
8= displacement angle of rotor with respect to a reference
axis rotating at synchronous speed
M= inertia constant of machine
Pa= accelerating output power
t = time


24
The solution of Eq.(4.1) gives 6 as function of t, which is known
as "swing curve". Examining the swing curve of the two machine
system will show whether the machines will remain in synchronism
after a disturbance.
In a two-finite machine system, the output and hence the
acceleration power of each machine depend upon the angular posi-
tion 6 and d£/dt. Thus, for a two machine system, the two simul-
taneous differential equations are:
2
Ml = Pml Pel( Si, 62 )
dt2
2
M2 ---2 = Pm2 Pe2( 61, 62)
dt2
(4.2)
(4.3)
Analytical solution of these equations are not possible and
solved numerically. The method is called step-by-step or point-
by-point solution.
In deriving the swing curve of a system under study, step-
by-step solution becomes very handy comparing to some of the
methods recommended for digital computers. In this method, the
change in the angular position of the rotor during a short
interval of time is computed by making the following assumptions:
1. The accelerating power(Pa) computed at the middle
instead of the beginning of the interval to eliminate the error
caused by Pa which is constant during a time interval At.
2. The angular velocity (w)is constant throughout any
interval at a value computed for the middle of the interval.


25
These assumptions are justified by decreasing the time interval;
as the time interval is decreased, the computed swing curve ap-
proaches the ture curve. Figure 4.1 demonstrates the above
assumptions.
If Eq.(4.1) is divided by M and integrated twice with respect to
t, Pa being treated as a constant, we obtain the following
d£ dt = W = wo Pa + M * (4.4)
O II + uo , Pa .2 + M 1 (4.5)
Then A W(n ~ 1/2 A t >- M Pa and wen - 1 /2 ) = W(n 3/2) + A W(n -1/2) (4.7)
also A 8n = At. w A <$n = At. (At)2 W(n 3/2)+ jj Pa A 8n = A 8( 2 n l ) + ^ ^ Pa(n l) (4.10)
6n = 6 Eq.(4.10) is the main form of the step-by-step solution of the
swing equation. The equation shows how to calculate the change
in 8 during an interval if the change in 8 for the previous
interval and the accelerating power for the interval in question
are known. The accelerating power is calculated at the middle of
each new interval. The solution progresses through enough
interval to obtain points for plotting the swing curve, which is
sufficiently accurate[2].


26
Pa
and 6 as Function of Time


27
4.3 Numerical Methods for Solutions of Differential Equations
The methods most commonly used for the solution of the
differential equations are: Euler method, the modified-Euler
method, Runge-Kutta and the trapezoidal method. Each of these has
advantages and disadvantages which are associated with numerical
stability, time-step size, computational effort per integration
step and accuracy of the obtained solutions.
Appendix A shows the computational printout for plotting
the swing curves of the two-finite machine system. For clearing
at 0.1 second the solution is obtained by use of Turbo Pascal
program with the modified-Euler method procedure. So the
numerical solution of the swing equation for the two generator,
three-bus power system is made by digital computer for 2.0 second
of simlated real time, for the intervals of 0.05 sec. Figure 4.2
shows the rotor angles of the two machine vs. time. A plot of the
rotor angle differences is shown in Figure 4.3 and the fault is
cleared in six cycles. It follows that the system is stable. The
rotor angle difference reach value of 11.67 and then decrease.
This is the value of £12 at t= 0.2sec. Note that the solution is
carried out for three swings to show that the subsequent swings
are not greater than the first so that the system appears to be
stable. But in the case of the angle differences increase
indefinitely, the system is unstable because both machine will
lose synchronism.


TORQUE ANGLE 5= .degree
Cycles
Figure 4.2 Plot of (Si and 62 vs. Time
to
CO


TORQUE ANGLE 8 .degree
Cycles
Figure 4.3 Plot of Delta Differences (61-62) vs. Time
to


CHAPTER V
STABILITY DOMAIN
5.1 Introduction
As mentioned earlier, there are at least three different
methods which are commonly used to study the stability of a power
system during a fault. The first method is out-lined below, and
the others are discussed in the next two chapters. Appendix B
contains further development.
5.2 Power-angle Characteristics and the Equal-area Criterion
For the two-finite machine system under study, the power
angle characteristics are shown in Figure 5.1, for condition
before, during, and after a three phase fault. The horizontal
line denoted by P=Pml2 represents the equivalent mechanical power
input to the machine. Before occurrence of the fault, the two
machine were operating at synchronous speed with a rotor angle at
t=0 is 6o = -2.5 degree as indicated by the intersection of Pml2
with the prefault curve, this operating point(5o ,Peo) is des-
ignated by the letter a. Once the 3-phase fault has occurred, the
electrical power Pel2 of the system increases to a value corres-
ponding to point b, at which Pai2 = ( Pel2(^> Pml2), this
results in a decrease in both d£/dt and 6. As 6 continues to


REAL POWER Pei2 (In per unit)
TORQUE ANGLE (in degree)
O PRE. DURING U POST A Pm 12 =-0.145 p.U.
Figure 5.1 Classical Power-Angle Curve


32
decreases the system remains disturbed, the power angle trajectory
moves along the faulted power angle curve form point b toward
point c, the fault is cleared, at which time for this case, Pel2
decreases to 0.85 p.u ( Sc = -143 degrees as indicated by point
d in the figure). The area Al in the figure is proportional to
the kinetic energy (K.E.) of the system during the period when the
fault occurs and is cleared. As S continues to decrease, the
power angle trajectory moves along the postfaulted power angle
sine curve from point d toward point e for which 5max. is -176.8
degrees. When this point is reached,
O
= SQR 2 Pal 2 d£
dt l M J
= 0
(5.1)
At this point the K.E. of the equivalent finite machine has re-
turned to its prefaulted value, A K.E = 0 .
At point e, the area A2 bounded by Pel2 = Pml2 and
Pen = [ Pc + Pmax. sin (S + 7(post)) ], Sc < 6 <6 max.,
is equal to the area Al bounded by Pel2 = Pml2 and
Pei2= [Pc + Pmax. sin(5 + 7(during))] So < S <6c .
There exists a (6=5c), ( So < S < Smax ), such that in
terms of the equal-area criterion of stability, Area Al is equal
to A2. and the system is stable, the value of Sc is determined
as follows:


33
Ai =
60
n
[Pmi2 (Pc + Pmax. sin( 6 + 7(during))] d6 (5.2)
Let x = 6 + 7 ^ dx = d6
n(6c + 7)
Ai =
[-0.1452 + 0.0026 0.061 sin x ] dx
(5.3)
(So +7)
Al= -0.1426 Sc 0.0547 + 0.061 cos ( Sc + 0.33 ) (5.4)
p 6max.
A2 =
6c
A2 =
[(Pc + Pmax. sin ( S + 7(post>)) Pmi2] d6 (5.5)
p(6max.+ 7 )
[ -0.167 + 1.12 sin x + 0.1452 ] dx
(5.6)
(6c+ 7>J
A2 = 0.0218 6c + 1.177 + 1.12 cos ( 6c + 4.3 ) (5.7)
Al = A2
-0.1644 6c 1.056 cos 6c + 0.084 sin 6c 1.232 = 0 (5.8)
From the power-angle diagram shown in Figure (5.1) the
critical clearing angle is located between 6o and 6max. Since
Eq.(5.8) is nonlinear, 6c= -143 degrees (-2.495 rad ) was found by
using trial and error.


34
If the fault clearing is delayed long enough so that the
equality of the two area cannot be satisfied, the two-finite
machine speed will not decrease to a synchronous value as long as
the machines remain electrically tied to each other. The torque
angle 6 will decrease monotonically without bound beyond the
maximum value possible for a marginally stable swing, 5max. By
the time 6 reaches a value of -180 degree, synchronism is lost
and the machines must be disconnected[5].


CHAPTER VI
PHASE PLANE TRAJECTORIES AND THE STABILITY BOUNDARY
6.1 Introduction
This technique provides a useful tool for studying the
stability of a system which is described by a second-order differ-
ential equation or a group of such systems. The swing equation of
the power system ( prefault, during fault, and postfault) again is
a second-order differential equation of the form:
into two first-order differential equations with the time t
suppressed.
(6.1)
To form the phase plane, the above equation is converted
Let Xl = X2
(6.2)
X2
Pml2 Pel2
(6.3)
M
Where Xl = 8 and (f> = 6ss + 7
dt
2 2
From the steady-state values of Pel2, 6 and d 6/dt the
stable postdisturbance system is
Pm = Pe (steady-state)
1a 1 a
Pe = Pm = Pc + Pmax. sin (5ss + 7 )
1 A 1a
(6.4)


36
c r ' ( Pml2 Pc v -i ,
oss =[ arcsm ( ------------) 7 J (6.5)
Pmax.
Sss = -3.2 (-0.0558 rad.), and d2£ss/dt2 = 0 ( as found
previously for the post-fault power system model ).
6.2 Phase Plane Trajectory
A phase plane plot is a plot of X2 vs. Xl, as shown in
Figure 6.1. As the time t increases, the two-tuple [Xi(t), X2(t)]
describes the trajectory in the phase plane. Since stability of
the postfault system is the basic requirement, finding the state
equilibrium point ( the origin of the phase plane) is important,
since this point is the steady state value.
The state equation for the postfault system is :
Xl = X2 (6.6)
X2 = [ Pml2 Pc Pmax. sin (xl + 0 ) ]/ M (6.7)
Since
Xl = X2 => S ( xl ,x2 ) = x2 (6.8)
X2 = [ Pml2 Pc Pmax. sin ( xl + ) ]/ M =>
Q( xl,x2 ) = [Pml2 Pc Pmax. sin( xl + )/ M (6.9)
Then the state equilibrium singular point is:
S( xl ,x2 ) = 0 => x2 = 0 (6.10)
Q( xl,x2 ) =[Pml2 Pc Pmax. sin(xi + )/ M = 0 (6.11)
( xl + (j> ) = [ arcsin ( 12-----^) ] = n7r (6.12)
Pmax.
where n=0, 1, 2, .........


State Variable X2 In degree/sec.
(Thousands)
8
-190 -170 -150 -130 -110 -90 -70 -50 -30 -10 10
State Variable X1, degree
Figure 6.1 Phase Plane Plot Of X2 Vs. Xl
co
-o


38
the singular point is chosen for n = 0, which defines the state
variables, Xl and X2, such that the origin of the phase plane is
the stable equilibrium point of the postfault system.
The equilibrium states of a system are defined as:
dS dS
dX, dXn
1 2
dQ dQ
dX, dX
L i 2 J
It is necessary to evaluate the equilibrium states of the system at
singular points. This is done by testing the stability of the
system by determining the nature of the roots. If positive real
roots exist, the system is unstable for the given operation
conditions, if no positive real roots exist the system is stable.
0
1
Pmax. cos 6
0
J(0,0)
Now, the stability of the system is from
-A
Pmax.
1
= A2 Pmax. => A = VPmax.
"A J(0,0)
The result is two real roots with opposite signs, then point e on
the phase plane is an unstable singular point being a saddle
point singularity of the trajectories.
The state equations for each system are :


39
Pre-fault system
Xl = X2
(6.13)
X2 =[Pml2 (Pc + Pmax. sin(5 + 5ss + 7 ))]/M (6.14)

During fault system
Xl = X2
(6.15)
X2 =[Pml2 (Pc + Pmax. sin(5 + 5ss + (6.16)
Post-fault system
Xl = X2
(6.17)
X2 =[Pml2 (Pc + Pmax. sin(5 + 5ss +7 , ))]/M (6.18)
(post)
These equations are plotted in the phase plane with torque
angle 6 at its minimum value (5min = 19 degrees as indicated by
point f in Figure 6.2). Figure (6.2) illustrates the relationship
between power-angle and the phase plane trajectories for a
marginally stable case. The prefault operating point (5o,Peo) is
designated by the letter a in the figure, the power system is in a
state of equilibrium with (8= d5/dt = 0), this point in the phase
plane is a stable singular point called a stable node or a vortex.
When the 3-phase fault occurs at t = t*, the electrical output
power of the system Pel2 increases to a value corresponding to
point b, this results in a decrease in the speed deviation d5/dt
and the torque angle 8 This decrease is depicted in the phase
plane fault trajectory between 60 and 8c. In the marginally
stable case, fault clearing is delayed long enough to permit the
system torque angle 8 to decrease to the critical value 5c at


Stcrt* Variable X2 In
40
Figure 6.2 Power-Angle Curves And Phase Plane Trajectories


41
which time the faulted line is isolated from the system. For
fault clearing torque angles greater than 6c, the postfault system
will be unstable. At the time the fault is cleared (at t=0.1
sec. corresponding to 6=6c) the output power Pel2 decreases from
the value of -0.5 to -0.85 p.u, 6 continues to decrease along the
phase plane post-fault trajectories from point c and d (at which
6 = 6c) toward point e. Point e (6 = 6max.) is an unstable
singular point which is a saddle point of the trajectories. For
the conservative system under study, as 6 increases a path is then
formed along the phase plane maximum trajectory in a clockwise
direction from point e to point f as show in the Figure.
Referring to Figure (6.2), tfmin can be determined from the
equal-area criterion, using postfault power-angle curve such that:
Jss
[Pml2 ( Pc + Pmax. sin(6 +7 , )] d 6 (6.19)
(post)
6min
£max
[ (Pc + Pmax. sin(6 +7 , )
(post)
Pml2] d 6 (6.20)
Let x = 6 + 7 ^ dx = d<$
(£ss + 7 )
(5max + 7 >
[-0.167 + 1.12 sin x + 0.1452]d x = -2.16 p.u (6.21)
6min = 19 ( 0.332 rad.)
(6.22)


42
The critical clearing time t can be determined from
6 cr
Eq.(5.1), using numerical integration techniques:
t =
cr
So
t =
cr
So
Sc
d S
SQR
Sc
[ T ioj PaI2 iS ]
dS
SQR [ J [ 1,11112 "Pel2] dS ]
Where
(6.23)
(6.24)
Pel 2 = [ Pc + Pmax. sin (5 + 7, )]
(dur)
(6.25)
So
Pal 2 d5 =
[Pmi2 (Pc + Pmax. sin(5 + 7)]d5 (6.26)
JSc
t =
cr
So
d S
(6.27)
SQR[-§-[(Pini2 Pc)(5-5o) + Pmax.(cos5 cos5o)]
t = 0.1 second
cr
6.3 Stability Boundary
The system is considered stable as long as the trajec-
tories follow the separatrix determined by the phase plane (as
shown in Figure 6.1). The equation for the separatrix can be
determined from the post-fault system differential equation of
motion :


43
d 2S
Since
Then
dt2
d r d($i
dt dt.
d. r d5
dt
1
Pal2 = Pml2 Pei2
= 2
d 8
dt
dt
d<5
dt
1 2
= 2 d 6
d 6
d 2S
dt2
d 6
dt
d 2S
dt2
Substituting Eq.(6.30) for d26/dt2 in Eq.(6.28), the
differential one-form is obtained.
c-i 2
d 6
d 6
dt
d 6
2
dt
[Pml2 (Pc + Pmax. sin 5)] d5
(6.28)
(6.29)
(6.30)
= [Pml2 (Pc + Pmax. sin 5)] (6.31)
(6.32)
In term of the state variable, Xl and X2,
d(X2) = [Pml2 (Pc + Pmax. sin(Xl + ^))] dXl (6.33)
h 1f
Integration is performed to obtain the general solution for a
postfault system phase plane trajectory:
X2 [Pml2(Xl) Pc(Xl) + Pmax. cos(Xl + ^)] + C = 0
2 M
The constant of integration C is evaluated at the saddle point of
the separatrix, which is ( 6max 0 ).
C= ---t(Pml2 Pc)(6max ) + Pmax. cos 6max] (6.34)


44
The stability boundary of the postfault system, is
X2 [ Pmi2(Xl) Pc(Xi) + Pmax. cos( Xl + ) ]
Z t|
and
+
2
[(Pml2 Pc)(<$max ) + Pmax. cos 6max] = 0
M
(6.35)
X2 + ^ ["( Pral2 Pc )( tfmax Xl )
2 M L
+ Pmax. [cos 6max cos (Xl + <£)j = 0
( tfmin (j> ) < Xl < ( £max )
(6.36)
(6.37)
The results obtained by this method agrees with that deter-
mined by the equal-area method.
The postfault system trajectories illustrated in Figure
(6.3), were generated by varying the fault clearing time (tfc)-
If the fault is cleared in time at which t < t operation will
fc cr
be stable along a curve as indicated in the figure for the three
stable phase plane trajectories, being confined to the phase plane
region enclosed by the stability boundary. For a sustained fault
or for a longer delay in clearing time at which the
faulted trajectory will enter the unstable region as shown in the
figure for the two unstable phase plane trajectories by the dotted
curve.
Since the power system under investigation is a conserva-
tive system (no damping), the stable trajectories are closed tra-
jectories. For stable trajectories in a non conservation system
X(t) > 0 as t oo.


(0
CO
Figure 6.3 Plot of Stable and Unstable Phase Plane Trajectories


CHAPTER VII
ASYMPTOTIC STABILITY
7.1 Introduction
Liapunovs "second"(also called the "direct") method
has been used with varying degrees of success in analyzing both
linear and non-linear stability criteria. A desirable feature
of the second method is that it will permit the determination of
system stability without obtaining solution of the system
differential equations.
7.2 Liapunov stability
The Liapunov second method permits stability analysis of
system differential equation, of the form :
X. = Fi (X) 1 < i < n (7.1)
Where X is the postfault system state vector ( x ,x ,..., x )
12 n
determined from the differential equations of motion which
describe the behavior of the power system. There are four
conditions which need to be studied to analyse the stability of
a system :
1. V (X) and its first partials are continuous in R, where R
is a region about the origin, and V(X) is a scalar function.
2. V(0) = 0


47
3. V(X) >0 for all X G R, X 0
4. V(X) < 0 for all X G R, X 0
If the four conditions are satisfied, then Eq.(7.1) is stable.
If condition 4 is replaced by :
V(X) <0 for all X G R X^O,
then Eq. (7.1) is asymptotically stable in the region R If
there exists a scalar function V(X) which satisfies the above
conditions for all X in state space, then Eq.(7.1) is globally
asymptotically stable. If a Liapunov function does not exists
in R, the system is unstable.
For a particular system, there is no general method to
determine the existence of a suitable Liapunov function. A useful
method for determining the V-function for the system under study
is the Energy Metric Algorithm(see appendix C).
7.3 Energy Metric Algorithm
Applying this method to a second-order system( n = 2), we
define the state variables, as before, such that X (0) = £ .
Xi = 8 and = 8ss + 7
x2=^-
dt
The resulting two first-order differential equations are the
same as those developed in previous chapter (Eq.6.17 and 6.18).
Let Xl = X2 (7.2)
X2 =[Pml2 (Pc + Pmax. sin(6 + 6ss +7 . ))]/M (7.3)


48
Since n = 2, there is only one equation thus formed.
dXl
= M
X2
dX2 [ Pml2 (( Pc + Pmax. sin( Xl + 0 )) ]
After the cross multiplication and rearrangement of the terms,
there results:
(7.4)
HI = X2 dX2 + --[Pc + Pmax. sin(Xl+ ) Pml2] dXl (7.4)
A line integration is preformed to obtain the V-function.
V=
HI
n VV0)
V=
ou
n x2(Xr x1)
[Pc + Pmax. sin(Xl + ) Pml2] dXi+
X2 dX2 (7.5)
V= fpc(Xl) + Pmax. [cos cos(Xl + )~\- Pml2(Xl)l
M L J
1 2
+ X2 (7.6)
2
The system is stable for the equilibrium point under
consideration if the total derivative of V with respect to time
is at least negative semi-definite. ( A function is negative
semi-definite if it is equal to zero at the origin and less
than or equal to zero in a region about the origin).
V(X) = Xi + X2 (7.7)
dXl 3X2
V(X)= [Pc-Pml2 + Pmax. sin(Xl + ^)] Xl+ X2 X2
(7.8)


49
Since Xl= X2 and V= 0
Then V(X)= [Pc-Pml2 + Pmax. sin(Xl + <£)](Xl- X2) (7.9)
M
Note that V(X), and V(X) are continuous for all X in state space.
Furthermore, V(X) < 0 for all X 6 R and V(£) = 0.
The region of stability is found to be:
R : X e L V(X) < X ] (7.10)
Where
X = V(X)I r A n
1 ( xj= dmax

X = | Pc( tfmax -)- Pml2( 6max (f> )
M L
+ Pmax. ( cos (f> cos 5max )j+ 0 (7.11)
The stability boundary of the postfault system, is
V(X) X = 0
2 2 r
X2 + --(Pc Pml2)(Xl) + Pmax.
M L
(7.12)
[cos cos(Xl + ^)]j
j (Pc Pml2)(5max (j>) + Pmax. (cos (f> cos 6max)l= 0
M L J
X2 + ^ [ ( Pml2 Pc )(6max tf> Xl)
M L
+ Pmax. [ cos 5max cos ( Xl + ) ]J = 0
( £min (j) ) < Xl < ( 5max )
(7.13)
(7.14)
This result is identical to that obtained using the phase
plane technique, in the following sense the maximum trajectory,
which is the stability boundary for the system, has been found by


50
setting X2 = 0 and Xl = ( tfmax ) in Eq.(7.10). The result
defines the region of stability for the equilibrium at
6= Smax (saddle point). This saddle point determines the
value of V for the maximum stable trajectory, since the stability
boundary passes through the singularity. The same result was
found using Eq.(6.36) which was developed by the phase plane
approach.
Eq.(7.14) defines the exact boundary between stable and
unstable regions of the power system operation. Also from
Eq.(7.12), which describes the closed surface with the greatest X
is the value of the Liapunov function V(X) on which boundary of
the largest stability estimate. The value of X corresponds to the
shaded area as shown in Figure(7.1), is the maximum amount of
stored potential energy in the post-fault system, or a given power
system operating condition. The value of the constant X can be
found using postfault system power-angle curve
a. Shaded Area
X = ------------ or
X =
6 ss'
n <$max.
1
[Pc + Pmax. sin(6ss + 7 - Pml2] d6 (7.15)
(post;
Stability domains corresponding to the smaller of the
values of Pml2, the larger will be the value of X, and the larger
will be the phase plane region of stability. As shown in Figure
(7.2) for a given postdisturbance system condition.


REAL POWER Pei2 (in per unit)
POST TORQUE ANGLE (in degree) A Pml2 =-0.145 p.U.
Figure 7.1 Determination of the Stability Boundary V-Function
Value, K, form the Postfault System Power-Angle Curve


25
Figure 7.2 Relationship Between X and the Mechanical
Shaft Power Pm12


53
7.4 A Suitable Liapunov Function and an Estimation of X
As indicated previously, the estimated X is the value of
the Liapunov function V(X) on the stability boundary. From
Eq.(7.6) V(X) is :
V(X) = |pc(Xl) + Pmax.tcos cos(Xl + )]
M L
1 1 2
- Pml2(Xl) + X2 (7.16)
J 2
Eq.(7.16) is positive definite, if
(Pc Pml2)(Xi) + Pmax.[cos cos(Xl + )] > 0 (7.17)
setting, V(X)|(x =Q^ x =Q) in Eq.(7.16).
1 2
V(0,0) = Pmax. (cos -cos ) = 0 ^ a minimum
From Eq.(7.16), V has a turning point at a critical point of the
power system equation :
Next set, V(X)I o An
1 (Xj =6max.- g>, x2=0)
V(^max. -(f), 0 )=[ (Pc Pml2 )( £max.- (f> )
+Pmax. (cos -cos 5max.)] = Xmax. > 0 (7.18)
Therefore, the stability boundary (surface) of the postfault
system, given by a Liapunov function is :
V(X) Xmax. =0 or
2 2 r
X2 + --- (Pml2 Pc)(6max Xl)
M L
+ Pmax.tcos 6max cos (Xl + ) ] = 0, provided that
surface V = X = constant.
(7.19)


54
Xmax. > X > 0, are closed surfaces, since the Liapunov theorem
is satisfied within the region given by Eq.(7.14). Also,
V(X) = [Pc Pml2 + Pmax. sin(Xi + ^)](Xi - X2) is the
M
solution for the power system of Eq.(7.2 and 7.3). The region is
asympotically stable. Closedness can be proved by considering the
behavior of the surface V = X .
Let Xi= C = a constant, Eq.(7.16) becomes
2
V=X=^f(Pc Pmi2)(C) + Pmax. [cos cos(C + $)]] + ^-X2
M L J 2
1 2
V= X2 + f(C) (7.20)
2
Then the equation becomes
1 2
X2 = X f(C) (7.21)
2
Eq.(7.21) describes the intersection of the surfaces V=X
with the plane Xl= C. For any Xmax. > X )> 0, as C increased
from zero, to the right hand side of Eq.(7.21) is initially
positive, and the intersections are closed ellipsoidal curves.
As C increases, the plane Xi= X moves along the negative Xl-axis,
and [X f(C)] will eventually become zero, and the curve of
intersection vanishes, which indicates the closedness of the
surface V=X.
For the case when X= Xmax. closure is on the Xl-axis for
f(C)= Xmax i.e., when[X f(C)]= 0, or


55
[(Pc Pml2)(C) + Pmax.Lcos cos(C + 0)]1 = 0 (7.22)
M L J
The two roots of Eq.(7.22) give two points:
(a) Principal node or vortex point (0,0)
(b) Between the origin and the complementary saddle point
(6max. -(f), 0).
Therefore, Eq.(7.19) defines a useful Liapunov stability bound
and an operating limit, which produces the largest possible
stability domain. It should be emphasized that the application of
the Liapunov method depends strongly on the choice of the Liapunov
function. However since the Energy Metric Algorithm is unique in
power system stability studies for generating a Liapunov function
directly from the system equations, the search is made somewhat
easier and if desired this search could be performed with the aid
of a digital computer.
Figure 7.3 illustrates the effect of the equivalent
mechanical power output of the shaft for the prefault system
loading (Pml2) on the phase plane postfault system stability
domain. Stability domains corresponding to the smaller the value
of Pml2, the greater will be the region of phase plane enclosed by
the stability boundary.


C4
o
to
-H
CO
Figure 7.3 Effect of the Equivalent Mechanical Power (Pmi2)
on the Phase Plane Postfault Stability Boundary
C7l
05


c
CHAPTER VII
MANSOURS MODIFICATION OF THE EMA
8.1 Introduction
The direct method of Liapunov has been tried in many dif-
ferent types of applications with varying degrees of success.
In particular, the Energy Metric Algorithm was applied by Mansour
[10, 11] to power systems with great success. In analyzing the
transient stability of multi-machine power system, Mansour
employed a number of different methods in generating Liapunov
functions. From a comparison of the results produced by each
of the different methods Mansour concluded that the Energy-Metric-
Algorithm with a slight modification produced the largest region
of stability.
8.2 Power System Representation
The equivalent two-finite machine power system model of
Figure 3.1 will be used in the analysis for finding a larger
stability boundary. The system equations of motion are :
(8.1)
Eq= fl(Eq,, E^)
(8.2)
Pm= f2( Pm, 6 )
(8.3)


Where
58
6= Torque angle
M= Inertia constant
D= Damping coefficient
Pm= Mechanical power output of the shaft
Pe= Electrical input power
g
ex= Exciter voltage
Eq=Voltage proportional to field flux linkage
* ^
Eq=Flux decay
Pm= Governor action 43
6= Angular speed deviation
8.3 System Equation of Motion
The power system is simplified for convenience by assuming
constant (independent of 6) uniform ( D/M ratio is the same for
each machine) damping, no transient saliency, no flux decay and
negligible governor action. This will yield the postfault system
model described by the state Eq.(6.17) and (6.18), with constant
uniform damping introduced into the system model. A more detailed
study of the power system transient stability problem was recently
accomplished by Roemish[5]. The system equations become:
----= Pm12 Pel2
dt
(8.4)
Where
Pei2= Pc + Pmax. sin (6ss + 7)
(8.5)


59
The system equations in state-variable form are
Xi = 6 and = 6ss + 7
x2
dt
Let Xl = X2
X2M= -DX2+ [Pml2 (Pc + Pmax. sin (Xi+ (f>))]
(8.6)
(8.7)
8.4 The unmodified EMA
The Energy Metric Algorithm is applied to the system
under study to determine a suitable Liapunov candidate function
(V-function).
dXl __________________________X2___________________ m
dX2 D X2 + [Pmi2 (Pc + Pmax. sin( Xl + ))] ' J
Since n = 2, there is only one equation formed.
After the cross multiplication and rearrangement of the terms,
there results :
HI = D X2 dXl+ [ Pc + Pmax. sin ( Xl + ) Pml2 ] dXl
+ M X2 dX2 (8.9)
Next, a line integration is performed to obtain the V-function.
ft
V(X) =
HI
U
, xl(x2 = 0,
V(X)= D X2 dXl +[ Pc + Pmax.sin( Xl + ) Pmi2 ] dXl
o*
n
vw
+
M X2 dX2
ou
(8.10)


60
V(X)= £ Pc(Xl) + Pmax.tcos (j> cos(Xl + )] Pml2(Xl)J
i __ j
Periodic Terms
2
+ 0.5 M X2 (8.11)
J
Quadratic Term
The Liapunov function [V(X)] obtained by applying the EMA
to the system equations of a two-finite machine power system
produces an incomplete quadratic form given by Eq.(8.11). By
completing the square of the quadratic form, the Liapunov function
candidate results in a more accuratelly defined stability region
for the power system under investigation than can be obtained by
assuming a totally conservative system.
8.5 Mansour*s Modified EMA
An improvement in the stability domain within which a
power system can operate means a better response more accuracy.
This means that under load and fault conditions the system can be
maintained in synchronism and thus avoid the need for emergency
load dropping. The missing terms from the candidate function can
be supplied by examining Eq.(8.11). Considering only the
quadratic terms [denoted as Vq(X)] of the Liapunov function
presented in the equation, which means ignoring the periodic terms
in brackets for simplicity. Mansour proposed adding the terms aX?
and X1X2 to the Liapunov function originally determined from


61
the system equation. The resulting generalization of Eq.(8.11) is
2 2 r
V(X)= 0.5 M X2 + /? X1X2 + a Xl + |_[ Pc Pml2](Xl)
+ Pmax.t cos cos ( Xl + )]]
(8.12)
2 2
Where Vq(X)= 0.5 M X2 + /? X1X2 + a Xl
From Liapunovs theorem it is necessary that the
V-function to be positive definite for stability. The values of a
and /3 for which the quadratic form is positive definite can be
determined by examining the Hessian matrix of second partials.
32V 32V
3X? 3X2X1
= 32V 32V
3XiX2 3x1
A real and symmetric matrix Vt can be formed from the coefficients
of the quadratic form of
Vq(X)= XTJf X
(8.13)
Vq(X)= |xi X2]
a2v
3x1x2
32V
3X? 3X2Xi
32V 32V
3X2
Xi
X2
(8.14)
For the function V to be positive definite the following require-
ments must be satisfied
> 0 (8.15)
3X?
det 51^0
(8.16)


62
Evaluation of the Hes sian matrix yields:
dv = 2 aXi d2V r = 2 a (8.17)
dXl dX?
dv = M X2 + P Xl =*> d2V = M (8.18)
dX2 dX§
d2v d2V p (8.19)
dXlX2 dX2Xl
2 a p
= (8.20)
P M
Since det If y 0 4- a > 0 (8.21)
From this there results
det Jf = 2 a H P2 (8.22)
and a > - P2 2M (8.23)
2 2
Since Vq(X)= a Xl + P X1X2 + 0 5 M X2 (8.24)
The total derivative of Vq(X) with respect to time is
Vq(X)= -& Xl + -9va_ X2 (8.25)
dXl dX2
Vq(X) = ( 2a Xl + /? X 2 ) Xl + ( M X2 + p Xl )X2 (8.26)
Since Xl = X2 and from Eq. (8.7) X2 =------------------^7- X2
M
Vq(X) =(2a Xl + fi X2)X2 + (M X2 + /? Xl)(- -|[-)X2 (8.27)
2
Vq(X) = ( /? D ) X2 + ( 2a )Xl X2 (8.28)
a is selected so that V is negative definite. This can be


63
0 £ D (8.29)
2
Vq(X) = ( /? D ) X2 (8.30)
A resulting Liapunov function is then obtained by the sum of the
periodic and quadratic terms to give
V v + V
total quadratic periodic
2 2 r
V(X)= 0.5 M X2 + /? X1X2 + a Xl + Lt Pc Pml2](Xl)
+ Pmax. [ cos - cos ( Xl + )]] (8.31)
accomplished if
a = 0
2M
Therefore
Since a = , from Eq. (8.29)
V(X)= 0.5 M X2 + 0 XlX2 + -|jj- Xl + [[ Pc Pml2](Xl)
(8.32)
+ Pmax. [ cos (j) cos ( Xl + 0 )]]
The time derivative of ^aj(X) is given by
V(X) = [2-5- Xl + /?X2 + [Pc Pmi2 + Pmax. sin(Xl + <£)]lxi
L 2M J
+ [ /? Xl + M X2 ]X2
(8.33)
Since
Xl = X2
and X2 =
D
X2
V(X) = [0- Xl + /?X2 + [Pc Pml2 + Pmax. sin(Xl + <£)]1
L 2M J
X2
+ [ p Xl + M X2 ] -J- X2 (8.34)
M
2
V(X)=(/9 -D)X2 + [Pc Pml2 + Pmax. sin(Xl + ^)]X2 (8.35)
By varying the parameter /? the largest stability region within
a reasonable degree of accuracy can be determined.


CHAPTER IX
CONCLUSIONS
The resulting knowledge gained from a two-finite machine
system which has been studied in this thesis, allows a number of
general conclusions to be drawn concerning the effect on stability
of certain concepts used in power system design, apparatus design,
and power system operation. The effect of system modifications
must be analytically observed before a fault, during a fault, and
after fault clearance. Experience has shown that some design
changes improve stability during all three conditions, while other
modifications are helpful during one condition and detrimental
during other changes.
Out of the three methods studied in this thesis it is
found that the second method of the Liapunov is exact for two
machines, but conservative for three or more machines. This is
one advantage of this method over the other two methods( equal-
area criteria and phase-plane trajectory). These two methods,
however, are suitable for a two-machine system.
The methods developed and used in this thesis provide an
on-line digital computer approach to the power system stability
problem, which is functional and convenient. Because of the


65
complexity of the problem and the constant advancement in the
state-of-the-art, there is still much further research to be
considered in this area. For instance, for complete automatic
operation and security maintenance, greater precision in the
development of the largest stability region for a given system
must be considered. Such research requires the development of the
necessary as well as the sufficient conditions for specific
Liapunov functions.


BIBLIOGRAPHY
1. Anderson, P.M., and A.A. Fouad, Power System Control and
Stability, Vol.l. The Iowa State University Press, Ames.
Iowa, 1977.
2. Stevenson, D., Elements of Power System Analysis,
4th Edition McGraw-Hill, New York, 1982.
3. Gless, G.E., "The Direct Method Of Liapunov Applied To
Transient Power System". IEEE Transactions, PAS-85, 1966.
4. Wall, E.T.,"A Topological Approach To The Generation Of
Liapunov Functions." Acta Technica Csav, Prague,
Czachoslovakia, No.2, 1968.
5. Roemish, W.R., A New Synchronous Generator Out-Of-Step
Relay Scheme. Ph.D. Dissertation, University Of Colorado At
Boulder, 1981.
6. Meyer, J.C., Modern Control Principles And Application,
McGraw-Hill New York, 1968.
7. Jordan, D.W., and P.Smith, Non Linear Ordinary Differential
Equations 2nd edition. Oxfored University Press, New York,
1987.
8. Ogata, K., State-Space Analysis Of Control Systems,
Prentice-Hall, Englewood Cliffs, N.J. 1967.


67
9. Siddiqee, M.W., "Transient Stability of an A-C Generator
by Liapunovs Direct method." Int. J. Control Vol. 8, No.2
1968, PP.131-144.
10. Mansour, M., "Stability Analysis and Control of Power
System." in Real Time Control of Electric Power System,
Edmund Hardschin, Ed. Elsevier, New York 1972, PP. 169-190
11. Mansour, M., "Generalized Liapunov Functional for Power
Systems." IEEE Transactions on Automatic Control, Vol.
AC-19, PP. 247-248, 1974.
13. Kimbark, E.W., Power System Stability, Vol.l. Wiley,
New York, 1956.
14. Lewis, W.A., A Basic Analysis of Synchronous Machines,
Pt. 1. AIEE Transaction, PAS-77, PP. 436-55, 1958.
15. Adkins, Bernard, General Theory of Electrical Machines.
Chapman and Hall, Ltd. London, 1959.
16. Crary, B. Selden, Power System Stability, Vol.l. John
Wiley and sons, Inc. New York, 1945.
17. Wall, E.T., The Generation of Liapunov Function in
Automatic Control Theory by an Energy Metric Algorithm.
Ph.D. Dissertation University of Denver, Denver, Colorado,
1967.
18. Weiss, J.R., "Transient Asymptotic Stability of Power
Systems as Established with Liapunov Functions." IEEE
Transactions, 1976, PP.1480-86.


APPENDIX A
PROGRAM AND FLOWCHART
This appendix contains the main program, and the listing
of computed swing curve coordinates.
A.1 Main Program
This main program "DELTA", written in turbo Pascal, the
program is designed to implement the modified-Euler procedure for
the two generator, three-bus problem.
A.2 Description of the Input Quantities
Name in the Symbol in the Description
program analysis
Pal, Pa2 Pal, Pa2 The accelerating output power in p.u.
OMEGA1, OMEGA2 C4 3 H 3 The machine shaft speeds in rad./sec.
DELTAl, DELTA2 Si, 62 The machine torque angles in degrees.
DELTA_DIFF Sl2=(Sl-S2) The relative torque angle between the two machine in degrees.


69
T t Time.
dt At The time interval
Cl, C2 1/Mi, 1/M2 Inertia constant of the machine in rad./sec .
dOMEGAIJSTART The initial estimation of shaft speed of machine 1.
dOMEGAl_END The final estimation of shaft speed of machine 1.
dOMEGA2_STRAT The initial estimation of shaft speed of machine 2.
dOMEGA2_END The final estimation of shaft speed of machine 2.
dDl_STRAT The initial estimation of torque angle of machine 1
dDl_END The final estimation of torque angle of machine 1
dD2_STRAT The initial estimation of torque angle of machine 2
dD2_END The final estimation of
torque angle of machine 2.


70
A.3 Descriptii
Name in the
program
Pal, Pa2
0MEGA1,0MEGA2
DELTA1, DELTA2
DELTA1_2
of the Output Quantities
Symbol in the
analysis
Pal, Pa2
on, a;2
Si, 82
5l2=(5l-52)
t
Description
The accelerating output
power in p.u.
The machine shaft speeds
in rad./sec.
The machine torque angles
in degrees.
The relative torque angle
between the two machine in
degrees.
Time.
TIME


71
PROGRAM DELTA;
TYPE
LIST_250 = ARRAY [1..250] OF REAL;
VAR I:INTEGER;
d t,D1,D2,Cl,C2:REAL;
dOMEGA1_START,dOMEGA2_START,dOMEGAl_END,dOMEGA2_END:REAL;
dDl_START,dD2_START,dDl_END,dD2_END:REAL;
OUTPUT_FILE:TEXT;
TIME,DELTA1,DELTA2,DELTA_DIFF:LIST_250;
Pal,Pa2,OMEGAl,OMEGA2:LIST_250;
PROCEDURE INITIALIZE;
BEGIN
dt:=0.05;
Cl:=75.3984;
C2:=39.27;
FOR I := 1 TO 250 DO
TIME[I] := dt (I1);
DELTA1[1]:=10.0/57.2985;
DELTA2[1]:=12.65/57.2985;
DELTA_DIFF[1]:=DELTA1[1]-DELTA2[1];
OMEGAl[1]:=377;
OMEGA2[l]:=377;
END;
PROCEDURE COMPUTE_DELTA;
BEGIN
FOR I := 1 TO 110 DO
BEGIN
IF (I<=2) THEN
BEGIN
PaltI]:=0.2999-0.06*SIN(DELTA1[I]-DELTA2[I]+0.01745);
Pa211 I:=0.992+0.06*SIN(DELTA11 IJ-DELTA2[I 10.01745);
END
ELSE
BEGIN
Pal[I]:=0.2024-1.15*SIN(DELTA1[I]-DELTA2[I]+0.22689);
Pa2[I]:=0.314+1.15*SIN(DELTA1[I]-DELTA2[I]-0.22689);
END;
dOMEGAl_START:=C1*Pal[I];
dOMEGA2_START:=C2*Pa2[I];
IF 1=1 THEN
BEGIN
dDl_START:=0.0;
dD2_START:=0.0
END;
IF I>1 THEN
BEGIN
dDl_START:=dDl_END;


72
dD2_START:=dD2_END;
END;
D1:=DELTA1[I]+dDl_START*dt;
D2:=DELTA2[I]+dD2_START*dt;
IF (I<=2) THEN
BEGIN
Pal[1+1]:=0.2999-0.06*SIN(Dl-D2+0.01745);
Pa2[1+1]:=0.992+0.06*SIN(Dl-D2-0.01745);
END
. ELSE
BEGIN
Pal[1+1]:=0.2024-1.15*SIN(Dl-D2+0.22689);
Pa2[1+1]:=0.314+1.15*SIN(Dl-D2-0.22689);
END;
dOMEGAl_END:=C1*Pal[1+1];
dOMEGA2_END:=C2*Pa2[1+1];
OMEGAlt1+1]:=0MEGA1[I]+(dOMEGAl_START+dOMEGAl_END)/2*dt;
OMEGA2C1+1]:=0MEGA2[I]+(dOMEGA2_START+dOMEGA2_END)/2*dt;
dDl_END:=0MEGA1[1+1]-377;
dD2_END:=OMEGA2[I+1]-377;
DELTA1[I+1]:=DELTA1[I] + (dDl_END+Ddl_START)/2 *d t;
DELTA2[1+1]:=DELTA2[I]+(dD2_END+dD2_START)/2*dt;
DELTA_DIFF[1+1]:=DELTA1[1+1]-DELTA2[1+1];
END;
END;
PROCEDURE WRITE_RESULTS;
BEGIN
ASSIGN(OUTPUT_FILE,A:D_OUTPUT);
REWRITE(OUTPUT_FILE);
WRITELN(OUTPUT_FILE,TIME PalPa2 OMEGAl DELTA1 OMEGA2 DELTA2
WRITELN(OUTPUT_FILE);
FOR I := 1 TO 41 DO
WRITELN(OUTPUT_FILE,TIME[I]:10:2,Pal[I]:9:3,Pa2[I]:10:3,
DELTA1[I]*57.2985:10:2,OMEGA2[I]:11:2,
DELTA2[I]*57.2985:9:2,DELTA_DIFF[I]*57.2:10:2);
CLOSE(OUTPUT_FILE);
END;
PROCEDURE MAIN;
BEGIN
INITIALIZE;
COMPUTE_DELTA;
WRITE_RESULTS
END;
BEGIN
MAIN
END.


73
LISTING OF COMPUTED SWING CURVE COORDINATES
TIME Pa! Pa2 OMEGA1 DELTA1 OMEGA2 DELTA2
0. 00 0. 302 0. 988 377. 00 10. 00 377. 00 12. 65
0. 05 0. 303 0. 987 378. 14 11. 63 378. 94 15. 43
0. 10 0. 087 -0. 084 379. 28 16. 53 380. 88 23. ,76
0. 15 0. 157 -0. 149 379. 78 23. 79 380. 63 34. ,51
0. 20 0. 176 -0. 166 380. 47 32. 74 380. 29 44. ,43
0. 25 0. 139 -0. 132 381. 11 43. 59 379. 98 53. .41
0. 30 0. 056 -0. 055 381. 51 55. 94 379. 78 61. .64
0. .35 -0. 043 0. 042 381. .54 68. 90 379. 76 69. ,58
0. 40 -0. 127 0. 127 381. 19 81. 40 379. 94 77. ,75
0. .45 -0. 170 0. 172 380. 58 92. 53 380. 26 86. .64
0, 50 -0. 159 0. 161 379. 91 101. 83 380. 62 90. .50
0. .55 -0. 098 0. 098 379. 38 109. 41 380. 89 107. ,26
0. .60 -0. 005 0. 004 379. 17 115. 93 381. 00 118. .57
0, ,65 0. 092 -0. 089 379. ,35 122. 41 380. 91 129. ,91
0. .70 0. 160 -0. 151 379. .87 129. .88 380. 66 140. ,75
0. .75 0. 176 -0. ,166 380. .55 139. 08 380. 32 150. ,74
0, ,80 0. 135 -0. .128 381. ,19 150. .17 380. 01 159. ,80
0, .85 0. 051 -0. ,049 381. .57 162. .72 379. ,82 168, ,14
0, .90 -0. ,049 0. .048 381. .57 175. .82 379. ,81 176. .21
0 .95 -0. .131 0. ,131 381. .21 188. .40 380. .01 184, ,54
1 .00 -0. .171 0 .173 380. .59 199. .56 380. .33 193 .62
1 .05 -0. .157 0 .158 379 .92 208 .89 380 ; 68 203 .67
1 .10 -0. .093 0 .093 379 .41 216 .51 380, .95 214 .61
1 .15 0. .001 -0 .002 379 .22 223 .13 381, .05 226 .07
1 .20 0 .097 -0 .093 379 .42 229 .77 380, .95 237 .52
1 .25 0 .162 -0 .154 379 .95 237 .47 380 .68 248 .46
1 .30 0 .175 -0 .165 380 .64 246 .91 380 .35 258 .53
1 .35 0 .131 -0 .125 381 .27 258 .24 380 .04 267 .67
1 .40 0 .045 -0 .044 381 .63 270 .98 379 .86 276 .11
1 .45 -0 .054 0 .053 381 .61 284 .22 379 .87 284 .31
1 .50 -0 .135 0 .135 381 .22 296 .87 380 .07 292 .82
1 .55 -0 .172 0 .174 380 .60 308 .07 380 .40 302 .08
1 .60 -0 . 155 0 .156 379 .93 317 .41 380 .75 312 .32
1 . 65 -0 .088 0 .088 379 .43 325 .09 381 .01 323 .43
1 .70. 0 .007 -0 .008 379 .26 331 .81 381 .10 335 .05
1 .75 0 .102 -0 .098 379 .49 338 .62 380 .98 346 .62
1 .80 0 .165 -0 .156 380 .04 346 .53 380 .71 357 .64
1 .85 0 . 174 -0 .164 , 380 .73 356 .22 380 .37 367 .79
1 .90 0 .127 -0 .121 381 .35 367 .79 380 .07 377 .02
1 .95 0 .039 -0 .038 381 .69 380 .72 379 .90 385 .56
2 .00 -0 .060 0 .059 381 .64 394 .09 379 .92 393 .90
DELTA1-2
-2.65
-3.79
-7.22
-10.71
-11.67
-9.80
-5.70
-0.68
3.64
5.88
5.33
2.15
-2.63
-7.48
-10.84
-11.64
-9.61
-5.41
-0.39
3.85
5.93
5.21
1.90
-2.93
-7.74
-10.97
-11.60
-9.42
-5.12
-0.10
4.04
5.98
5.08
1.65
-3.23
-7.99
-11.09
-11.55
-9.21
-4.83
0.19


74
A.4 Flowchart


75
Flowchart Continuous


APPENDIX B
LISTING OF COMPUTED PHASE-PLANE TRAJECTORY COORDINATES
The main program used to solve the non-linear differential
equations called Simnon. Simnon is a program that can solve dif-
ferential, and difference equations using a fourth order Runge-
Kutta integration algorithm, it can also be used to simulate
dynamical systems that are composed of subsystems. The solution
to the problem can be divided in the following steps.
1. Enter system descriptions
2. Simulate
3. Analyze the results
4. Change parameters and initial condition.
This appendix contains the system descriptions program
list, subroutines (Macros) which perform a whole sequence of
commands, and listing of computed phase-plane trajectory
coordinates.


77
Simnon nn Intcrnctive Simulation Progrnm for
Nonlinear Systems. Version 1.00
Copyright (c) Department of Automatic Control
Lund Institute of Technology, Lund, SWEDEN 1986
All Itights Reserved
--861020011371----
CONTINUOUS SYSTEM IIUSDl
STATE XI X2
DEIi DX1 DX2
TIME T
DXl=x2
DX2=((A-l>SIN(xlMM/180+D))180/PI)
A= IE T D=IF T U=IF T K: .1
pi:3.1416926
END
MACRO IIUSDM1
SYST IIUSDl
STORE xl x2
INIT xl:0.7
I NIT x2:0
SI MU 0 .1
AXES II -200 25 V -8000 8000
SHOW x2(xl)
STORE xl x2
INIT xl:xl
INIT x2:x2
SI MU .1 .2
SHOW x2(xl)
text PHASE PLANE PLOT OF X2 Vs.Xl
mark a 15 1.2
mark "State Variable Xl In Degree
mark n 1.2 15
mark v"Statc VAR X2 DEG/S
END
Print of Simnon store file:STORE For the clearing time t = 0.1 sec.
TIME Xl X2
0. 0.7 0 .
2.E-3 0.632745 -67.2490
4 ,E-3 0.431045 -134.434
6.E-3 9.50969E-2 -201.407
0 E-3 -0.374773 .-260.345
1 .E-2 -0.978113 -334.947
1.2E-2 -1.71435 -401.233
1.4E-2 -2.5028 -467.148
1.6E-2 -3.58266 -532.639
1.0E-2 -4.71304 -597.66
2 .E-2 -5.97296 -662.17
2.2E-2 -7.36135 -726.132
2.4 E-2 -8.8771 -789.510
2.6E-2 -10.519 -852.308
2.8E-2 -12.2859 -914.488
3.E-2 -14.1766 -976.055
3.2E-2 -16.1097 -1037.01
3.4 E-2 -18.3242 -1097.38
3.6E-2 -20.5789 -1157.19
3.8E-2 -22.9526 -1216.46
4 E-2 -25.4444 -1275.27
4.2E-2 -28.0534 -1333.66
4.4E-2 -30.7708 -1391.71


5.8E-2
6.E-2
6.2E-2
6.4E-2
6.6E-2
6.8E-2
7 i'E-2
7.2E-2
7.4E-2
7.6E-2
7.8E-2
8 .E-2
8.2E-2
8.4E-2
8.6E-2
8.8E-2
9.E-2
9.2E-2
9.4E-2
9.6E-2
9.8E-2
9.99999E-2
0.100027
0.100055
0.100082
0.100386
0.102
0.104
0.106
0.108
0.11
0.112
0.114
0.116
0.118
0.12
0.122
0.124
0.126
0.128
0.13
0.132
0.134
0.136
0.138
0.14
0.142
0.144
0.146
0.148
0.15
0.152
0.154
0.156
0.158
0.16
0.162
0.164
0.166
0.168
0.17
0.172
0.174
0.176
0.178
0.18
0.182
0.184
0.186
0.188
0.19
0.192
-53.0963
-56.7523
-60.5289
-64.4284
-68.4534
-72.6072
-76.8933
-81.3159
-85.8799
-90.5904
-95.4534
-100.475
-105.663
-111.024
-116.567
-122.3
-128.233
-134.374
-140.734
-147.322
-154.149
-161.224
-161.323
-161.421
-161.518
-162.582
-167.57
-172.188
-175.082
-176.255
-175.713
-173.454
-169.471
-163.76
-156.33
-147.218
-136.503
-124.331
-110.927
-96.5966
-81.7202
-66.7242
-52.0439
-38.082
-25.1769
-13.585
-3.4791
5.04026
11.9273
17.172
20.7829
22.7741
23.1564
21.9319
19.0936
14.6283
8.52321
0.777713
-8.58175
-19.4831
-31.7898
-45.2861
-59.6713
-74.5703
-89.5583
-104.2
-118.088
-130.879
-142.31
-152.199
-160.437
-166.968
-1798.14
-1857.98
-1918.83
-1980.89
-2044.41
-2109.64
-2176.82
-2246.24
-2318.16
-2392.87
-2470.66
-2551.81
-2636.61
-2725.34
-2818.26
-2915.6
-3017.59
-3124.41
-3236.18
-3352.98
-3474.8
-3601.56
-3590.92
-3579.22
-3567.51
-3437.23
-2741.42
-1877.52
-1016.21
-157.811
700.042
1560.13
2423.38
3286.59
4140.15
4966.03
5736.75
6416.12
6962.64
7335.88
7504.68
7454.62
7191.87
6741.65
6142.01
5435.58
4662.08
3853.62
3032.73
2212.75
1399.39
592.787
-210.352
-1014.76
-1824.66
-2641.81
-3463.42
-4280.01
-5073.58
-5816.5
-6472.26
-6998.82
-7354.92
-7508.26
-7443.07
-7164.14
-6695.42
-6074.
-5341.8
-4538.16
-3694.95
-2834.66


Siranon nn Interactive Simulation Program for
Nonlinear Systems. Version 1.00
Copyright (c) Department of Automatic Control
Lund Institute of Technology, Lund, SWEDEN 1986
All Rights Reserved
---861020011371----
CONTINUOUS SYSTEM IIUSD2
STATE XI X2
DIOR DX1 DX2
TIME T
DXl=x2
DX2=((A-B*SIN(xl*PI/180+D))*180/PI)
A=IF T D=IF t B=IF T Z:0.125
pi :3.1415926
END
MACRO IIUSDM2
SYST IIUSD2
STORE xl x2
INIT xl:0.7
INIT x2:0
SIMU 0 0.18/Z1
AXES H -700 25 V -11000 8000
SHOW x2(xl)/Zl
STORE xl x2
INIT xl:0.7
INIT x2:0
PAR Z:0.08
SI MU 0 0.19/Z2
SHOW x2(xl)/Z2
STORE xl x2
INIT xl:0.7
INIT x2:0
PAR Z:0.09
SI MU 0 0.2/Z3
SHOW x2(xl)/Z3
STORE xl x2
INIT xl:0.7
INIT x2:0
PAR Z:0.1
SI MU 0 0.21/Z4
SHOW x2(xl)/Z4
STORE xl x2
INIT xl:0.7
INIT x2:0
PAR Z:0.115
SI MU 0 0.21/Z5
SHOW x2(xl)/Z5
STORE xl x2
INIT xl:0.7
INIT x2:0
PAR Z:0.125
SI MU 0 0.21/Z6
SHOW x2(xl)/Z6
text PHASE PLANE PLOT OF X2 VS. Xl
mark a 15 1.2
mark "STATE VARIABLE Xl In DEG
mnrk n 1.2 15
mnrk v"STATE VAR X2 DEG/S
END


80
Print
TIME
0.
2 .E-3
4 .E-3
6. E-3
8 .E-3
1. E-2
1.2E-2
1.4E-2
1.6E-2
1.8E-2
2. E-2
2. 2E-2
2.4E-2
2.6E-2
2.8E-2 '
3. E-2
3.2E-2
3.4E-2
3.6F.-2
i. ii i-:
4 .E-2
4.2E-2
4.4E-2
4.6E-2
4.8E-2
5. E-2
5.2E-2
5.4E-2
5.6E-2
5.8E-2
6. E-2
6.2E-2
6.4E-2
6.6E-2
6.8E-2
7 .E-2
7.2E-2
7.4E-2
7.6E-2
7.8E-2
8. E-2
8.2E-2
8.4E-2
8.6E-2
8.8E-2
9. E-2
9.2E-2
9.4E-2
9.6E-2
9.8E-2
9.99999E
0.102
0.104
0.106
0.108
0.11
0.112
0.114
0.116
0.118
0.12
0.122
0.124
0.126
0.128
0.13
0.132
0.134
0.136
0.138
0.14
0.142
0.144
0.146
0.148
0.15
0.152
of Simnon store file:STORE For the clearing time t = 0.07 sec.
XI X2
0.7 0.
0.632745 -67.2498
0.431045 -134.434
9.50 969E-2 -201.487
-0.374773 -268.345
-0.978113 -334.947
-1.71435 -401.233
-2.5828 -467.148
-3.58266 -532.639
-4.71304 -597.66
-5.97296 -662.17
-7.36135 -726.132
-8.8771 -789.518
-10.519 -852.308
-12.2859 -914.488
-14.1766 -976.055
-16.1897 -1037.01
-18.3242 -1097.38
-20.5789 -1157.19
22.9526 -1216.46
-25.4444 -1275.27
-28.0534 -1333.66
-30.7788 -1391.71
-33.6201 -1449.53
-36.5769 -1507.2
-39.6489 -1564.85
-42.8363 -1622.62
-46.1395 -1680.64
-49.5592 -1739.09
-53.0963 -1798.14
-56.7523 -1857.98
-60.5289 -1918.83
-64.4284 -1980.89
-68.4534 -2044.41
-72.6072 -2109.64
-76.8933 -2176.82
-81.2454 -2160.13
-85.4915 -2075.63
-89.5074 -1930.59
-93.1766 -1729.81
-96.3938 -1479.76
-99.068 -1188.22
-101.125 -863.856
-102.508 -515.958
-103.179 -154.14
-103.122 211.815
-102.336 572.082
-100.844 916.896
-98.685 1236.68
-95.9196 1522.23
-92.6246 1764.98
-88.8933 1957.35
-84.833 2093.1
-80.5617 2167.73
-76.2046 2178.76
-71.8894 2125.94
-67.7422 2011.2
-63.8832 1838.58
-60.4226 1613.85
-57.4578 1344.2
-55.0703 1037.82
-53.3251 703.567
-52.2685 350.677
-51.9285 -11.4194
-52.314 -373.219
-53.4148 -725.215
-55.202 -1057.99
-57.6278 -1362.32
-60.6264 -1629.38
-64.1151 -1851.03
-67.9955 -2020.17
-72.1569 -2131.11
-76.4786 -2179.96
-80.8341 -2164.89
-85.0958 -2086.33
-89.1388 -1946.87
-92.8457 -1751.13


81
" Print of Simnon store file:STORE For the clearing time t
TIME XI X2
0. 0.7 0.
2.E-3 0.632745 -67.2498
4 .E-3 0.431045 -134.434
6.E-3 9.50969E-2 -201.487
8.E-3 -0.374773 -268.345
1 .E-2 -0.978113 -334.947
1.2E-2 -1.71435 -401.233
1.4E-2 -2.5828 -467.148
1.6E-2 -3.58266 -532.639
1 8E-2 -4.71304 -597.66
2 .E-2 -5.97296 -662.17
2.2E-2 -7.36135 -726.132
2.4E-2 -8.8771 -789.518
2.6E-2 -10.519 -852.308
2.8E-2 -12.2859 -914.488
3.E-2 -14.1766 -976.055
3.2E-2 -16.1897 -1037.01
3.4E-2 -18.3242 -1097.38
3.6E-2 -20.5789 -1157.19
3.8E-2 -22.9526 -1216.46
4.E-2 -25.4444 -1275.27
4.2E-2 -28.0534 -1333.66
4.4E-2 -30.7788 -1391.71
4.6E-2 -33.6201 -1449.53
4.8E-2 -36.5769 -1507.2
5.E-2 -39.6489 -1564.85
5.2E-2 -42.8363 -1622.62
5.4E-2 -46.1395 -1680.64
5.6E-2 -49.5592 -1739.09
5.8E-2 -53.0963 -1798.14
6.E-2 -56.7523 -1857.98
6.2E-2 -60.5289 -1918.83
6.4E-2 -64.4284 -1980.89
6.6E-2 -68.4534 -2044.41
6.8E-2 -72.6072 -2109.64
7.E-2 -76.8933 -2176.B2
7.2E-2 -81.3159 -2246.24
7.4E-2 -85.8799 -2318.16
7.6E-2 -90.5904 -2392.87
7.8E-2 -95.4534 -2470.66
8 .E-2 -100.475 -2551.B1
8.2E-2 -105.291 -2221.58
8.39476E-2 . -109.228 -1812.4
8.4E-2 -109.322 -1800.73
8.43141E-2 -109.877 -1730.12
B.6E-2 -112.464 -1334.1
8.8E-2 -114.637 -834.679
9.E-2 -115.789 -315.251
9.2E-2 -115.893 211.624
9.4E-2 -114.946 733.507
9.6E-2 -112.971 1237.88
9.8E-2 -110.015 1711.99
9.99999E-2 -106.152 2142.89
0.102 -101.481 2517.68
0.104 -96.1265 2824.08
0.106 -90.2372 3051.22
0.100 -83.9803 3190.52
0.11 -77.5374 3236.58
0.112 -71.0973 3187.75
0.114 -64.8482 3046.31
0.116 -58.9699 2818.22
0.118 -53.6272 2512.4
0.12 -48.9648 2139.91
0.122 -45.1039 1712.96
0.124 -42.1409 1244.19
0.126 -40.1467 746.109
0.128 -39.1678 230.917
0.13 -39.2265 -289.499
0.132 -40.3215 -803.36
0.134. -42.4277 -1298.75
0.136 -45.496 -1763.44
0.138 -49.4526 -2184.88
0.14 -54.1983 -2550.45
0.142 -59.6091 -2848.02
0.144 -65.5378 -3066.74
0.146 -71.8177 -3197.97
0.140 -78.2676 -3236.14
0.08 sec.


82
" Print of Simnon store fil
TIME XI
0 . 0.7
2 E-3 0.632745
4.E-3 0.431045
6.E-3 9.50969E-2
8.E-3 -0.374773
l.E-2 -0.978113
1.2E-2 -1.71435
1 4Er-2 -2.5828
1.6E-2 -3.58266
1.8E-2 -4.71304
2 E-2 -5.97296
2.2E-2 -7.36135
2.4E-2 -8.8771
2.6E-2 -10.519
2.0F.-2 -12.2059
3.E-2 -14.1766
3.2E-2 -16.1897
3.4E-2 -18.3242
3.6E-2 -20.5789
3.8E-2 -22.9526
4.E-2 -25.4444
4.2E-2 -28.0534
4.4E-2 -30.7788
4.6E-2 -33.6201
4.8E-2 -36.5769
5.E-2 -39.6489
5.2E-2 -42.8363
5.4E-2 -46.1395
5.6E-2 -49.5592
5.8E-2 -53.0963
6. E-2 -56.7523
6.2E-2 -60.5289
6.4E-2 -64.4284
6.6E-2 -68.4534
6.8E-2 -72.6072
7.E-2 -76.8933
7.2E-2 -81.3159
7.4E-2 -85.8799
7.6E-2 -90.5904
7.8E-2 -95.4534
8 .E-2 -100.475
8.2E-2 -105.663
8.4E-2 -111.024
8.6E-2 -116.567
8.8E-2 -122.3
9.E-2 -128.233
9.00277E-2 -128.316
9.00555E-2 -128.4
9.00833E-2 -128.483
9.03S16E-2 -129.361
9.2E-2 -133.593
9.4E-2 -137.54
9.6E-2 -140.01
9.8E-2 -140.966
9.99999E-2 -140.394
0.102 -138.302
0.104 -134.721
0.106 -129.708
0.108 -123.353
0.11 -115.781
0.112 -107.163
0.114 -97.712
0.116 -87.6824
0.118 -77.3594
6.12 -67.0444
0.122 -57.0369
0.124 -47.6169
0.126 -39.031
0.128 -31.484
0.13 -25.1363
0.132 -20.1065
0.134 -16.4769
0.136 -14.3002
0.138 -13.6051
0.14 -14.4004
0.142 -16.6761
:STORE For the clearing time t =
X2
0.
-67.2498
-134.434
-201.487
-268.345
-334.947
-401.233
-467.148
-532.639
-597.66
-662.17
-726.132
-789.518
-852.308
-914.488
-976.055
-1037.01
-1097.38
-1157.19
-1216.46
-1275.27
-1333.66
-1391.71
-1449.53
-1507.2
-1564.85
-1622.62
-1680.64
-1739.09
-1798.14
-1857.98
-1918.83
-1980.89
-2044.4 1
-2109.64
-2176.82
-2246.24
-2318.16
-2392.87
-2470.66
-2551.81
-2636.61
-2725.34
-2818.26
-2915.6
-3017.59
-3009.26
-3000.08
-2990.9
-2891.62
-2333.53
-1608.52
-858.528
-96.1357
667.604
1421.81
2154.42
2850.98
3494.09
4063.41
4536.88
4892.84
5112.87
5184.84
5105.07
4879.03
4520.37
4048.42
3485.23
2852.65
2170.37
1454.87
719.449
-25.1134
-769.364
-1503.81
0.09 sec.


83
Print
TIME
0.
2 E-3
.4 .E-3
6. E-3
0 .E-3
l.E-2
1.2E-2
1. AE-2
1.6E-2
1.8E-2
2. E-2
2.2E-2
2 4E-2
2.6E-2
2.8E-2
3. E-2
3.2E-2
3.AE-2
3.6E-2
3.BE-2
A. E-2
A.2E-2
A.AE-2
4.6E-2
4.8E-2
5. E-2
5.2E-2
"5.4E-2
5.6E-2
5.8E-2
6. E-2
6.2E-2
6.4E-2
6.6E-2
6.8E-2
7. E-2
7.2E-2
7.4E-2
7.6E-2
7.8E-2
8. E-2
8.2E-2
8.4E-2
8.6E-2
8.8E-2
9. E-2
9.2E-2
9.4E-2
9.6E-2
9.8E-2
9.99999E-
0:100027
0.100055
0.100082
0.100386
0.102
0.104
0.106
0.108
0.11
0.112
0.114
0.116
0.118
0.12
0.122
0.124
0.126
0.128
0.13
0.132
0.134
0.136
0.138
of Simnon store file: STOKE For the clearing time t 0.1 r.oc
XI X2
0.7 0.
0.632745 -67.2498
0.431045 -134.434
9.50969E-2 -201.487
-0.374773 -.268.345
-0.970113 -334.947
-1.71435 -401.233
-2.5820 -467.140
-3.50266 -532.639
-4.71304 -597.66
-5.97296 -662.17
-7.36135 -726.132
-8 8771 -709.518
-10.519 -052.300
-12.2059 -914.480
-14.1766 -976.055
-16.1097 -1037.01
-18.3242 -1097.38
-20.5709 -1157.19
-22.9526 -1216.46
-25.4444 -1275.27
-28.0534 -1333.66
-30.7788 -1391.71
-33.6201 -1449.53
-36.5769 -1507.2
-39.6489 -1564.05
-42.8363 -1622.62
-46.1395 -1600.64
-49.5592 -1739.09
-53.0963 -1798.14
-56.7523 -1857.98
-60.5209 -1910.03
-64.4204 -1980.09
-68.4534 -2044.41
-72.6072 -2109.64
-76.8933 -2176.82
-01.3159 -2246.24
-05.0799 -2310.16
-90.5904 -2392.. 07
-95.4534 -2470.66
-100.475 -2551.81
10 5.663 -2636.61
-111.024 -2725.34
-116.567 -201B.26
-122.3 -2915.6
-128.233 -3017.59
-134.374 -3124.41
-140.734 -3236.10
-147.322 -3352.98
-154.149 -3474.8
-161.224 -3601.56
-161.323 -3590.92
-161.421 -3579.22
-161.518 -3567.51
-162.582 -3437.23
-167.57 -2741.42
-172.188 -1877.52
-175.082 -1016.21
-176.255 -157.811
-175.713 700.042
-173.454 1560.13
-169.471 2423.38
-163.76 3286.59
-156.33 4140.15
-147.218 4966.03
-136.503 5736.75
-124.331 6416.12
-110.927 6962.64
-96.5966 7335.88
-81.7202 7504.68
-66.7242 7454.62
-52.0439 7191.87
-38.082 6741.65
-25.1769 6142.01


84
" Print of Simnon store file:STORE For the clearing time t == 0.115 sec.
TIME XI X2
0. 0.7 0.
2.5E-3 0.594918 -84.0545
5.E-3 0.279834 -167.981
7.5E-3 -0.244773 -251.652
l.E-2 -0.978112 -334.947
1.25E-2 -1.9191 -417.749
1.5E-2 -3.06635 -499.949
1.75E-2 -4.41826 -581.451
2.E-2 -5.97296 -662.17
2.25E-2 -7.72839 -742.033
2.5E-2 -9.68236 -820.989
2.75E-2 -11.8325 . -899.
3.E-2 -14.1766 -976.055
3.25E-2 -16.712 -1052.16
3.5E-2 -19.4366 -1127.35
3.75E-2 -22.3481 -1201.69
4 E-2 -25.4444 -1275.27
4.25E-2 -28.7239 -1348.2
4.5E-2 -32.185 -1420.65
4.75E-2 -35.8268 -1492.79
5. E-2 -39.6489 -1564.85
5.25E-2 -43.6513 -1637.09
5.5E-2 -47.8347 -1709.8
5.75E-2 -52.2009 -1783.31
6. E-2 -56.7523 -1857.98
6.25E-2 -61.4922 -1934.22
6.5E-2 -66.425 -2012.46
6.75E-2 -71.5565 -2093.16
7.E-2 -76.8933 -2176.82
7.25E-2 -82.4435 -2263.97
7.5E-2 -88.2165 -2355.15
7.75E-2 -94.223 -2450.91
8 .E-2 -100.475 -2551.81
8.25E-2 -106.987 -2658.42
8.5E-2 -113.773 -2771.26
8.75E-2 -120.849 -2890.84
9.E-2 -128.233 -3017.6
9.25E-2 -135.943 -3151.89
9.5E-2 -143.999 -3293.95
9.75E-2 -152.419 -3443.88
0.1 -161.224 -3601.57
0.1025 -170.433 -3766.69
0.105 -180.064 -3938.63
0.1075 -190.131 -4116.48
0.11 -200.65 -4298.. 97
0.1125 -211.629 -4484.47
0.115 -223.073 -4671.02
0.117371 -233.631 -4199.1
0.1175 -234.17 -4174.64
0.118788 -239.402 -3955.41
0.12 -244.091 -3790.73
0.1225 -253.264 -3575.3
0.125 -262.104 -3523.71
0.1275 -271.018 -3634.23
0.13 -280.412 -3908.45
0.1325 -290.699 -4350.02
0.135 -302.303 -4961.34
0.1375 -315.644 -5737.5
0.14 -331.111 -6656.5
0.1425 -349.003 -7666.26
0.145 -369.44 -8672.53
0.1475 -392.249 -9537.98
0.15 -416.881 -1.01058E4
0.1525 -442.423 -1.0251E4
0.155 -467.749 -9935.05
0.1575 -491.768 -9224.51
0.16 -513.658 -8257.85
0.1625 -532.971 -7187.41
0.165 -549.611 -6135.47
0.1675 -563.729 -5179.17
0.17 -575.619 -4356.14
0.1725 -585.631 -3677.37
0.175 -594.123 -3139.24


85
" Print of Simnon store file:STORE For the clearing time t 0.125 sec.
TIME XI X2
0 . 0.7 0.
2.5E-3 0.594918 -84.0545
5.E-3 0.279834 -167.981
1.5E-3 -0.244773 -251.652
1 .E-2 -0.978112 -334.947
1.25E-2 -1.9191 -417.749
1.5E-2 -3.06635 -499.949
1.75E-2 -4.41826 -581.451
2.E-2 -5.97296 -662.17
2.25E-2 -7.72839 ' -742.033
2.5E-2- -9.68236 -820.989
2.75E-2 -11.8325 -899.
3.E-2 -14.1766 -976.055
3.25E-2 -16.712 -1052.16
3.5E-2 -19.4366 -1127.35
3.75E-2 -22.3481 -1201.69
4. E-2 -25.4444 -1275.27
4.25E-2 -28.7239 -1348.2
4.5E-2 -32.185 -1420.65
4.75E-2 -35.0268 -1492.79
5.E-2 -39.6489 . -1564.85
5.25E-2 -43.6513 -1637.09
5.5E-2 -47.8347 -1709.8
5.75E-2 -52.2009 -1783.31
6.E-2 -56.7523 -1857.98
6.25E-2 -61.4922 -1934.22
6.5E-2 -66.425 -2012.46
6.75E-2 -71.5565 -2093.16
7 .E-2 -76.8933 -2176.82
7.25E-2 -82.4435 -2263.97
7.5E-2 -88.2165 -2355.15
7.75E-2 -'94.223 -2450.91
8 .E-2 -100.475 -2551.81
8.25E-2 -106.987 -2658.42
8.5E-2 -113.773 -2771.26
8.75E-2 -120.849 -2890.84
9. E-2 -128.233 -3017.6
9.25E-2 -135.943 -3151.89
9.5E-2 -143.999 -3293.95
9.75E-2 -152.419 -3443.88
0.1 -161.224 -3601.57
0.1025 -170.433 -3766.69
0.105 -180.064 -3938.63
0.1075 -190.131 -4116.48
0.11 -200.65 -4298.97
0.1125 -211.629 -4484.47
0.115 -223.073 -4671.02
0.1175 -234.983 -4856.33
0.12 -247.352 -5037.9
0.1225 -260.167 -5213.14
0.125 -273.41 -5379.5
0.1275 -287.237 -5727.32
0.13 -302.239 -6311.97
0.1325 -318.972 -7106.05
0.135 -337.904 -8060.29
0.1375 -359.331 -9081.69
0.14 -383.249 -1.00253E4
0.1425 -409.242 -1.07111E4
0.145 -436.452 -1.09774E4
0.1475 -463.715 -1.07507E4
0.15 -489.836 -1.00833E4
0.1525 -513.888 -9124.79
0.155 -535.367 -8053.48
0.1575 -554.187 -7018.17
0.16 -570.57 -6115.91
0.1625 -584.92 -5397.5
0.165 -597.728 -4883.26
0.1675 -609.511 -4577.78
0.17 -620.79 -4479.89
0.1725 -632.082 -4588.06
0.175 -643.901 -4901.76
0.1775 -656.761 -5419.53
0.18 -671.163 -6133.15
6.1825 -687.57 -7017.47


Simnon an Interactive Simulation Program for
Nonlinear Systems. Version 1.00
Copyright (c) Department of Automatic Control
Lund Institute of Technology, Lund, SWEDEN 1986
All Rights Reserved
---861020011371----
CONTINUOUS SYSTEM IICP
STATE XI X2
DER DX1 DX2
TIME T
DXl=x2
DX2=((P-B*SIN(xl*PI/180+D))*180/PI)
D: 1.38
B:7358
P:-151.2
PI:3.1415926
END
MACRO IICPM
SYST IICP
STORE xl x2
INIT xl:7
INIT x2:0
SIMU 0 ,2/Cl
AXES H 250 75 V -11000 10000
SHOW x2(xl)/Cl
STORE xl x2
INIT xl:-190
INIT x2:0
PAR P:111.69
SIMU 0 .210,2
SHOW x2(xl)/C2
STROE xl x2
INIT xl:50
INIT x2:0
PAR P:-19.7
SIMU 0 .25/C3
SHOW x2(xl)/C3
STROE xl x2
INIT xl:70
INIT x2:0
PAR P:151.2
SIMU 0 .25/C4
SHOW x2(xl)/C4
text PHASE PLANE PLOT OF X2 Vs.Xl
mark a 15 1.2
mark "State Variable Xl In Degree
mark a 1.2 15
mark v"State VAR X2 DEG/S
END


87
Foe Pml2
TIME
0.1
0.102
0.104
0.106
0.108
0.11
0; 112
0.114
0.116
0.118
0.12
0.122
0.124
0.126
0.128
0.13
0.132
0.134
0.136
0.138
0.14
0.142
0.144
0.146
0.148
0.15
0.152
0.154
0.156
0.158
0.16
0.162
0.164
0.166
0.168
0.17
0.172
0.174
0.176
0.178
0.18
0.182
0.184
0.186
0.188
0.19
0.192
0.194
0.196
0.198
0.2
0.202
0.204
0.206
0.208
0.21
0.212
0.214
0.216
0.218
0.22
0.222
0.224
0.226
0.228
0.23
0.232
0.234
0.236
0.238
0.24
0.242
-0.16 p.u
XI X2
-190. 0.
-189.199 801.708
-186.788 1611.39
-182.744 2435.18
-177.036 3275.42
-169.634 4128.3
-160.522 4981.47
-149.723 5811.86
-137.313 6584.69
-123.452 7255.08
-108.394 7773.16
-92.4921 8093.01
-76.1762 8183.15
-59.919 8034.96
-44.1845 7665.15
-29.3813 7111.22
-15.8299 6421.92
-3.7506 5646.89
6.72928 4828.78
15.5572 3999.26
22.7323 3178.38
28.2833 2375.97
32.25 1593.95
34.6703 828.703
35.5713 73.1355
34.9635 -681.547
32.8397 -1444.28
29.1758 -2222.63
23.9357 -3021.03
17.0786 -3838.81
8.57264 -4667.88
-1.58876 -5490.31
-13.3649 -6276.73
-26.6448 ' -6986.42
-41.227 -7570.6
-56.8102 -7979.59
-73.0012 -8172.86
-89.3432 -8129.1
-105.362 -7851.99
-120.614 -7369.09
-134.73 -6724.35
-147.437 -5967.63
-158.557 .-5145.33
-167.999 -4294.41
-175.733 -3440.2
-181.767 -2596.97
-186.132 -1770.01
-188.858 -957.981
-189.97 -155.318
-189.48 645.743
-187.383 1453.31
-183.658 2274.04
-178.275 3111.16
-171.204 3962.26
-162.424 4816.83
-151.948 5653.99
-139.841 6441.24
-126.245 7135.35
-111.394 7686.85
-95.6241 8048.22
-79.352 8184.4
-63.0461 8081.9
-47.176 7752.68
-32.1641 7230.86
-18.3494 6563.85
-5.97183 5801.81
4.82447 4989.22
13.9743 4159.99
21.4685 3336.38
27.3315 2530.02
31.6024 1744.19
34.3202 976.16
For Pml2 = -0.12 p.u
TIME XI
0 . 1 7 .
0 . 102 6.20 184
0 .104 3.8093
0 . 106 -0.1 69865
0. .108 -5.7 1633
0, .11 -12. 79
0, .112 -21. 3178
0 .114 -31 . 1013
0 .116 -42. 2046
0 .118 -54 . 1486
0 . 12 -66. 7133
0 . 122 -79. 5509
0, .124 -92. 2888
0. .126 -104 . 558
0 .128 -116 . 02
0, .13 -126 .380
0 .132 -135 .435
0 .134 -142 .998
0 .136 -148 . 963
0 .138 -153 . 261
0 . 1 4 -155 . 053
0 .142 -156 .718
0 .144 -155 . 851
0 .146 -153 .258
0 .148 -148 .958
0 .15 -142 .991
0 .152 -135 .428
0 .154 -126 .379
0 . 156 -116 .01
0 .158 -104 .547
0 .16 -92. 2774
0 .162 -79. 5392
0 .164 -66. 7016
0 .166 -54. 1373
0 .168 -42. 194
0 .17 -31. 1716
0 .172 -21. 3093
0 . 174 -12. 7827
0 .176 -5.7 1043
0 .178 -0.1 65393
0 .18 3.81 23
0 .182 6.20334
0 .184 6.99999
0 .186 6.20031
0 .108 3.80624
0 .19 -0.1 74454
0 .192 -5.72244
0 .194 -12 . 7975
0, .196 -21. 3268
0 , .198 -31. 1916
0 , .2 -42 . 2159
0. .202 -54 . 1607
0, .204 -66. 7259
0. .206 -79. 5636
0. .208 -92. 3013
0 .21 -104 .57
0, .212 -116 .031
0, .214 -126 .398
0 , .216 -135 .444
0. .218 -143 .004
0. ,22 -148 . 968
0. ,222 -153 .264
0. .224 -155 . 854
0. ,226 -156 .718
0. .228 -155 .849
0. 23 -153 .255
0. 232 -148 . 953
0. 234 -142 . 985
0. 236 -135 .42
0. 238 -126 .369
0. 24 -115 . 999
0. 242 -104 .535
X2
0 .
-798.015
-1593.93
-2383.71
-3159.52
-3908.15
-4609.96
-5238.95
-5764.42
-6154.47
-6381.03
-6425.28
-6281.67
-5959.13
-5478.93
-4870.04
-4163.78
-3389.19
-2570.25
-1724.82
-865.217
0.390772
865.996
1725.59
2571 .
3389.92
4164.45
4870.64
5479.42
5959.5
6281.88
6425.32
6380.91
6154.18
5764.
5238.41
4609.34
3907.47
3158.81
2382.97
1593.18
797.259
-0.75876
-798.775
-1594.69
-2384.46
-3160.26
-3908.86
-4610.61
-5239.52
-5764.88
-6154.78
-6381.17
-6425.23
-6281.43
-5958.72
-5478.38
-4869.38
-4163.04
-3388.4
-2569.42
-1723.97
-864.356
1.25462
866.857
1726.44
2571.83
3390.71
4165.19
4871.3
5479.97
5959.9


88
For Pral2 =
TIME
0.1
0.102
0.104
0.106
0.108
0.11
0.112
0.114
0.116
0.118
0.12
0.122
0.124
0.126
0.128
0.13
0.132
0.134
0.136
0.338
0.14
0.142
0.144
0.146
0.148
0.15
0.152
0.154
0.156
0.158
0.16
0.162
0.164
0.166
0.168
0.17
0.172
0.174
0.176
0.178
0.18
0.182
0.184
0.186
0.108
0.19
0.192
0.194
0.196
0.198
0.2
0.202
0.204
0.206
0.208
0.21
0.212
0.214
0.216
0.218
0.22
0.222
0.224
0.226
0.228
0.23
0.232
0.234
0.236
0.238
0.24
0.242
.18 p.u For Pml2 = -0.2 p.u
XI X2 TIME XI X2
50. 0. 0.1 70. 0.
49.3421 -658.912 0.103 68.9811 -682.49
47.3564 -1329.8 0.106 65.8668 -1403.28
44.0075 -2023.77 0.109 60.4855 -2199.91
39.2397 -2749.96 0.112 52.5569 -3106.68
32.9822 -3514.09 0.115 41.7098 -4148.32
25.1578 -4316.38 0.118 27.5286 -5326.92
15.6963 -5148.96 0.121 9.6511 -6600.34
4.55403 -5993. 0.124 -12.0594 -7056.86
-8.26209 -6816.24 0.127 -37.2844 -8906.23
-22.6665 -7572.59 0.13 -65.0603 -9521.16
-38.4701 -8205.36 0.133 -93.8087 -9539.25
-55.366 -8655.61 0.136 -121.694 -8961.39
-72.9377 -8874.57 0.139 -147.144 -7952.33
-90.6923 -8836.49 0.142 -169.209 -6742.78
-108.115 -8546.2 0.145 -187.597 -5527.22
-124.731 -8037.67 0.148 -202.485 -4422.08
-140.156 -7364.29 0.151 -214.289 -3474.67
-154.118 -6585.62 . 0.154 -223.496 -2688.91
-166.464 -5755.89 0.157 -230.567 -2047.43
-177.136 -4917.26 0.16 -235.099 -1525.05
-186.147 -4097.86 0.163 -239.81 -1095.7
-193.551 -3313.03 0.166 -242.542 -735.061
-199.425 -2568.02 0.169 -244.268 -421.268
-203.848 -1860.B2 0.172 -245.097 -13 4.561
-206.889 -1184.61 0.175 -245.084 143.461
-208.601 -529.656 0.170 -244.227 430.743
-209.015 115.431 0.181 -242.472 745.715
-208.138 762.687 0.184 -239.705 1108.19
-205.955 1424.01 0.187 -235.753 1540.12
-202.425 2110.33 0.19 -230.371 2065.88
-197.491 2830.51 0.193 -223.238 2711.59
-191.077 3589.88 0.196 -213.956 3502,29
-183.105 4388.16 0.199 -202.061 4454.95
-173.504 5216.84 0.202 -187.068 5564.64
-162.23 6056.24 0.205 -168.563 6782.23
-149.293 6873.23 0.208 -146.382 7988.77
-134.783 7620.87 0.211 -120.837 8987.62
-118.894 8241.96 0.214 -92.897 9548.36
-101.939 8677.66 0.217 -64.1512 9510.28
-84.3395- 8880.08 0.22 -36.4349 8878.35
-66.5911 8824.91 0.223 -11.3105 7818.86
-49.2079 8518.64 0.226 10.2799 6559.25
-32.6608 7996.53 0.229 28.0358 5287.55
-17.3297 7312.56 0.232 42.1047 4112.96
-3.47866 6526.14 0.235 52.8524 3075.81
8.74262 5690.83 0.238 60.6946 2173.01
19.. 2 8 4847.93 0.241 66. 1379.35
28.1477 4024.67 0.244 69.0455 660.371
35.4017 3235.64 0.247 69.9993 -21.5226
41.1164 2485.47 0.25 68.9152 -704.638
45.368 1771.7 0.253 65.732 -1427.27
4 8.2227 1087.18 0.256 60.2745 -2226.89
49.7296 421.927 0.259 52.2592 -3137.64
49.9156 -235.59 0.262 41.3126 -4183.77
48.7843 -897.445 0.265 27.0108 -5366.34
46.3151 -1575.41 0.268 9.01984 -6641.36
42.4649 -2280.02 0.271 -12.8103 -7894.6
37.1717 -3019.4 0.274 -38.1348 -8933.62
30.3613 -3797.54 0.277 -65.9684 -9531.41
21.957 -4612.11 0.28 -94.7175 -9529.55
11.8959 -5451.63 0.283 -122.547 -8934.8
0.149153 -6292.66 0.286 -147.9 -7915.8
-13.2512 -7097.89 0.289 -169.85 -6703.45
-28.1851 -7816.8 0.292 -188.122 -5490.04
-44 4218 -8390.68 0.295 -202.905 -4389.48
-61.613 -8762.66 0.298 -214.619 -3447.31
-79.3094 -8890.91 0.301 -223.751 2666.47
-97.0033 -8760.2 0.304 -230.761 -2029.17
-114.188 -8386.47 0.307 -236.044 -1510.15
-130.415 -7812.09 0.31 -239.914 -1083.35
-145.34 -7094.19 0.313 -242.612 -724.523


APPENDIX C

THE ENERGY METRIC ALGORITHM
This topic is from the article "The generation of Liapunov
function in control theory by an energy metric algorithm" by,
Edward T. Wall, Professor. Department of Electrical Engineering,
University of Colorado at Denver.
The basic procedure involved consists of the following six
steps:
Step 1:
Describe the system as a set of first-order nonlinear
differential equations.
X.= F.(X), 1 < i < n (1)
Where X is the vector (X,, X ....., X ).
12 n
Step 2:
Create a set of differential equations of integral curves
by forming quotients of first-order nonlinear differential
equations and eliminating dt in each case.
dX. F.(X)
----- = ------ (j>i), and (j^i) (2)
dXj F.(X)
There will be n^n ^ such equations.
2


90
Step 3:
The differential equations of integral curves are converted to
the same number of differential one-forms by multiplying and
clearing the denominator terms.
F.(X) dXj F.(X) dXj = 0, (j>i) (3)
Step 4:
By addition and substitution these differential one-forms
are reduced to a single one-form.
HI = Mi,(X)dX, + HI (X)dX + ... + HI (X)dX (4)
ll 22 n n
where the His are determined by the additions and substitutions.
Step 5:
Perform a line integration of the one-forms. The path
of integration is chosen along the state coordinates to form and
"energy-like" scalar function directly from the system equation.
HI =

Hl1(r1> 0...0)dri+
ou
W V 0--->dV
Vw
..,x ,,t )dr
n-l n n
0
The result of the line integration is the V function.
Step 6:
(5)
The total derivative of V with respect to t is taken to
obtain V. The proper theorems are then considered to determine
stability or a region of stability.


Full Text

PAGE 1

AN INVESTIGATION OF THE POWER SYSTEM STABILITY DOMAINS FOR A TWO-FINITE-MACHINE SYSTEM by Hussein Abdullatif Dghayes University of Colorado, 1985 A thesis submitted to the Faculty of the Graduate School of the University of Colorado in partial fulfillment of the requirements for the degree of Master of Science Department of Electrical Engineering and Computer Science 1988

PAGE 2

This thesis for the Master of Science degree by Hussein Abdullatif Dghayes has been approved for the Department of Electrical Engineering and Computer Science by R. Roemish

PAGE 3

iii Dghayes Abdul latif Hussein. (M.S., Electrical Engineering) An Investigation of the Power System Stability Domains for System Thesis directed by professor Edward T. Wall Abstract-Analysis of the stability of power system following a transient disturbance involves the study of a large set of non-linear differential equations. This thesis presents a systematic procedure for investigation of stability domains of a two-finite machine, interconnected by passive loads, power system model in which the classical representation is used for each machine and resistance is neglected. The direct method of Liapunov is used to determine the stability boundary of the power system. The method is first compared with the equal-area criterion and the phase plane technique, with all three approaches giving identical results for the simplified model With constant damping introduced into the system model, a larger, more accurate, region of stability is obtained by use of Mansour's modification of the Energy-Metric-Algorithm.

PAGE 4

iv ACKNOWLEDGEMENTS This thesis is submitted to the University of Colorado as the final requirement for the degree of Master of Science in Electrical Engineering studies at the Denver campus. This study would not have been possible without the cooperation and support of many people to whom I will always be grateful. I especially wish to express my sincere appreciation to Dr. Edward T Wall, for his cheerful support and encouragement during the preparation of this thesis. Also, the enthusiasm of Dr. William R. Roemish for his suggestions and encouragement throughout this study. I would further like to thank the faculty of the University of Colorado at Denver for their efforts in providing an excellent education. Finally, I would like to thank my family, for their moral support and inspiration.

PAGE 5

OONTENTS CHAPTER I INTRODUCTION. . . . . . . . . . 1 I I. STABILITY.......................................... 4 2.1 Introduction.................................. 4 2.2 Steady-State Stability......... ............... 4 2.3 Dynamic Stabi li ty ............................. 5 2.4 Transient Stability........................... 6 2.5 Methods of Simulation 7 III. CLASSICAL MODEL 3.1 Introduction 9 3.2 Classical Stability Study Assumption... ... 9 3.3 Power System Representation ...... ............. 10 3.4 System Equation of Motion....... .............. 10 3.5 Pre-Fault Power System Model...... ..... ....... 14 3.6 During Fault Power System Model. ..... ..... .... 18 3.7 Post-Fault Power System Model................. 20 IV. NUMERICAL SOLUTION OF THE SWING EQUATION 4.1 Introduction.................................. 23 4.2 Point-by-Point Solution of the Swing Equation ...................................... 23 4.3 Numerical" Methods for Solution of Differential Equation ...................................... 27

PAGE 6

vi V. STABILITY DOMAIN 5.1 Introduction.................................. 30 5 2 Power-Angle Characteristics...... .......... ... 30 VI APPLICATION OF THE PHASE PLANE TRAJECTORIES AND THE STABILITY BOUNDARY 6.1 Introduction.... . . . . . . . 35 6.2 Phase Plane Traj ectory ........................ 36 6.3 Stability Boundary....... ..................... 42 VI I ASYMPT9TIC STABILITY 7.1 Introduction.................................. 46 7.2 Liapunov Stabi I i ty ............................ 46 7.3 Energy Metric Algorithm(EMA) ........... ....... 47 7.4 A Suitable Liapunov Function and an Estimation of )( ............................ 53 VIII. MANSOUR'S MODIFICATION OF THE EMA 8.1 Introduction.................................. 57 8.2 Power System Representation.. ...... ........... 57 8.3 System Equation of Motion ..... ........... ..... 58 8.4 The Unmodified EMA ............................ 59 8.5 Mansour's Modified EMA ... .......... .... ....... 60 IX. CONCLUSIONS. . . . . . . . . . 64

PAGE 7

vii BIBLIOGRAPHY ............................................... 66 APPENDIX A. LISTING OF COMPUTER PROGRAM FOR THE SOLUTION OF THE SWING EQUATION. . . . . . . . . 68 B. LISTING OF COMPUTED PHASE-PLANE TRAJECTORY COORDINATES ............................................ 76 C. THE ENERGY METRIC ALGORITHM ............ ...... ...... .... 89 D. NOTATION............................................... 91

PAGE 8

Table 3.1 TABLES Generator Data On p.u to Generator MVA Base} 3 2 Transmission Line Data On 100 MVA Base viii 11 (Resistance Neglected) ............................ 11 3.3 Generator Data (in p.U to 100 MVA Base) .. .... ........ 12 3.4 Load-Flow Data...................................... 12 3.5 Load-Flow Data (on 100 MVA Base) ...... ............. 12

PAGE 9

FIGURES Figure 3.1 Single-Line Diagram Of Two-Machine Power System.. 11 3.2 Pre-Fault System Reactance Diagram ............... 15 3.3 Pre-Fault Admittance Diagram ... .... ...... ........ 15 3.4 Pre-Fault Power-Angle Curve .... ...... ........... 17 3.5 During Fault System Reactance Diagram ............ 19 3.6 During Fault Admittance Diagram ....... ........... 19 3.7 Post-Fault System Reactance Diagram .............. 21 3.8 Post-Fault Admittance Diagram ........ ............ 21 4.1 Actual and Assumed Values of Pa, W, and 8 as Function of Time .............................. 26 4.2 Plot of 81 and 82 vs. Time ...... ;. ...... ......... 28 4.3 Plot of Delta Differences (81-82) vs. Time. ...... 29 5.1 Classical Power-Angle Curve .......... .... ........ 31 6.1 Phase Plane Plot Of X2 Vs. Xl ..... .............. 37 6.2 Power-Angle Curves And Phase Plane Trajectories ..................................... 40 6.3 Plot of Stable and Unstable Phase Plane Traj ectories ..................................... 45 7.1 Determination of the Stability Boundary V-Function Value, form the Postfault System Power-Angle Curve 51

PAGE 10

x 7.2 Relationship Between X and the Mechanical Shaft Power Pm12 ............................... 52 7.3 Effect of the Equivalent Mechanical Power (Pm12) on the Phase Plane Postfault Stability Boundary.. 57

PAGE 11

CHAPTER I INTRODUCTION The dependence of modern society on electrical energy requires that major power failures be avoided. The loss of synchronism between generators and concentrated loads, caused by power surges or system faults, presents a potential cause of power failures. In order to evaluate the hazards and to take steps to prevent loss of synchronism, accurate methods of analysis of on-line system stabil ity must be developed and put into use. Current practice usually involves analytical studies which result in system design or operating procedure modifications. Even with the use of modern digital computer modeling, this approach does not fully protect the modern large scale system during emergency fault conditions. This thesis considers the simulation of a proven system, which is modelled by a two-machine equivalent system, during possible fault and erratic operating conditions. A goal is to seek, by analytical investigation, the range of operating limita tions for continual safe and reliable operation of such a power system. The basic approach used involves a computer stability simulation of the composite generator groups and equivalent loads

PAGE 12

2 A comprehensive presentation of the class of problems involved and the analytical models currently in use, together with selected models involving a stability analytical tool known as the second method of Liapunov, is given. This method together with the complex problem of developing Liapunov functions is developed and compared by example with the classical semi-graphical methods in current use, for relative accuracy and convenience in application. Transient stability involves the study of the system's capacity to recover, after a disturbance, to synchronous steadystate operation. Such a study should consider the action of voltage and speed regulators which would normally require the use of Park's transformation[14] and also a consideration of variations in mechanical power. In addition, damping parameters which complicate the analysis, should be included. However in this thesis a first approach to the problem which gives conservative results is used. This is an extension of the standard procedure used in transient stability studies known as the step-by-step method. The problem may be formulated as follows: given a system initially in a steady operation and assume a disturbance at time to. Then the question: is there a stable equilibrium position for the system after the disturbance is cleared, and if so, what is the critical clearing time, that is, the maximum time that the perdisturbances may remain, before the system loses its capability

PAGE 13

3 to return to steady-state. The advantage of the use of Liapunov's second method as developed in this thesis, is that the extensive computing time required in the.last phase, is replaced by a stability criterion (which is dependent on the construction of a Liapunov function) which allows the determination, in the state space, of the stability domain which surrounds the stable asymptotic solution of the post-fault system. This method is discussed in this study. The success of the approach depends upon the development of the Liapunov function, which is a main concern of this thesis. It is essential to develop the best possible Liapunov function for the model under study. This function is developed to produce the largest possible asymptotic stability domain. The method selected for the generation of the Liapunov function for the two-machine equivalent system is the Energy Metric Algorithm as extended by Mansour This algorithm permits the Liapunov function to be generated directly from the system equations. This approach is very valuable in power system stability studies.

PAGE 14

CHAPTER II STABILITY 2.1 Introduction Power system stability concerns the operating condition in which the various synchronous machines in a power system remain in synchronism, or in equilibrium with one another. A stable electric power system returns to normal operation after a disturbance or change in the operating condition, the most severe type of disturbance to which a power system is subjected is a three phase short circuit. Stability studies will be divided into three different categories depending on the extent of the disturbance on the system : 2.2 Steady-State Stability Steady-state stability is the ability of the synchronous machines of a system to remain in synchronism with each other subsequent to the occurrence of small disturbances. These small disturbances include very small gradual load changes which can occur on a continuous basis during the system's normal operations. This type of stability is realized by the use of manual control devices, since the majority of the excitation systems in use have

PAGE 15

5 sufficient deadband and long enough time delays so that the voltage regulator action ( E fd= 0 ) can be neglected for small disturbances. The steady-state stability limit refers to the maximum flow of power possible through a particular operating point without the loss of synchronism and when the power is increased very gradually. 2.3 Dynamic Stability Dynamic stability concerns the ability of the generators of a large interconnected system to remain in synchronism over long periods of time following small changes or normal random disturbances in the dynamic system. The stability limit is greater than the steady state limit. This type of stability can be controlled by the use of automatic control devices and voltage regulators. Since the excitation systems are used without deadband, the voltage regulator action 0) cannot be neglected for either small or severe disturbances. In dynamic stability analysis the system equations may be linearized about the operating point. Stabil ity of the linear time-invariant systems will be determined by use of the following methods: RouthHurwitz criterion, root-locus, and Nyquist cri terion. Using the state space form X=AX+Bu (2.1)

PAGE 16

6 the stability characteristics of a system may be determined by examining the eigenvalues of the! matrix. In the state space equation above, X is an n vector denoting the states of the system and A is a coefficient matrix The system inputs are represented by the rth vector u, and these inputs will be related to the differential equations by an n x r matrix B. This description has the advantage that! may be time-varying and that may be used to represent several jnputs including system load responses. 2.4 Transient Stability Transient stability is the ability of the machines to remain in synchronism subsequent to the occurrence of large disturbances .in the system. The disturbance may be a sudden increase in load, or switching operations which result in the loss of a major system component, such as a large generator or a long transmission However, the most severe type of disturbance is the three phase fault. The limits of stability in these cases are defined to be the maximum flow of power possible through a given operating point without the loss of synchronism at the time of the disturbance. The system equations for a transient stability study are nonlinear in nature and the study of the first swing cycle is most important. The system is described by a set of coupled non-linear differential equations of the form : X = f( t) (2.2)

PAGE 17

7 Where f is an n vector of nonlinear functions. The determination of the dynamic behavior of the system described by Eq. (2.2), is more difficult than that of a linearized system such as for Eq. (2.1). Digital computers are extensively used to analyze the stability off line for reliability studies but difficult to use on line for emergency decisions during faults. 2 5 Methods of Simulation The first step in a stability study is to make a mathematical model of the system during the fault. The elements included in the model are those affecting the dynamics of the machines (acceleration or deceleration of the machine rotors). Generally, the elements of the power system that influence the electrical and mechanical torques of the machines should be included in the model. These elements are: 1. The network before, during, and after disturbance. 2. The parameters of the synchronous machines. 3. The loads and their characteristics. 4. The excitation systems. 5. The mechanical turbines and speed governors. 6. System components such as transformers, and capacitors. The complexity of the model depends upon the type of transient and system being investigated. Using the techniques of modern control theory, the system stability of a large nonlinear system can be determined without obtaining the solution of the

PAGE 18

8 system of differential equations. that is, the stability limit can be calculated directly from the system equations and used during system operation. In the subsequent chapters, the close relationship of three different methods of analyzing stability will be discussed. These are the classical equal-area criterion, the phase plane trajectories technique, and the second method of Liapunov. These will be demonstrated for the case of two connected finite machines. Of the three methods, it will be shown that the second method of Liapunov is exact for a two machine system but at present conservative for three or more machines. These will be carried out for an example involving the transient stability region for the system using the three different methods.

PAGE 19

CHAPTER III CLASSICAL MODEL 3.1 Introduction Transient stability studies are performed to determine if a system will remain in synchronism subsequent to the occurrence of a major disturbance. That is, will the machine rotors remain in synchronism and will they return to a constant speed of operation following the disturbance. 3.2 Classical Stability Study Assumption The basic assumptions usually made in a classical stability study are: 1. Mechanical power input to each machine is held constant. 2. Damping is neglected. 3. Loads are represented by passive impedance. 4. Each machine has a constant-voltage behind transient reactance. 5. The mechanical rotor angle of each machine coincides with the angle of the voltage behind transient reactance. 6. Machine saturation is neglected. The main reason for making the previous assumptions, is that the calculations are simplified and, thus, attention can be focused on the main problem of determining the stability region

PAGE 20

10 within a reasonable degree of accuracy. References [1] and [13] consider the detailed reasoning involved. 3.3 Power System Representation. The equivalent two-machine power system of Figure 3.1 will be used as an example (problem 2 .18 [1]) in the comparative analysis of the stability boundaries for the three methods. The computations are based on a two-finite machine system model in which the classical representation is used for each machine. In all cases resistances are neglected. The initial internal voltages of the generators, the initial torque angles, and the load equivalent passive admittance matrix are determined from the load-flow data. In transient stability studies a load-flow calculation is first made to obtain system conditions prior to the disturbance In this calculation all power system data is converted to a 100 MVA Base for the. two generators transmission lines, and system components. The load-flow data are given in Tables 3.2, 3.3, and 3.5 respectively. 3.4 System Equation Of Motion During transients in a two-finite machine system, the system goes through the following stages: (a) System before the disturbance or fault; prefault system (b) System during the disturbance or fault; faulted system. (c) System after the removal of the disturbance or fault; post-fault system.

PAGE 21

11 SUS1 BUS2 6 suss ZLS Figure 3.1 Single-Line Diagram of Two-Machine Power System Table 3.1 Generator Data (in p.u to Generator MVA Base) GENERATOR NO. H(MW.S/MVA) RATING(MVA) 1 0.28 0.080 5 0 50.0 2 0.250 0.070 4.0 120.0 Table 3.2 Transmission Line Data on 100 MVA Base (Resistance Neglected) Line Number 3 4 5 6 X p.u(100MVA) 0.08 0.06 0 .06 0.13

PAGE 22

12 Table 3.3 Generator Data (in p.u to 100 MVA Base) GENERATOR NO. XT(P.U) H(MW.S/MVA) RATING(MVA) 1 0.56 0.16 2.5 0.5 2 0.21 0.05 4.8 1.2 Table 3.4 Load-Flow Data VOLTAGE LOAD GENERATOR -BUS NO MAG. ANGLE MW MVAR MW MVAR P.U DEGREE 1 1.03 0 0 0 30 23.1 2 1.02 -0.5 80 40 100 37.8 3 1.018 -1.0 50 20 0 0 Table 3.5 Load-Flow Data (on 100 MVA Base) VOLTAGE LOAD GENERATOR -BUS NO MAG. ANGLE MW MVAR MW MVAR P.U DEGREE 1 1.03 0 0 0 0.30 0.23 2 1.02 -0.5 0.80 0.40 1.0 0.37 3 1.018 -1.0 0.50 0.20 0 0

PAGE 23

13 The form of the swing equation of the two-finite machine the system during the previous three stages will be the same as that given by Eq (3.1), except for the difference in the input output power from one stage to another. Thus, the swing equation for machine one is 2 M1 d 01 = Pal = Pm1 Pel dt2 The swing equation for machine two is d202 M2 --= Pa2 = Pm2 Pe 2 dt2 (3.1) (3.2) (3.3) The relative torque angle0 = 01-02 will be used, because its value is significant in showing the stability or instability of the two machine system. d20 Pal Pa2 dt2 ="Ml M2 (3.4) Next, multiply each side of Eq.(3.4) by the equivalent inertia constant. M = M1M2 M1 + M2 2 M = Pa1M2 dt2 M1 + M2 PUM1 M1 + M2 M d20 = M2Pm1 -M1Pm2 dt2 M1 + M2 M2Pe1 M1Pe2 M1 + M2 Which may be written more simply as in Eq. (3.1) d20 M --= Pm12 Pe12 dt2 (3.5) (3.6) (3.7) (3.8)

PAGE 24

14 Where Pal2= the equivalent acceleration output power Pm12= the equivalent mechanical input power Pel2= the equivalent electrical output power M = the equivalent inertia constant o = the relative internal voltage angle between the two machines(ol 02). 3.5 Prefault Power System Model The single-line diagram of Figure 3.3 is a further simpli-fication of the system of Figure 3.2. The elements of the bus admittance matrix for the network are reduced by delta-wye trans-formation to two nodes in addition to the reference node. From the load-flow study of the pretransient network, the following phasor values are determined: '100 El 1.21 eJ P.U E2 = 1.15 e j12.650 P.U 1. 789 + j 0.716 P.U ZL3= 1.039 + j 0.521 P.U Pml= 0 3 P.U Pm2= 1. 0 P. U The Y matrix for the reduced network is o 064 -j 1. 001 0.182 + j 0813 Y t ma rIX = 0.182 + j 0.813 O. 528 -j 1. 358

PAGE 25

15 Pt'1 __ .... X12 .--Pt'2 E IL1>...L I Ml'23 XU-5 I M2''f:! XI,,-5 ZL3 X'dI-0.56 XTl.0.16 XI2-0.13 X'.I2'0.21 XT2oB.B:S9 XB210.08 Xls-e.0e ZL2tl.9S fLall.IS (ALL VALUES IN PER UNIT ON R lee MVR BRSE) Figure 3.2 Pre-Fault Reactance Diagram -0.182-J0.813 0.710 El E2 -J0.188 -J0.545 Figure 3.3 Pre-Fault Admittance Diagram

PAGE 26

16 The equation of electrical power output from a machine is given by n P ei = G .. + [ E. E. Y .. cos( 6 .. C. + C. ) 1 11 -1 J 1J 1J 1 J i=l j;ti So that, 2 Pe1 = E1 G11+ E 1 E 2 Y12 cos ( 612 C1 + C2 ) Pel = 0.094 + 1.16 cos ( 77.4C1 + C2 ) 2 Pe2 = E2 G22+ E 1 E 2 Y12 cos ( 612 C2 + C1 ) Pe = 0.698 + 1.16 cos ( 77.4C2 + C1 ) 2 Since P PelM2 Pe2MI e12 -M1 + M2 Where -5 2 -5 2 M1=23X10 sec ./elec deg., M2=45X10 sec ./elec deg. Then Pe12 = -0.177 + 1.135 sin ( C12 + 4.1 ) (3.9) (3.10) (3.11) (3.12) (3.12) (3.13) (3.14) The maximum power that can be transmitted under rated operating condition occurs at a torque angle of( C12 + 1 ) = _90, where 1 = 4.1 and C12 = -94.1. Notice that the characteristic is a sine function of C between the two rotors as plotted in Fig-ure 3.4. This sine curve is displaced from the origin vertically by an amount Pc, which represents the power dissipation in the equivalent network, and horizontally by the angle 1, which is the real component of the transfer admittance Y12. The system is stable only if the torque angle C is in the range from -900to +90, in which the slope dP/dc is positive, so that an increase in torque angle results in an increase in transmitted power.

PAGE 27

1 0.5 -c: ::J .... Q) 0.. :. 0 N Pc Q) a.. a: w -0.5 3= 0 a.. -l c:x: -1 w a: -1.5 r
PAGE 28

18 So the maximum power resulting is Pe12 = -0.177 + 1.135 sin(-94.10) = -1.31 p.u max. (3.15) Pe12 = -1.31 p.U occurs at D12 = -94.10 ( -1.642 rad. ) max. And the steady state synchronous speed, Pm12= Pe12, is Pm12 = = -0.1452 p.U -0.1452 = -0.177 + 1.135 sin ( DO +4.1) arcsin ( DO +4.1 ) =0. 0318 DO = _2.5 (-0.0436 rad.) 3.6 Power System Model During Fault (3.16) The single-line diagram for the faulted power system is shown in Figure 3.5 and 3 6. The disturbance to be considered is a three-phase short circuit which occurs near bus number 3 at the end of line 5. The line-relaying schemes will sense the fault on the line and the fault will be cleared in six cycles by simulta-neously tripping of the line terminal breakers. Using the same power output equation given in (3.9), and the bus admittance matrix for the reduced network during the dis-turbance the Y matrix for the reduced network is 0.0001j 1. 295 Y t ma rIX = 0.001 + j 0 044 The new system equations are: 2 0.001 + j 0.044 0.006 -j 3.419 Pe1= El GIl + E 1 E 2Y12 cos (012 Dl + D2) (3.17)

PAGE 29

1 -....... E1L1>-L X'dl.56 x'd 2' 0.21 X13=0.08 (ALL VALUES X12 ZLa XTl.16 XI2'0.13 XT2S.S38 X3280.06 ZL2. 93 ZL3z1. 16 IN PER UNIT ON A 100 MVA BASE) 19 ... Xd2 E2L12-M2''I5 Xls-5 Figure 3.5 During Fault Sys tem Reactance Diagram -0.00077-J0.0Lf'i 0.0060 El E2 -Jl.25 -J3.38 E 1 = 1 2 1/1 P. U E 2 = 1 1 5 It 2. P. U Figure 3.6 During Fault Admittance Diagram

PAGE 30

Since Then Pe1= 0.000015 + 0.0606 cos ( 890 01 + 02) 2 Pe2= E2 G22+ E 1 E 2 Y12 cos ( 812 02 + 01 ) Pe = 0.0079 2 o + 0.0606 cos ( 89 02 + 01 ) Pe = Pe1M2 Pe2M1 12 M1 + M2 Pe12 = -0.0026 + 0.061 sin ( 012 + 1.00 ) 3.7 Post Fault Power System Model 20 (3.18) (3.19) (3.20) (3.21) (3.22) The post-fault system single-line before and after network reduction, is indicated in Figures 3.7 and 3.8. The mag-nitudes of E1 and E2 computed for the prefault system operation condition are used in the postfault, and fault model, to find the electrical output power Pe12, according to the previous assump-tions. 01 and 02 are time dependent variables which are changed for each time interval. When the fault is cleared by simultaneous opening of the line end breakers, the power-angle equation applies, since a net-work change has occurred. This requires the use of Eq.(3.9) and the bus admittance matrix for the reduced network post-disturb-ance. The Y matrix for the reduced network is : 0.0667j 0.999 0.186 + j 0.807 Y t ma rlX = 0.186 + j 0.807 0.519 -j 1.347

PAGE 31

21 1 __ X12 __ PII' 2 El& I Mla23 X18-5 E1 XTl.1S X12.13 X32.0S ZL2.93 ZL3-1.1S X'diI0.5S x'd 2 10.21 XI3.08 (ALL VALUES IN PER UNIT ON R 199 MVA BASE) Figure 3.7 Post-Fault System Reactance Diagram -0. 185 -J 0. 8 10 0. 70Lt -J0.19 -J0.5Lt E 1=1. 21l1i P.U. E.2 = 1 1 5 It 2 S P. U Figure 3.8 post-Fault Admittance Diagram E2

PAGE 32

And the system equations are 2 ( Pe1= E1 G1l+ E 1 E 2 Y12 cos Pel= 0.09760 + 1.15 cos ( Pe2= G22+ E 1 E 2 Y12 cos ( Pe2= 0.6860 + 1.15 cos ( Since Pel2= PelM2 -Pe2Ml Ml + M2 Then Pe12 = -0.167 + 1.120 sin (Jl2 81+ 82 ) 77 81 + 82) (Jl2 82 + 81 ) 77_ 82 + 81 ) ( 812 + 4.3 ) 22 (3.23) (3.24) (3.25) (3.26) (3.27) (3 28) The steady-state values of Pe12, 8 and d 28/dt2 for the stable postdisturbance system are Pm12 = Pe12(steady-state) = -0.1452 p.u o -0.1452 = -0.167 + 1.12 sin ( 8ss + 4.3 ) arcsin ( 8ss + 4.3 ) = 0.0218 8ss = -3.2 (-0.0558 rad.), and d 28ss/dt2 = 0 From the power-angle diagram shown in Figure 5.1. 8max. = 8ss 8max. = -176.8 (3 0857 rad.)

PAGE 33

CHAP1'ER IV NUMERICAL SOLUTION OF THE SWING EQUATION 4.1 Introduction The analysis of first-swing classical transient stability constitute the important tools for judging system performance. The reason for its relative importance is that if the system is stable on first swing, it will for most cases be stable on the subsequent swings. The torque angle 0 is calculated as a function of time over a period long enough to determine whether 0 will increase without limit or reach a maximum and start to decrease. 4.2 Point-by-Point Solution of the Swing Equation The swing equation governing the motion of each machine of the power system, is (4.1) Where 0= displacement angle of rotor with respect to a reference axis rotating at synchronous speed M= inertia constant of machine Pa= accelerating output power t = time

PAGE 34

24 The solution of Eq.(4.1) gives 8 as function of t, which is known as "swing curve". Examining the' swing curve of the two machine system will show whether the machines will remain in synchronism after a disturbance. In a two-finite machine system, the output and hence the acceleration power of each machine depend upon the angular posi-tion 8 and d8/dt. Thus, for a two machine system, the two simul-taneous differential equations are: d281 M1 --= Pm1 Pe 1 ( 81, 82 ) dt2 2 M2 d 82 = Pm2 Pe2( 81, 82) dt2 Analytical solution of these equations are not possible and (4.2) (4.3) solved numerically. The method is called step-by-step or point-by-point solution. In deriving the swing curve of a system under study, step-by-step solution becomes very handy comparing to some of the methods recommended for digital computers. In this method, the change in the angular position of the rotor during a short interval of time is computed by making the following assumptions: 1. The accelerating power(Pa) computed at the middle instead of the beginning of the interval to eliminate the error caused by Pa which is constant during a time interval dt. 2. The angular velocity (w)is constant throughout any interval at a value computed for the middle of the interval.

PAGE 35

25 These assumptions are justified by deereasing the time interval; as the time interval is decreased, the eomputed swing eurve ap-proaehes the ture eurve. Figure 4.1 demonstrates the above assumptions. If Eq.(4.1) is divided by M and integrated twice with respect to t, Pa being treated as a constant, we obtain the following Then and also de Pa t dt = W = Wo + M Pa 2 e = eo + wot + M t fl t fl wen -1/2)= M Pa(n 1) wen -1 /2 ) = wen -3 /2) + fl wen -1 /2 ) fl en = flt. wen -1/2) (flt)2 fl en = flt. wen -3/2)' + M Paen -1) (flt)2 fl en = flee n -1 ) + M Paen -1) en = een 1) + flen (4.4) (4.5) (4.6) (4.7) (4.8) (4.9) (4.10) (4.11) Eq.(4.10) is the main form of the step-by-step solution of the swing equation. The equation shows how to calculate the change in e during an interval if the change in e for the previous interval and the aceelerating power for the interval in question are known. The accelerating power is calculated at the middle of each new interval. The solution progresses through enough interval to obtain points for plotting the swing eurve, whieh is suffieientlyaccurate[2].

PAGE 36

Pa Pan-2 (On -1 /2 cq,-3/2 On -1 On-2 I I --_,_ --1__ -L----("J..::---, I I I I I -"\ -"1 --"1 -+---0,---. I I I n-2 n-3/2 n-1/2 n -1 n t .1t t .1t Figure 4.1 Actual and Assumed Values of Pa, W, and C as Function of Time 26

PAGE 37

27 4.3 Numerical Methods for Solutions of Differential Equations The methods most commonly used for the solution of the differential equations are: Euler method, the modified-Euler method, Runge-Kutta and the trapezoidal method Each of these has advantages and disadvantages which are associated with numerical stability, time-step size, computational effort per integration step and accuracy of the obtained solutions. Appendix A shows the computational printout for plotting the swing curves of the two-finite machine system. For clearing at 0.1 second the solution is obtained by use of Turbo Pascal program with the modified-Euler method procedure. So the numerical solution of the swing equation for the two generator, three-bus power system is made by digital computer for 2.0 second of simlated real time, for the intervals of 0.05 sec. Figure 4.2 shows the rotor angles of the two machine vs time A plot of the rotor angle is shown in Figure 4.3 and the fault is cleared in six cycles. It follows that the system is stable. The rotor angle difference reach value of -11.67 and then decrease. This is the value of 812 at t= 0.2sec. Note that the solution is carried out for three swings to show that the subsequent swings are not greater than the first so that the system appears to be stable. But in the case of the angle differences increase indefinitely, the system is unstable because both machine will lose synchronism.

PAGE 38

400 (/) Q) .--l 0 350 0 300 c .-1 Q) 'd 00

PAGE 39

6 rJ) Q) .-i CJ 4 CJ 2 -" Q) Q) ... 0 as' "0 -c.o 2 W ..J C) 4 z w 6 a a: 0 I-8 -1 0 1 2 0 0.2 I 0.4 0.6 0.8 1 1.2 1 .4 TIME, second I I I I I 1.6 1.8 2 I 1 2 24 36 48 60 72 84 96 108 120 Cycles Figure 4.3 Plot of Delta Differences (01-02) VS. Time

PAGE 40

CHAPTER V STABILITY DOMAIN 5.1 Introduction As mentioned earlier, there are at least three different methods which are commonly used to study the stability of a power system during a fault. The first method is out-lined below, and the others are discussed in the next two chapters. Appendix B contains further development. 5.2 Power-angle Characteristics and the Equal-area Criterion For the two-finite machine system under study, the power angle characteristics are shown in Figure 5.1, for condition before,during, and after a three phase fault. The horizontal line denoted by P=Pm12 represents the equivalent mechanical power input to the machine. Before occurrence of the fault, the two machine were operating at synchronous speed with a rotor angle at t=O is CO = -2.5 degree as indicated by the intersection of Pm12 with the prefault curve, this operating point(Co ,Peo) is designated by the letter a. Once the a-phase fault has occurred, the electrical power Pe12 of the system increases to a value corresponding to point b, at which Pa12 = (Pe12(b) Pm12), this results in a decrease in both d8/dt and C. As C continues to

PAGE 41

0.2 0 ....... ... s:::: L Q) -0.4 a. s:::: ........ N -0.6 ,.... Q) a.. w 3: -0.8 0 a.. --I Omex w -1 -1.2 c -1.4 -180 -160 -140 o PRE. b ,. " s:: u 55' I I -120 -100 -80 -60 -40 -20 o TORQUE ANGLE (in degree) o DURING POST ... Pm 12 =-0 .1,45 p.U. Figure 5.1 Classical Power-Angle Curve

PAGE 42

32 decreases the system remains disturbed, the power angle trajectory moves along the faulted power angle curve form point b toward point c, the fault is cleared, at which time for this case, Pel2 decreases to 0.85 p.U ( oc = -143 degrees as indicated by point d in the figure). The area Al in the figure is proportional to the kinetic energy (K.E.) of the system during the period when the fault occurs and is cleared. As 0 continues to decrease, the power angle trajectory moves along the postfaulted power angle sine curve from point d toward point e for which Omax. is -176.8 degrees. When this point is reached, do SQR[ -!-J :: .. dol = 0 (5.1) dt At this point the K.E. of the equivalent finite machine has returned to its prefaulted value, K.E = 0 At point e, the area A2 bounded by Pel2 = Pml2 and Pel2 = Pc + Pmax. sin (0 + ,(post)], 0 0 max., is equal to the area Al bounded by Pel2 = Pml2 and Pe12= [Pc + Pmax; sin(o + ,(during)] 00 0 Oc There exists a (o=Oc), (00 0 omax), such that in terms of the equal-area criterion of stability, Area Al is equal to A2. and the system is stable, the value of oc is determined as follows:

PAGE 43

33 Oc Al = i[Pm12 -(Pc + Pm.x. sine + 7(during,) d. (5.2) 00 Let x = 0 + '"/ dx = do (Oc + '"/) Al = J[-0.1452 + 0.0026 0.061 sin x ) dx (00 +'"/) (5.3) AI= -0.1426 Oc 0.0547 + 0.061 cos ( Oc + 0.330 ) (5.4) Omax. A2 =i [(Pc + Pm.x. sin ( + 7(post, Pm .. ) dO (5.5) Oc (omax.+ '"/ ) A2 = i [ -0.167 + 1.12 sin x + 0.1452 ) dx (oc+ '"/) (5.6) 0 A2 = 0.0218 Oc + 1.177 + 1 12 cos ( Oc + 4 .3) (5.7) Al = A2 :::} -0.1644 Oc 1.056 cos Oc + 0 084 sin oc -1.232 = 0 (5.8) From the power-angle diagram shown in Figure (5.1) the critical clearing angle is located between 00 and omax. Since Eq.(5.8) is nonlinear, oc= -143 degrees (-2.495 rad ) was found by using trial and error.

PAGE 44

" 34 If the fault clearing is delayed long enough so that the equality of the two area cannot be satisfied, the two-finite machine speed will not decrease to a synchronous value as long as the machines remain electrically tied to each other. The torque angle S will decrease monotonically without bound beyond the maximum value possible for a marginally stable swing, Smax. By the time S reaches a value of -180 degree, synchronism is lost and the machines must be disconnected[51.

PAGE 45

CHAPl'ER VI PHASE PLANE TRAJECTORIES AND THE STABILITY BOUNDARY 6.1 Introduction This technique provides a useful tool for studying the stability of a system which is described by a second-order differ-ential equation or a group of such systems. The swing equation of the power system ( prefault, during fault, and postfault) again is a second-order differential equation of the form: d26 M --+ Pe12 ::: Pm12 dt2 (6.1) To form the phase plane, the above equation is converted into two first-order differential equations with the time t suppressed. Let Where Xl ::: X2 X2 Pm12 Xl ::: 6 do X2 --dt Pe12 M if> and if> ::: 6ss + "f (6.2) (6.3) From the steady-state values of Pe12, 0 and d20/dt2 the stable postdisturbance system is Pm12 ::: Pe12 (steady-state) Pe12 ::: Pm12 ::: Pc + Pmax. sin (6ss + "f ) (6.4)

PAGE 46

oss =[ arcsin ( Pml2 Pc ) Pmax. 'Y ] 36 (6.5) o oss = -3.2 (-0.0558 rad.), and d 2oss/dt2 = 0 ( as found previously for the post-fault power system model ). 6 2 Phase Plane Trajectory A phase plane plot is a plot of X2 vs. Xl, as shown in Figure 6.1. As the time t increases, the two-tuple [XI(t), X2(t)] describes the trajectory in the phase plane. Since stability of the postfault system is the basic requirement, finding the state equilibrium point ( the origin of the phase plane) is important, since this point is the steady state value. The state equation for the postfault system is Xl = X2 (6.6) X2 = [ Pml2 Pc -Pmax. sin (Xl + ) ]/ M (6.7) Since Xl = X2 S ( XI,X2 ) = X2 (6.8) X2 = [ Pml2 Pc -Pmax. sin ( Xl + ) ]/ M Q( XI,X2 ) = [Pm12 Pc -Pmax. sine Xl + )/ M Then the state equilibrium singular point is: S( XI,X2 ) = 0 X2 = 0 Q( XI,X2 ) Pc -Pmax. sin(xi + M= 0 ( Xl + ) = [ arcsin ( Pml2 Pc ) ] = Pmax. where n = 0, 1, 2, (6.9) (6.10) (6.11) (6.12)

PAGE 47

8 7 6 5 U 4 Q) CI) 3 Q) e 01 2 Q)-.. "OCI) 1 c"O C NO 0 Q) 0 =: J _.c: .!It 0 1: Q) -3 I +' 0 -4l +' Ul "a,b Point Fault TT.iectory\ -51 -6 j -7 -8 -190 -170 -150 -130 -110 -90 -70' -50 -30 -10 10 State Variable X1. degree Figure 6.1 Phase Plane Plot Of X2 Vs. Xl

PAGE 48

38 the singular point is chosen for n = 0, which defines the state variables, Xl and X2, such that the origin of the phase plane is the stable equilibrium point of the postfault system. The equilibrium states of a system are defined as: as as aXI ax2 axl aX2 It is necessary to evaluate the equilibrium states of the system at singular points. This is done by testing the stability of the system by determining the nature of the roots. If positive real roots exist, the system is unstable for the given operation conditions. if no positive real roots exist the system is stable. o 1 Pmax. cos D o (0,0) Now, the stability of the system is from -,\ 1 = ,\2 Pmax. ,\ = VPmax. Pmax. -,\ (0,0) The result is two real roots with opposite signs, then point e on the phase plane is an unstable singular point being a saddle point singularity of the trajectories. The state equations for each system are

PAGE 49

Pre-fault system Xl = X2 X2 =[Pm12 (Pc + Pmax. sin(o + Oss + 1 ]/M
PAGE 50

,.... c :;, c C N -l. t K. I C "-I,.... t i 40 0.1 b 0 -0.1 c -0.2 -0.3 -0.-4 -0.7' -0.8 -0.0 .,.1 -1.1 -1.2 -1.3 I -1.-4 1 -I 11-1 I -180 -leo -140 -120 -100 -eo -eo -40 -20 0 20 I pocl TOfqlHl Angle, + I OIJRINO <> POST I 1 A I Dmnx. \,. 69s t.. 6min. 8 7 I II I II I -4 I J I I 2 Saddle e I 0 f Point -1 I -2 I Faul t T.raj ectory -3 I --4 I -II -8 -7 --------_8 -100 -170 -I!!o -130 -110 -PO -70 -!!O -30 -10 10 state Variable X 1. Figure 6.2 Power-Angle Curves And Phase Plane Trajectories

PAGE 51

41 which time the faulted line is isolated from the system. For fault clearing torque angles greater than oc, the postfault system will be unstable. At the time the fault is cleared (at t=0.1 sec. corresponding to o=oc) the output power Pe12 decreases from the value of -0.5 to -0.85 p.u, 0 continues to decrease along the phase plane post-fault trajectories from point c and d (at which o = oc) toward point e. Point e (0 = Omax.) is an unstable singular point which is a saddle point of the trajectories. For the conservative system study, 0 increases a path is then formed along the phase plane maximum trajectory in a clockwise direction from point e to point f as show in the Figure. Referring to Figure (6.2), omin can be determined from the equal-area criterion, using postfault power-angle curve such that: Let (OSS oss omin -( Pc + Pmax. sin(o + 1 t)] d 0 (POs ) omax J[ (Pc oss x = 0+1 (max + [-0.167 + 1 ) + Pmax. sin(o + 1 ) (post) dx = do 1 + 1.12 sin x + 0.1452]d x omin = 0 19 (0.332 rad.) -Pm12] d 0 = -2.16 p.u (6.19) (6.20) (6.21) (6.22)

PAGE 52

The critical clearing time t can be determined from cr 42 Eq.(5.1), using numerical integration techniques: Where Pm12 -Pe12] d8 ] Pe12 = [ Pc + Pmax. sin ( 8 + 1 d )] ( ur) 8 8 JPa12 do = J[Pm12 (Pc + Pmax. sin(o + 80 8 (6.23) (6.24) (6.25) (6.26) t _j_8C ______ ----""d.;;..8 _________ (6 27) ::[(Pm12 Pc) (0-00) + Pmax.(coso cosoo)] t 0.1 second cr 6.3 Stability Boundary The system is considered stable as long as the trajec-tories follow the separatrix determined by the phase plane (as shown in Figure 6.1). The equation for the separatrix can be determined from the post-fault system differential equation of motion :

PAGE 53

43 d 2 6 Pal2 M--= dt2 = Pml2 Pel2 (6.28) Since -L.[..MJ' = 2 [ :!, ] d [ :! ] dt dt dt 2 d 2 6 d. [:!] = 2 d6 dt2 (6.29) Then 1 d [!: r = d 2 6 2 d6 dt2 (6.30) Substituting Eq.(6.30) for d 26/dt2 in Eq.(6.28), the differential one-form is obtained. -.JL _d_ [ d6] 2 = [Pm12 (Pc + Pmax. sin 6)] 2 d6 dt (6.31) d. [ d6]2 = [Pm12 (Pc + Pmax. sin 6)] d6 dt M (6.32) In term of the state variable, Xl and X2, d(X2) = (Pc + Pmax. sin(XI + dXI 2 M (6.33) Integration is performed to obtain the general solution for a postfault system phase plane trajectory: X2 PC(XI) + Pmax. COS(XI + + C = 0 2 M The constant of integration C is evaluated at the saddle point of the separatrix, which is ( 6max 0 ). C= Pc)(6max + Pmax. cos 6max] M (6.34)

PAGE 54

44 The stability boundary of thepostfault system, is X2 PmI2(XI) PC(XI) + Pmax. cos( Xl + )] 2 M + ) + Pmax. cos = 0 (6.35) M and X2 + [( Pml2 Pc )( Xl ) 2 M + Pmax. [cos cos (Xl + = 0 (6.36) ( ) Xl ( ) (6.37) The results obtained by this method agrees with that deter-mined by the equal-area method. The postfault system trajectories illustrated in Figure (6.3), were generated by varying the fault clearing time (tfc). If the fault is cleared in time at which t f < t I operation will c -cr be stable along a curve as indicated in the figure for the three stable phase plane trajectories, being confined to the phase plane region enclosed by the stability boundary. For a sustained fault or for a longer delay in clearing time at which tf)t I the c cr faulted trajectory will enter the unstable region as shown in the figure for the two unstable phase plane trajectories by the dotted curve. Since the power syst,em under investigation is a conserva-tive system (no damping), the stable trajectories are closed tra-jectories. For stable trajectories in a non conservation system X(t) 0 as t 00.

PAGE 55

('of o 00 ,N cp-CU 'C) cuo ""'00 o "0 No X o CUI -. ,-0 coN stable { unstable {t fc = 0.115 sec. t f c = 0.12 sec. t f c = t c r = 0.1 sec. t fc = 0.09 sec. t fc = 0.08 sec. t fc = 0.07 sec. State Var I ab I e Xl, In dearee *101 Figure 6.3 Plot of Stable and Unstable Phase Plane Trajectories

PAGE 56

CHAPTER VII ASYMPTOTIC STABILITY 7.1 Introduction Liapunov's "second"(also called the "direct") method has been used with varying degrees of success in analyzing both linear and non-linear stability criteria. A desirable feature of the second method is that it will permit the determination of system stability without obtaining solution of the system differential equations. 7.2 Liapunov stability The Liapunov second method permits stability analysis of system differential equation, of the form : X. = F. (X) -1 1(7.1) Where X is the postfault system state vector ( x1,x2"'" xn ) determined from the differential equations of motion which describe the behavior of the power system. There are four conditions which need to be studied to analyse the stability of a system : 1. V (!) and its first partials are continuous in where is a region about the origin, and V(!) is a scalar function. 2. V(Q) = 0

PAGE 57

47 3. vex) > 0 for all! E R, X 0 4. vex) 0 for all! E R, X 0 If the four conditions are satisfied, then Eq.(7.1) is stable. If condition 4 is replaced by : V(!) < 0 for all! E R, 0, then Eq. (7.1) is asymptotically stable in the region R. If there exists a scalar function V(!) which satisfies the above conditions for all X in state space, then Eq.(7.1) is globally asymptotically stable. If a Liapunov function does not exists in R, the system is unstable. For a particular system, there is no general method to determine the existence of a suitable Liapunov function. A useful method for determining the V-function for the system under study is the Energy Metric Algorithm(see appendix C). 7.3 Energy Met.ric Algorithm Applying this method to a second-order system( n = 2), we define the state variables, as before, such that X (0) = 0 Xl = 8 if> d8 X2 --dt and if> = 8ss + "f The resulting two first-order differential equations are the same as those developed in previous chapter (Eq.6.17 and 6.18). Let Xl = X2 (7.2) X2 =[PmI2 (Pc + Pmax. sin(8 + 8ss + "f t ]/M (7.3) (POs )

PAGE 58

48 Since n = 2, there is only one equation thus formed. dXI X2 dX2 = M [Pm12 Pc + Pmax. sinC Xl + c/J ] (7.4) After the cross multiplication and rearrangement of the terms, there results: m = X2 dX2 + __ l __ [pc + Pmax. sin(XI+ c/J) -Pm12] dXI (7 4) M A line integration is preformed to obtain the V-function. v= J J x (x= X ) 2 I I sin(XI + c/J) -Pm12] dXI+ X2 dX2 (7.5) o V= --l--[pc(Xi) + Pmax.[cos c/J COS(XI + c/J)]PmI2(XI)] M 1 2 +-X2 2 The system is stable for the equilibrium point under (7.6) consideration if the total derivative of V with respect to time is at least negative semi-definite. (A function is negative semi-definite if it is equal to zero at the origin and less than or equal to zero in a region about the origin). av' av VeX) = --Xl + X2 (7.7) axl aX2 1 V(X)= -[PC-PmI2 + Pmax. sin(XI + c/J)] XI+ X2 X2 (7.8) M

PAGE 59

49 Since XI= X2 and v= 0 Then 1 V(X)= ---[PC-PmI2 + Pmax. sin(XI + X2) (7.9) M Note that V(X), and V(!) are continuous for all X in state space. Furthermore, V(!) 0 for all X E and = o. Where The region of stability is found to be: : X E [ veX) < J --= V(!) I ( x = t:max _..I. 0 I v 0/ ,x2= ) = __ 1 __ [ Pc( 8max ) -Pm12( 8max ) M + Pmax. ( cos cos 8max )]+ 0 The stability boundary of the postfault system, is V(!) = 0 (7.10) (7.11) (7.12) 2 X2 + PmI2)(XI) + Pmax. [cos COS(XI + M 2 [ ] ---(Pc PmI2)(8max + Pmax. (cos cos 8max) = 0 M 2 2 [ X2 + --(PmI2':" Pc )(8max Xl) M + Pmax. [ cos 8max cos ( Xl + ) J] = 0 (7.13) ( 8min ) Xl ( 8max ) (7.14) This result is identical to that obtained using the phase plane technique, in the following sense the maximum trajectory, which is the stability boundary for the system, has been found by

PAGE 60

50 setting X2 = 0 and Xl = ( omax ) in Eq.(7.10). The result defines the region of stability for the equilibrium at 0= omax (saddle point). This saddle point determines the value of V for the maximum stable trajectory, since the stability boundary passes through the singularity. The same result was found using Eq.(6.36) which was developed by the phase plane approach. Eq.(7.14) defines the exact boundary between stable and unstable regions of the power system operation. Also from Eq.(7.12), which describes the closed surface with the greatest is the value of the Liapunov function VeX) on which boundary of the largest stability estimate. The value of corresponds to the shaded area as shown in Figure(7.1), is the maximum amount of stored potential energy in the post-fault system, or a given power system operating condition. The value of the constant can be found using postfault system power-angle curve = Shaded Area M or omax. l: ;J [Pc oss + Pmax. sin(oss + 1 t -Pm12] do (POs (7.15) Stability domains corresponding to the smaller of the values of Pm12, the larger will be the value of and the larger will be the phase plane region of stability. As shown in Figure (7.2) for a given postdisturbance system condition.

PAGE 61

8mex I "...., -0. .... c !.... Q) -0. c.. c '-'" C\I -0. ID c... w 3: -0. D c... ...J W -1. -1. -180 -160 -140 -120 -100 -80 -60 -40 -20 o POST TORQUE ANGLE (in degree) A Pm 12 :-0. 145 p.U. Figure 7.1 Determination of the Stability Boundary V-Function Value, X, form the Postfault System Power-Angle Curve

PAGE 62

25 20 15 C\I_ 10 (/)"0 Z C 5 ttl -(I) O::l a: J:: 0 -t:. z -5 -1 0 -1 5 20 2 -1.5 -1 -0.5 o REAL POWER PrTl1 2. IN PER UNIT Figure 7.2 Relationship Between X and the Mechanical Shaft Power Pm12

PAGE 63

53 7.4 A Suitable Liapunov Function and an Estimation of As indicated previously, the estimated is the value of the Liapunov function VeX) on the stability boundary. From Eq.(7.6) vex) is : vex) = --1--[PC(XI) + Pmax.[cos COS(XI + M PmI2(XI)] + __ 1 __ 2 Eq.(7.16) is positive definite, if (Pc PmI2)(XI) + Pmax.[cos COS(XI + > 0 setting, V(!)I(x =0 x =0) in Eq.(7.16). I 2 (7.16) (7.17) V(O,O) = Pmax. (cos -cos = 0 a minimum From Eq.(7.16), V has a turning pOInt at a critical point of the power system equation : Next set, V(!) I (X x =0) I 2 V(8max. 0 )=[(Pc -Pml2 )( 8max.) +Pmax. (cos -cos 8max.)] > 0 (7.18) Therefore, the stability boundary (surface) of the postfault system, given by a Liapunov function is : V(!) = 0 or + Pc) (8max Xl) M + Pmax.[cos 8max cos (Xl + 0, provided that surface V = = constant. (7.19)

PAGE 64

54 > > 0, are closed surfaces, since the Liapunov theorem is satisfied within the region given by Eq.(7.14). Also, 1 VeX) = ---[Pc -Pml2 + Pmax. sin(Xl + -X2) is the M solution for the power system of Eq.(7.2 and 7.3). The region is asympotically stable. Closedness can be proved by considering the behavior of the surface V = Let Xl= C = a constant, Eq.(7.16) becomes 2 Pml2)(C) + Pmax.[cos cos(C + --1--X2 M 2 1 2 V= ---X2 + fCC) 2 Then the equation becomes 1 2 ---X2 2 = fCC) (7.20) (7.21) Eq. (7.21) describes the intersection of the surfaces with the plane Xl= C. For any > > 0, as C increased from zero, to the right hand side of Eq.(7.21) is initially positive, and the intersections are closed ellipsoidal curves. As C increases, the plane Xl= moves along the negative Xl-axis, and f(C)] wi 11 eventually become zero, and the curve of intersection vanishes, which indicates the closedn'ess of the surface For the case when closure is on the Xl-axis for f(C)= i.e., f(C)]= 0, or

PAGE 65

55 __ l __ [(PC Pm12)(C) + Pmax.[cos cos(C + (7.22) M The two roots of Eq.(7. 22) give two points: (a) Principal node or vortex point (0,0) (b) Between the origin and the complementary saddle point ( Omax. 0). Therefore, Eq.(7.19) defines a useful Liapunov stability bound and an operating limit, which produces the largest possible stability domain. It should be emphasized that the application of the Liapunov method depends strongly on the choice of the Liapunov function. However since the Energy Metric Algorithm is unique in power system stability studies for generating a Liapunov function directly from the system equations, the search is made somewhat easier and if desired this search could be performed with the aid of a digital computer. Figure 7.3 illustrates the effect of the equivalent mechanical power output of the shaft for the prefault system loading (Pm12) on the phase plane postfault system stability domain. Stability domains corresponding to the smaller the value of Pm12, the greater will be the region of phase plane enclosed by the stability boundary.

PAGE 66

N o NO X o (Xl UI -. Pm 1 2 = -0.18 p.u. Pm 1 2 = -0.2 p.u. Pm 12 = -0.12 p.u._ ....... Pm 12 = -0.16p.u. > 1-250.00 -200.00 -150.00 -100.00 State Variable -50.00 Xl. I n 0.00 degree 50.00 Figure 7.3 Effect of the Equivalent Mechanical Power (Pm12) on the Phase Plane Postfault Stability Boundary 100.00

PAGE 67

( CHAPl'ER VI I MANSOUR'S MODIFICATION OF THE EMA 8.1 Introduction The direct method of Liapunov has been tried in many dif-ferent types of applications with varying degrees of success. In particular, the Energy Metric Algorithm was applied by Mansour [10, 11] to power systems with great success. In analyzing the transient stability of multi-machine power system, Mansour employed a number of different methods in generating Liapunov functions. From a comparison of the results produced by each of the different methods Mansour concluded that the Energy-MetricAlgorithm with a slight modification produced the largest region of s tabil i ty 8.2 Power System Representation The equivalent two-finite machine power system model of Figure 3.1 will be used in the analysis for finding a larger stability boundary. The system equations of motion are M d2S + D...EL = Pm(S) Pe (F4,S) dt2 dt (8.1) Eq= fl(F4,S, Eex) (8.2) Pm= f2( Pm, S ) (8.3)

PAGE 68

Where 0= Torque angle M= Inertia constant D= Damping coefficient Pm= Mechanical power output of the shaft Pe= Electrical input power Eex= Exciter voltage EQ=Voltage proportional to field flux linkage EQ=Flux decay Pm= Governor action 0= Angular speed deviation 8.3 System Equation of Motion 58 The power system is simplified for convenience by assuming constant (independent of 0) uniform ( DIM ratio is the same for each machine) damping, no transient saliency, no flux decay and negligible governor action. This will yield the postfault system model described by the state Eq.(6.17) and (6.18), with constant uniform damping introduced into the system model. A more detailed study of the power system transient stability problem was recently accomplished by Roemish[5]. The system equations become: Where d28 M--+ dt2 Pm12 Pe12 dt Pe12= Pc + Pmax. sin (8ss + (8.4) (8.5)

PAGE 69

59 The system equations in state-variable form are Xl = S 4> and 4> = Sss + "f X2 dS --dt Let Xl = X2 (8.6) X2M= -DX2+ [Pml2 (Pc + Pmax. sin (Xl+ 4)] (8.7) 8.4 The unmodified EMA The Energy Metric Algorithm is applied to the system under study to determine a suitable Liapunov candidate function (V-function) dXl X2 dX2 =M D X2 + [Pml2 (Pc + Pmax. sine Xl + 4> ))] (8.8) Since n = 2, there is only one equation formed. After the cross multiplication and rearrangement of the terms, there results: m = D X2 dXl+ [ Pc + Pmax. sin ( Xl + 4> ) -Pml2 ] dXl + M X2 dX2 (8.9) Next, a line integration is performed to obtain the V-function. V(X)= J 111 J Xl(X2 = 0) V(X)= D X2 dXl +[ o J x2 (Xl=x l ) + M X2 dX2 o Pc + Pmax.sin( Xl + 4> ) -Pml2 ] dXl (8 10)

PAGE 70

60 V(X)= [ PC(Xl) + Pmax.[cos COS(XI + Pm12(Xl)] 2 + 0.5 M X2 Quadratic Term Periodic Terms (B.11) The Liapunov function [VeX)] obtained by applying the EMA to the system equations of a two-finite machine power system produces an incomplete quadratic form given by Eq.(B.11). By completing the square of the quadratic form, the Liapunov function candidate results in a more accuratelly defined stability region for the power system under investigation than can be obtained by assuming a totally conservative system. B.5 Mansour's Modified EMA An improvement in the stability domain within which a power system can operate means a better response more accuracy. This means that under load and fault conditions the system can be maintained in synchronism and thus avoid the need for emergency load dropping. The missing terms from the candidate function can be supplied by examining Eq.(B.11). Considering only the quadratic terms [denoted as Vq(X)] of the Liapunov function presented in the equation, which means ignoring the periodic terms in brackets for simplicity. Mansour proposed adding the terms axi and p XIX2 to the Liapunov function originally determined from

PAGE 71

61 the system equation. The resulting generalization of Eq.(8.11) is 2 2 V(X)= 0.5M X2 + P X1X2 + a Xl + [[ Pc Pm12](Xl) + Pmax.[ cos cos ( Xl + )]] (8.12) 2 2 Where Vq(X)= 0.5 M X2 + P X1X2 + a Xl From Liapunov's theorem it is necessary that the V-function to be positive definite for stability. The values of a and P for which the quadratic form is positive definite can be determined by the Hessian matrix of second partials. a 2 v a 2 v ad aX2Xl *= a 2 v a 2 v aX1X2 ad A real and symmetric matrix can be formed from the coefficients of the quadratic form of T Vq(X)= X X (8.13) a 2 v a 2 v Xl ad aX2Xl Vq(X)= [Xl X2 ] a 2 v a 2 v (8.14) aX1X2 X2 For the function V to be positive definite the following require-ments must be satisfied 2 a V > 0 ad (8.15) det > 0 (8.16)

PAGE 72

Evaluation of the Hessian matrix yields: av :::} a 2 v --= 2 aXl --= 2 a ad axl av + P Xl :::} a 2 v --= M X2 --= M aX2 ad a 2 v a 2 v = P --axlX2 ax2XI = [ 2 a P (J M Since det > 0 a > 0 From this there results and Since det = 2 a M p2 2 a>...L -2M 2 2 Vq(X)= a Xl + (J XIX2 + 0.5 M X2 The total derivative of Vq(X) with respect av aVg Vq(X)= Xl + X2 axl aX2 to time Vq(X) = ( 2a Xl + (J X2 ) Xl + ( M X2 + (J Xl is )X2 Since Xl = X2 and from Eq. (8.7) D X2 = M X2 Vq(X) =(2a Xl + {J X2)X2 + (M X2 + P XI)()X2 Vq(X) = ( (J -D ) + ( 2a )XI X2 a is selected so that V is negative definite. This can be 62 (8 17) (8.18) (8.19) (8.20) (8.21) (8.22) (8.23) (8.24) (8.25) (8.26) (8.27) (8.28)

PAGE 73

63 accomplished if ID2 a = 2M (8.29) Therefore 2 Vq(X) = ( p -D ) X2 (8.30) A resulting Liapunov function is then obtained by the sum of the periodic and quadratic terms to give V t t 1= V d t' + V d' o a qua ra IC perlO IC 2 2 V(X)= 0.5 M X2 + P XIX2 + a Xl + [[ Pc PmI2](XI) + Pmax. [ cos cos ( Xl + )]] (8.31) Since a = ,from Eq. (8.29) 2 2 V(X)= 0.5 M X2 + P XIX2 + Xl + [[ Pc Pm12](XI) + Pmax. [ cos cos ( Xl + )]] (8.32) The time derivative of Vtotal(X) is given by Since + PX2 + [Pc -Pm12 + Pmax. sin(XI + 2M + [ P Xl + M X2 ]X2 (8.33) Xl = X2 D and X2 = X2 Xl + PX2 + [Pc -Pml2 + Pmax. sin(XI + 2M + [ P Xl + M X2 ] X2 (8.34) 2 V(X)=(P -D)X2 + [Pc -Pml2 + Pmax. sin(XI + (8.35) By varying the parameter P the largest stability region within a reasonable degree of accuracy can be determined.

PAGE 74

CHAPTER IX CONCLUSIONS The resulting knowledge gained from a two-finite machine system which has been studied in this thesis, allows a number of general conclusions to be drawn concerning the effect on stability of certain concepts used in power system design, apparatus design, and power system operation. The effect of system modifications must be analytically observed before a fault, during a fault, and after fault clearance. Experience has shown that some design changes improve stability during all three conditions, while other modifications are helpful during one condition and detrimental during other changes. Out of the three methods studied in this thesis, it is found that the second method of the Liapunov is exact for two machines, conservative for three or more machines. This is one advantage of this method over the other two methodsC equalarea criteria and phase-plane trajectory). These two methods, however, are suitable for a two-machine system. The methods developed and used in this thesis provide an on-line digital computer approach to the power system stability problem, which is functional and convenient. Because of the

PAGE 75

65 complexity of the problem and the constant advancement in the state-of-the-art, there is still much further research to be considered in this area. For instance, for complete automatic operation and security maintenance, greater precision in the development of the largest stability region for a given system must be considered. Such research requires the development of the necessary as well as the sufficient conditions for specific Liapunov functions.

PAGE 76

BIBLIOGRAPHY 1. Anderson, P.M., and A.A. Fouad, Power System Control and Stability, Vol.l. The Iowa State University Press, Ames. Iowa, 1977. 2. Stevenson, D., Elements of Power System Analysis, 4th Edition McGraw-Hill, New York, 1982. 3. Gless, G.E., "The Direct Method Of Liapunov Applied To Transient Power System". IEEE Transactions, PAS-85, 1966. 4. Wall, E.T., "A Topological Approach To The Generation Of Liapunov Functions." Acta Technica Csav, Prague, Czachoslovakia, No.2, 1968. 5. Roemish, W.R., A New Synchronous Generator Out-0f-Step Relay Scheme. Ph.D. Dissertation, University Of Colorado At Boulder, 1981. 6. Meyer, J.e., Modern Control Principles And Application, McGraw-Hill New York, 1968. 7. Jordan, D.W., and P.Smith, Non Linear Ordinary Differential Equations 2nd edition. Oxfored University Press, New York, 1987. 8. Ogata, K., State-Space Analysis Of Control Systems, Prentice-Hall, Englewood Cliffs, N.J. 1967.

PAGE 77

67 9. Siddiqee, M.W., "Transient Stability of an A-C Generator by Liapunov's Direct method." Int. J. Control Vol. 8, No.2 1968, PP.131-144. 10. Mansour, M., "Stability Analysis and Control of Power System." in' Real Time Control of Electric Power System, Edmund Hardschin, Ed. Elsevier, New York 1972, PP. 11. Mansour, M., "Generalized Liapunov Functional for Power Systems." IEEE Transactions on Automatic Control, Vol. AC-19 PP. 247-248, 1974. 13. Kimbark, E.W., Power System Stability, Vol.l. Wiley, New York, 1956. 14. Lewis, W.A., A Basic Analysis of Synchronous Machines Pt. 1. AlEE Transaction, PAS-77, PP. 436-55, 1958. 15. Adkins, Bernard, General Theory of Electrical Machines. Chapman and Hall, Ltd. London, 1959. 16. Crary, B Selden, Power System Stability, Vol.l. John Wiley and sons, Inc. New York, 1945. 17. Wall, E. T., The Generation of Liapunov Function in Automatic Control Theory by an Energy Metric Algorithm. Ph.D. Dissertation University of Denver, Denver, Colorado, 1967. 18. Weiss, J.R., "Transient Asymptotic Stability of Power Systems as Established with Liapunov Functions." IEEE Transactions, 1976, PP.1480-86

PAGE 78

APPENDIX A PROGRAM AND FLOWCHART This appendix contains the main program, and the listing of computed swing curve coordinates. A.I Main Program This main program "DELTA", written in turbo Pascal, the program is designed to implement the modified-Euler procedure for the two generator, three-bus problem. A 2 Description of the Input Quantities Name in the program Pal, Pa2 OMEGAI, OMEGA2 DELTA I DELTA2 Symbol in the analysis Pal, Pa2 WI, W2 Description The accelerating output power in p.U. The machine shaft speeds in rad./sec. The machine torque angles in degrees. The relative torque angle between the two machine in degrees.

PAGE 79

T dt C1, C2 t at 11Ml, 11M2 69 Time. The time interval Inertia constant of the machine in rad./sec2 The in'itial estimation of shaft speed of machine 1. The final estimation of shaft speed of machine 1. The initial estimation of snaft speed of machine 2. The final estimation of shaft speed of machine 2. The initial estimation of torque angle of machine 1. The final estimation of torque angle of machine 1. The initial estimation of torque angle of machine 2. The final estimation of torque angle of machine 2.

PAGE 80

70 A.3 Description of the Output Quantities Name in the program Pal, Pa2 OMEGAl,OMEGA2 DELTAl, DELTA2 TIME Symbol in the analysis Pal, Pa2 WI, W2 t Description The accelerating output power in p.u. The machine shaft speeds in rad./sec. The machine torque angles in degrees. The relative torque angle between the two machine in degrees. Time.

PAGE 81

71 PROGRAM DELTA; TYPE = ARRAY [1 .. 250] OF REAL; VAR I: INTEGER; dt,Dl,D2,Cl,C2:REAL; dOMEGAl_START,dOMEGA2_START,dOMEGAl_END,dOMEGA2_END:REAL; dDl_START,dD2_START,dDl_END,dD2_END:REAL; OUTPUT_FILE: TEXT; TIME,DELTAl,DELTA2,DELTA_DIFF:LIST_250; Pal,Pa2,OMEGAl,OMEGA2:LIST_250; PROCEDURE INITIALIZE; BEGIN END; dt:=0.05; Cl:=75.3984; C2:=39.27; FOR I : = 1 TO 250 00 TlME[!] := dt (I-I); DELTAl[1]:=10.0/57.2985; DELTA2[1]:=12.65/57.2985; DELTA_DIFF[1]:=DELTAl[1]-DELTA2[1]; OMEGAl[1]:=377; OMEGA2[ 1] : =377; PROCEDURE COMPUTE_DELTA; BEGIN FOR I : = 1 TO 11 0 00 BEGIN IF (1<=2) THEN BEGIN Pnl[I]:=0.2999-0.06*SIN(DELTAl[I]-DELTA2[1]+0.01745); Pa2L 1:1: =0.992+0. 06*SIN(DEL'l'All I I 1-0. 017t1G); END ELSE BEGIN Pal[I]:=0.2024-1.15*SIN(DELTAl[I]-DELTA2[I]+0.22689); Pa2[1):=0.314+1 15*SIN(DELTAl[I)-DELTA2[1)-0.22689); END; dOMEGAl_START:=Cl*Pal[I]; dOMEGA2_START:=C2*Pa2[1]; IF 1=1 THEN BEGIN dDl_START: =0.0; dD2_START:=0.0 END; IF 1>1 THEN BEGIN dDl_START:=dDl_END;

PAGE 82

dD2_START:=dD2_END; END; Dl:=DELTAl[I]+dDl_START*dt; D2:=DELTA2[1]+dD2_START*dt; IF (1<=2) THEN BEGIN Pal[l+l]:=O.2999-0.06*SIN(Dl-D2+0.01745); Pa2[I+l]:=O.992+0.06*SIN(Dl-D2-0.01745); END ELSE BEGIN Pal[l+l]:=O.2024-1.15*SIN(Dl-D2+0.22689); Pa2[I+l]:=O.314+1.15*SIN(Dl-D2-0.22689); END; dOMEGAl_END:=Cl*Pal[l+l]; dOMEGA2_END:=C2*Pa2[I+I]; 72 OMEGAI[I+l]:=OMEGAl[I]+(dOMEGAl_START+dOMEGAl_END)/2*dt; OMEGA2[I+l]:=O!ffiGA2[1]+(dOMEGA2_START+dOMEGA2_END)/2*dt; dDI_END:=OMEGAI[I+1]-377; dD2_END:=OMEGA2[I+l]-377; DELTAl[I+l]:=DELTAl[I]+(dDl_END+Ddl_START)/2*dt; DELTA_DIFF[I+l]:=DELTA1[I+1]-DELTA2[I+1]; END; END; PROCEDURE WRITE_RESULTS; BEGIN ASSIGN(OUTPUT_FILE, 'A:D_OUTPUT'); REWRITE(OUTPUT_FILE); WRITELN(OUTPUT_FILE,'TIME Pal Pa2 OMEGAl DELTAl OMEGA2 DELTA 2 WRITELN(OUTPUT_FILE); FOR I : = 1 TO 41 00 WRITELN(OUTPUT_FILE,TIME[I]:lO:2,Pal[I]:9:3,Pa2[I]:10:3, DELTAl[I]*57.2985:10:2,OMEGA2[I]:11:2, DELTA2[I]*57.2985:9:2,DELTA_DIFF[I]*57.2:10:2); CLOSE(OUTPUT_FILE); END; PROCEDURE MAIN; BEGIN INITIALIZE; COMPUTE_DELTA; WRITE_RESULTS END; BEGIN MAIN END.

PAGE 83

73 LISTING OF COMPUTED SWING CURVE COORDINATES TIME Pal P a 2 OMEGA 1 DELTA 1 OMEGA2 DELTA2 DELTAl-2 0 00 0 302 0 988 377 00 10. 00 377.00 12 .65 -2. 6 5 0.05 0.303 0 987 378.14 11.63 378 .94 15 .43 -3.79 0 10 0 087 -0. 08 4 379.28 16.53 380.88 23 .76 -7.22 0 15 0.157 -0. 1 4 9 379 .78 23.79 380.63 34 .51 -10.71 0.20 0 176 -0. 166 380 .47 32 .74 380 .29 44.43 -11 .67 0.25 0 139 -0.132 381.11 43.59 379 .98 53.41 -9.80 0 30 0 056 -0. 055 381. 51 55.94 379.78 61. 64 -5.70 0 .35 -0. 043 0 042 381.54 68.90 379 .76 69 58 -0.68 0 4 0 -0. 127 0.127 381.19 81.40 379 .94 77 75 3 6 4 0 .45 -0.170 0 172 380.58 92.53 380.26 86 .64 5.88 0.50 -(J .lt;9 0.161 379.91 101.83 380.62 96 .GO G .3:1 0.55 -0.098 0 098 379 38 109.41 380 .89 107 .26 2.15 0 .60 -0.005 0.004 379 .17 115 .93 381. 00 118 .57 -2.63 0 .65 0.092 -0. 089 379 .35 122.41 380.91 129 .91 -7.48 0 .70 0.160 -0.151 379.87 129 .88 380 .66 140 .75 -10. 84 0 .75 0 176 -0. 166 380.55 13-9.08 380 .32 150 .74 -11.64 0 80 0 135 -0.128 381 .19 150 .17 3 80 .01 159 80 -9.61 0 .85 0 051 -0. 049 381.57 162 .72 379 .82 168 14 -5.41 0 .90 -0. 049 0 048 381.57 175.82 379.81 176.21 -0.39 0 .95 -0. 131 0 131 381 .21 188 .40 380 .01 184 54 3 85 1.00 -0. 171 0 173 380 .59 199 56 380.33 193 .62 5 .93 1.05 0 158 379 .92 208 .89 380 ,68 203 .67 5 .21 1.10 -0. 093 0.093 379.41 216 .51 380 .95 214 .61 1. 90 1.15 0 001 -0.002 379.22 223 .13 381. 05 226 .07 -2.93 1.20 0 097 -0. 093 379.42 229.77 380 .95 : 237 52 7 .74 1.25 0 162 -0. 154 379.95 237.47 380 .68 248 .46 -10.97 1. 30 0 175 -0. 165 380.64 246.91 380 .35 258 .53 -11. 6 0 1. 35 0 131 -0. 125 381 .27 258 .24 380.04 267 .67 -9.4 2 1.40 0 '.04 5 -0. 044 381. 63 270.98 379 .86 276 .11 -5. 12 1.45 -0. 05 4 0 053 381. 61 284 22 379.87 284 .31 -0. 10 1.50 -0. 135 0 135 381.22 296 .87 380.07 292.82 4.04 1. 55 -0.172 0.174 380.60 308 .07 380 4 0 302 .08 5 9 8 1. 60 -0. 155 0 156 379 .93 317.41 380 .75 312 32 5 .08 1. 6 5 0 088 0 088 379 .43 325.09 381.01 323.43 1.65 1 70. 0 007 -0.008 379.26 331. 81 381.10 335.05 -3.23 1. 75 0 102 -0.098 379.49 338.62 380.98 346 .62 -7.99 1. 80 0.165 -0.156 380.04 346 .53 380 .71 357 .64 -11.09 1.85 0.174 -0.164 380.73 356 .22 380.37 367 .79 -11.5 5 1. 90 0 127 -0. 121 381. 35 367 79 380 .07 377 .02 -9.21 1.95 0 039 -0. 038 381.69 380.72 379.90 385 56 -4.83 2 .00 -0. 0 6 0 0.059 381.64 394 .09 379.92 393 .90 0 1 9

PAGE 84

A.4 Flowchart INITIALIZE dt. C1. C2, TIME DELTA1, DELTA2, OMEG1, OMEGA2, DELTA_DIFF 74 CLOSE OlITPlIT FILE

PAGE 85

Flowchart Continuous COMPUTE Pa1, Pa2 (during f.ault) 75 COMPUTE dOMEGA1_END, dOMEGA2_END,OMEGA1, OMEGA2,dD1_END.dD2_END DELTA1. DELTA2. DELTA-DIFF

PAGE 86

APPENDIX B LISTING OF OOMPUTED PHASE-PLANE TRAJECTORY OOORDINATES The main program used to solve the non-linear differential equations called Simnon. Simnon is a program that can solve dif-ferential, and difference equations using a fourth order Runge-Kutta integration algorithm, it can also be used to simulate dynamical systems that are composed of subsystems. The solution to the problem can be divided in the following steps. 1. Enter system descriptions 2. Simulate 3. Analyze the results 4. Change parameters and initial condition. This appendix contains the system descriptions program 'list, subroutines (Macros) which perform a whole sequence of commands, and listing of computed phase-plane trajectory coordinates.

PAGE 87

Simnon on Interactive Simulation Program for Nonlinear Systems. Version 1.00 Copyright (c) Deportmcnt of Automatic Control Lund InAtitute of Technology. Lund. SWEDF:N 19a6 All R !ghts Rcscrved ---a61020--011371---lXlN'I'INtJOtJS SYSTEM IIlJSIH STATE Xl X2 DEn. DXl DX2 TIME T IlXl=x2 DX2=AIl+SIN(xl+PI/J80+J))+JaO/I'I) A= I Jo' 1 '
PAGE 88

5.8E-2 6.E-2 6.2E-2 6.4E-2 6.6E-2 6.8E-2 7 ; ,E-2 7 .2E-2 7 7.6E-2 7.8E-2 8.E-2 8.2E-2. 8.4E-2 0.68-2 0.OE-2 9.E-2 9.2E-2 9.4E-2 9.6E-2 9.8E-2 9.99999E-2 0.100027 0.100055 0.1000B2 0.1003B6 0.102 0.104 0.106 0.108 0.11 0.112 0 .114 0.116 0 0.12 0.122 0.124 0.126 0.12B 0.13 0.132 0.134 0.136 0.13B 0 0.142 0.144 0.146 0.148 0.15 0.152 0 .154 0.156 0.158 0.16 0.162 0.164 0.166 0.16B 0.17 0.172 0 .174 0.176 b .17B 0.1 B 0.182 0.lB4 0..1 B 6 O.lBB 0.19 0.192 -53.0963 -56.7523 -60.52B9 -64.42B4 -68.4534 -72.6072 -76.8933 -81.3159 -B5.8799 -90.5904 -95.4534 -100.475 -105.663 -111.024 -116.567 '-122.3 374 -140.73. 4 -.147.322 -154.149 -161.224 -161.323 -161.421 -161.518 -162.5B2 -167.57 -172.188 -175.082 -176.255 -175.713 -173.454 -169.47l -163.76 -156.33 -147.21B -136.503 -124.331 -110.927 -96.5966 -81.7202 -66.724 2 -52.0439 -38.082 -25.1769 -13.585 -3.4791 5.04026 11.9273 17.172 20.7829 22.7741 23.1564 21.9319 19.0936 14.6283 8.52321 0.7777l3 -8.58175 -19.4831 -31.7898 -45.2861 -59.67l3 -74.5703 -89.5583 -104.2 -118.088 -130.879 -142.31 .-152.199 -160.437 -166.968 -1798.14 -1857.9B -191B.B3 -19BO.B9 -2044.41 -2109.64 -2176.82 -2246.24 -231B.16 -2392.87 -2470.66 -2551. 81 -2636. 61 -2725.34 -2B1B.26 -2915.6 -3017.59 -3124.41 -3236.18 -3352.98 -3474.8 -3601.56 -3590.92 -3579.22 -3567.51 -3437.23 -2741.42 -1877.52 -1016.21 -157.811 700.042 1560.13 2423.38 3286.59 4140.15 4966.03 5736.75 6416.12 6962.64 7335.88 7504.68 7454.62 7191.87 6741.65 6142.01 5435.58 4662.08 3853.62 3032.73 2212.75 1399.39 592.787 -210.352 -1014.76 .-1824.66 -2641.81 -3463.42 -4280.01 -5073.58 -5816.5 -6472.26 -6998.82 -7354.92 -7508.26 -7443.07 -7l64.14 -6695.42 -6074. -5341.8 -4538.16 -3694.95 -2834.66 78

PAGE 89

Simnon -nn Interactive Simulation Program for Nonlinear Systems. Version 1.00 Copyright (c) Department of Automatic Control Lund Institute of Technology, Lund, SWEDEN 1986 All Rights Reserved ---861020--011371---CX)NTlNUOUS SYSTEM I1USD2 HTA'I'I'; X I X2 DX I DX2 TIME T DX1=x2 DX2=A-BSIN{x1PI/180+D/PI) A=IF T
PAGE 90

80 Print of Simnon store file:STORE For the clearing time t 0.07 sec. TIME Xl X2 O. 0.7 O. 2.E-3 0.632745 -67.2498 4.E-3 0.431045 -134.434 6.E-3 9.50969E-2 -201.467 8.E-3 -0.374773 -266.345 I:E-2 -0.978113 -334.947 1.2E-2 -1.71435 -401.233 1.4E-2 -2.5828 -467.148 1.6E-2 -3.58266 -532.639 1. 6E-2 -4.71304 -597.66 2.E-2 -5.97296 -662.17 2.2E-2 -7.36135 -726.132 2.4E-2 -8.8771 -769.518 2.6E-2 -10.519 -852.308 2.6E-2 -12.2859 -914.488 3.E-2 -14.1766 -976.055 3.2E-2 -16.1897 -1037.01 3.48-2 -18.3242 -1097.38 ) 6r-:-2 -20.5799 -1157.19 1.111-:":' : () I, :-; I ; -I? I (I. I ri 4.1::2 -25.11444 -12'15.2'1 4.2E-2 -28.0534 -1333.66 4.4E-2 -30.7788 -1391.71 1.6E-2 -33.6201 -1449.53 4.8E-2 -36.5769 -1507.2 5.E-2 -39.6489 -1564.85 5.2E-2 -42.8363 -1622.62 5.4E-2 -46.1395 -1680.64 5.68-2 -49.5592 -1739.09 5.6E-2 -53.0963 -1798.14 6.E-2 -56.7523 -1857.98 6.2E-2 -60.5289 -1918.83 6.4E-2 -64.4264 -1960.69 6.6E-2 -66.4534 -2044.41 6.6E-2 -72.6072 -2109.64 7.E-2 -76.8933 -2176.82 7.2E-2 -81.2454 -2160.13 7.4 E-2 -85.4915 -2075.63 7.6E-2 -89.5074 -1930.59 7.8E-2 -93.1766 -1729.81 8.E-2 -96.3938 -1479.76 8.2E-2 -99.068 -1188.22 8.4E-2 -101.125 -863.856 8.6E-2 -102.506 -515.958 6.88-2 -103.179 -154.14 9.E-2 -103.122 211.815 9.2E-2 -102.336 572.062 9.4E-2 -100.844 916.696 9.6E-2 -98.685 1236.66 9.6E-2 -95.9196 1522.23 9 .99999E-2 -92.6246 1764.98 0.102 -88.8933 1957.35 0.104 -84.833 2093.1 0.106 -80.5617 2167.73 0.106 -76.2046 2178.76 0.11 -71. 8894 2125.94 0.112 -67.7422 2011.2 0.114 -63.8832 1838.58 0.116 -60.4226 1613.85 0.118 -57.4578 1344.2 0.12 -55.0703 1037.82 0.122 -53.3251 703.567 0.124 -52.2685 350.677 0.126 -51.9285 -11.4194 0.128 -52.314 -373.219 0.13 -53.4148 -725.215 0.132 -55.202 -1057.99 0.134 -57.6278 -1362.32 0.136 -60.6264 -1629.38 0.138 -64.1151 -1851.03 0 .14 -67.9955 -2020.17 0.142 -72.1569 -2131.11 0.144 -76.4786 -2179.96 0.146 -80.8311 -2164.89 0.148 -85.0958 -2086.33 0 .15 -89.1388 -1946.87 0.152 -92.8157 -1751.13

PAGE 91

81 .. Print of Simnon store file: STORE For the clearing time t 0.08 sec. TIME Xl X2 O. 0.7 O. 2.E-3 4.E-3 0.431045 -134.134 6.E-3 9.50969E-2 -201.487 8.E-3 -0.374773 -268.345 1.E-2 -0.978113 -334.947 1.2F.:-2 -1.71435 -401.233 1.11':-2 -2.5828 -467.118 1.61':-2 -3.56266 -532.639 1.8E-2 -4.71304 -597.66 2.E-2 -5.97296 -662.17 2.2E-2 -7.36135 -726.132 2.4E-2 -8.8771 -789.518 2.6E-2 -10.519 -852.308 2.0E-2 -12.2859 -914.468 3.E-2 -14.1766 -976.055 3.2E-2 -16.1897 -1037.01 3.4E-2 -18.3242 -1097.38 3.6E-2 -20.5789 -1157.19 3.0-2 -22.9526 -1216.46 4.-2 -25.4414 -1275.27 4.2E2 -28.0534 -1333.66 4.4-2 -30.7788 -1391.71 4.6E-2 -33.6201 -1449.53 4.8E-2 -36.5769 -1507.2 5.E-2 -39.6489 -1564.85 5.2-2 -42.8363 -1622.62 5.4E-2 -46.1395 -1680.64 5.6E-2 -49.5592 -1739.09 5.8E-2 -53.0963 -1798.14 6.E-2 -56.7523 -1857.98 6.2E-2 -60.5289 -1918.83 6.4E-2 -'64.4284 -1980.89 6.6E-2 -68.4534 -2044.41 6.8E-2 -72.6072 -2109.64 7.E-2 -76.8933 -2176.82 7.2E-2 -81.3159 -2246.24 7 .4E-2 -85.B799 -2318.16 7.6E-2 -90.5904 -2392.87 7.8E-2 -95.4534 -2470.66 8.E-2 -100.475 -2551.B1 B.2E-2 -105.291 -2221.5B 8.39476E-2 -109.228 -lB12.4 8 .4E-2 -109.322 -1800.73 8.43141E-2 -109.877 -1730.12 B.6E-2 -112.464 -1334.1 8.81':-2 -114.637 -B34.679 9.E-2 -115.789 -315.251 9.2E-2 -115.893 211.624 9.4E-2 -114.946 733.507 9.6E-2 -112.971 1237.88 9.8E-2 -110.015 1711.99 9.99999E-2 -106.152. 2142.89 0.102 -101.481 2517.68 0.104 -96.1265 2824.08 0.106 -90.2372 3051.22 0.108 -B3.9803 3190.52 0.11 -77.5374 3236.50 0.112 -71.0973 3187.75 0.114 -64.8482 3046.31 0.116 -5B.9699 2B1B.22 0.118 -53.6272 2512.4 0.12 -48.964B 2139.91 0.122 -45.1039 1712.96 0.124 -42.1409 1244.19 0.126 -40.1467 746.109 0 .128 -39.167B 230.917 0.13 -39.2265 -289.499 0.132 -40.3215 -B03.36 0.134. -42.4277 -1296.75 0.136 -45.496 -1763.44 O.13B -49.4526 -2184.88 0.14 -54.19B3 -2550.45 0.142 -59.6091 -2848.02 0 .144 -65.5378 -3066.74 0.146 -71. 6177 -3197.97 0.146 -78.2676 -3236.14

PAGE 92

82 Print of Simnon store file: STORE For the clearing time t 0.09 sec. TIME Xl X2 O. 0.7 O. 2.E-3 0.632745 -67.2498 4.E-3 0.431045 -134.434 6.E-3 9.50969E-2 -201.487 8.E-3 -0.374773 -268.345 1. E-2 -0.978113 -334.947 1.2E-2 -1.71435 -401.233 1.48,.2 -2.5828 -467.148 1.6E-2 -3.58266 -532.639 1.8E-2 -4.71304 -597.66 2.E-2 -5.97296 -662.17 2.2E-2" -7.36135 -726.132 2.4E-2 -8.8771 -789.518 2.68-2 -10.519 -852.308 2.88-2 -]2.2859 -914.488 3.E-2 -14.1766 -976.055 3.2E-2 -16.1897 -1037.01 3.4E-2 -18.3242 -1097.38 3.6E-2 -20.5789 -1157.19 3.88-2 -22.9526 -1216.46 4.E-2 -25.4441\ -1275.27 4.2E-2 -28.0534 -1333.66 4.4E-2 -30.7788 -1391.71 4.6E-2 -33.6201 -1449.53 4.813-2 -36.5769 -1507.2 5.8-2 -39.6489 -1564.85 5.2E-2 -42.8363 -1622.62 5.48-2 -46.1395 -1680.64 5.68-2 -49.5592 -1739.09 5.8E-2 -53.0963 -1798. 14 6.-2 -56.7523 -1857.98 6.2-2 -60.5289 -1918.83 6.4E-2 -64.4284 -1980.89 6.613-2 -68.4534 -2041.41 6.8E-2 -72.6072 -2109.64 7.-2 -76.8933 -2176.82 7.2-2 -81.3159 -2246.24 7.413-2 -85.8799 -2318.16 7.6E-2 -90.5904 -2392.87 7.8E-2 -95.4534 -2470.66 8.-2 -100.475 -2551.81 8.2E-2 -105.663 -2636.61 8.4E-2 -111.024 -2725.34 8.6E-2 -116.567 -2818.26 8.8E-2 -122.3 -2915.6 9.-2 -128.233 -3017.59 9.00277E-2 -128.316 -3009.26 9.00555E-2 -128.4 -3000.08 9.00833E-2 -128.483 -2990.9 9.03816E-2 -129.361 -2891.62 9.28-2 -133.593 -2333.53 9.4E-2 -137.54 -1608.52 9.6E-2 -140.01 -858.528 9.88-2 -}40.966 -96.1357 9.99999E-2 -HO.394 667.604 0.102 -138.302 1421.81 0.104 -134.721 2154.42 0.106 -129.708 2850. 98 0.108 -123.353 3494.09 0.11 -115.781 4063.41 0.112 -107.163 4536.88 0 .114 -97.712 4892. 84 0.116 -87.6824 5112.87 0.118 -77.3594 5184.84 0.12 -67.0444 5105.07 0 .122 -57.0369 4879.03 0.124 -47.6169 4520.37 0.126 -39.031 4048.42 0.120 -31. 484 3485.23 0.13 -25.1363 2852.65 0.132 -20.1065 2170.37 0 .134 -16.4769 1454.87 0.136 -14.3002 719.449 0.138 -13.6051 -25.1134 0.11\ -14.1\004 -769.364 0.142 -16.6761 -1503.81

PAGE 93

" E'rint of Silllnon store file:STOH TIME O. 2.E-3 .4.-3 6.E-3 0.E-3 1.-2 1.2E-:2 1.1-2 1.6E-2 1.8E-2 2.E-2 2.2E-2 2.41>-2 2.6E-2 2.6-2 3.E"'"2 3.2-2 3.4-2 3.6E-2 3.08-2 4.E-2 4.2E-2 4.4E-2 4.6-2 4.0E-2 .5.E-2 5.28-2 -'5.11:)-2 5.6-2 5.8Ei-2 6.-2 6.2E-2 6.48-2 6.65-2 6.0E-2 7.E-2 7.28-2 7.18-2 7.68-2 7.0E-2 0.E-2 8.2E-2 0.11,-2 0.682 0 .OE2 9.E-2 9.28-2 9.4E-2 9.68-2 9.08-2 9.99999-2 0;100027 0.100055 0.100062 0_100386 0_102 0 .104 0.106 0.106 0 .11 0.112 0.114 0.116 0.110 0.12 0.122 0.124 0.126 0.128 0.13 0.132 0.134 0.136 0.130 Xl 0.7 0.632745 0 .431045 9.509698-2 -0.374773 -0.970113 -1.71435 -2.5028 -3.58266 -4 71304 -5.97296 -7.36135 -0.0"171 -10.519 -12.2059 -11.1766 -16.1097 -10.3242 -20.5789 -22.9526 -25.4444 -20.0534 -30.7700 -33.6201 -36.5769 -39.6489 -12.8363 -16.1395 -49.5592 -53.0963 -56.7523 -60.5209 -60.4534 -72.6072 -76.8933 -01.3159 -05.0799 -90.5904 -95.4534 -100.475 -105.663 -111_024 -116.567 -122.3 -120.233 -134.3"14 -140.734 -147.322' -154.149 -161.224 -161. 323 -161.421 -161.518 -162.582 -167.57 -172.188 -175.082 -176.255 -175 :713 -173.454 -169.471 -163.76 -156.33 -147.210 -136.503 -124.331 -110.927 -96.5966 -61.7202 -66.7242 -52.0439 -38.082 -25.1769 X2 O. -67.2496 -134.434 -201.487 -.268.345 -331.947 -401.233 -467.140 -532.639 -597.66 -662.17 -726.132 -789.518 -052.308 -914.460 -976.055 -1037.01 -1097.30 -1157.19 -1216.46 -1275.27 -1333.66 -1391.71 -1449.53 -1507.2 -1564.05 -1622.62 -1600.61 -1739.09 -1790.14 -1057.90 -1910.03 -1900.09 -2014.41 -2109.64 -2176.82 -2216.24 -2310.16 -2392 07 -2470.66 -2551. 01 -2636.6' 1 -2725.34 -2010.26 -2915.6 -3017.59 -3124.41 -3236.10 -3352.90 -3474.8 -3601.56 -3590_92 -3579.22 -3567.51 -3437.23 -2741.42 -1877.52 -1016.21 -157.811 700.042 1560.13 2423.38 3206.59 1140.15 4966.03 5736.75 6416.12 6962.64 7335. 88 7504.68 7454.62 7191. 87 6741.65 6142_01 83 For the c1earil1g time t O.! !:'ec.

PAGE 94

84 Print of Simnon store file:STORE For the clearing time t n 0.115 sec. TIME Xl x2 o. 0 7 o. 2.5E3 0.594918 -84.0545 5.E-3 0.279834 -167.981 7.5E-3 -0.244773 -251.652 1. E-2 -0.978112 -334.947 1.25E-2 -1.9191 -417.749 1.51'-2 -3.06635 -499.949 1.75E-2 -4.41826 -581.451 2 .E-2 -5.97296 -662.17 2.25E-2 -7.72839 -742.033 2 .5E-2 -9.68236 -820.989 2.75E-2 -11.8325 -899. 3.E-2 -14.1766 -976.055 3.25E-2 -16.712 -1052.16 3.5E-2 -19.4366 -1127.35 3.751'-2 .-22.3.481 -1201.69 4.E-2 -25.4444 -1275.27 4.25E-2 -28.7239 -1348.2 4.5E-2 -32.185 -1420.65 4.75E-2 -35.8260 -1192.79 5 .E-2 -39.6489 -1564.05 5.25E-2 -43.6513 -1637.09 5.5E-2 -47.8347 -1709.8 .5.75E-2 -52.2009 -1783.31 6.E-2 -56.7523 -1857.98 6.25E-2 -61.4922 -1934.22 6.5E-2 -66.425 -2012.46 6.751'-2 -71.5565 -2093.16 7.E-2 -76.8933 -2176.82 7.25E-2 -82.4435 -2263.97 7.5E-2 -88.2165 -:-2355.15 7.751'-2 -94.223 -2450.91 8.E-2 -100.475 -2551.81 8.251'-2 -106.987 -2658.42 8.51'-2 -113.773 -2771.26 8.75E-2 -120.849 -2890.84 9.1'-2 -128.233 -3017.6 9.25E-2 -135.943 -3151.89 9 .5E-2 -143.999 -3293.95 9.751'-2 -152.419 -3443.88 0.1 -161.224 -3601.57 0.1025 -170.433 -3766.69 0.105 -180.064 -3938.63 0.1075 -190.131 -4116.40 0 11 -200.65 -4298..97 0.1125 -211.629 -4484.47 0.115 -223.073 -4671.02 0.117371 -233 631 -4199.1 0.1175 -234.17 -4174.64 0.118788 -239.402 -3955.41 0.12 -244.091 -3790.73 0.1225 -253.264 -3575.3 0.125 -262.104 -3523.71 0.1275 -271.018 -3634.23 0.13 -200.412 -3908.45 0.1325 -290.699 -4350.02 0.135 -302.303 -4961.34 0.1375 -315.644 -5737.5 0.14 -331.111 -6656.5 0.1425 -349.003 -7666.26 0.145 -369.44 -8672.53 0.1475 -392.249 -9537.98 0.15 -416.881 -1.01058E4 0.1525 -442.423 -1.0251E4 0.155 -467.749 -9935.05 0.1575 -491. 768 -9224.51 0.16 -513.658 -8257.85 0.1625 -532.971 -7187.41 0.165 -549.611 -6135.47 0.1675 -563.729 -5179.17 0.17 -575.619 -4356.14 0.1725 -585.631 -3677.37 0.175 -594.123 -3139.24

PAGE 95

85 Print of Simnon store file:STORE For the clearing time t a 0.125 sec. TIME Xl X2 O. 0.7 O. 2.5E-3 0.594918 -84.0545 5.E-3 0.279834 -167.981 1.5E-3 -0.244773 -251.652 1. E-2 -0.978112 -334.947 1.25E-2 -1.9191 -417.749 1. 5E-2 -3.06635 -499.949 1.75E-2 -4.41826 -581.451 2.E-2 -5.97296 -662.17 2.25E-2 -7.72839 -742.033 2.5E-2 -9.68236 -820.989 2.75E-2 -11.8325 -899. 3.E-2 -14.1766 -976.055 3.251::-2 -16.712 -1052.16 3.58-2 -19.4366 -1127.35 3.758-2 -22.3481 -1201.69 4.E-2 -25.4444 -1275.27 4.25E-2 -28.7239 -1348.2 4.51::-2 -32.185 -1420.65 4.758-2 -35.8268 -1492.79 5.8-2 -39.6489. -1564.85 5.25E-2 -43.6513 -1637.09 5.58-2 -47.8347 -1709.8 5.75E-2 -52.2009 -1783.31 6.E-2 -56.7523 -1057.98 6.25E-2 -61.4922 -1934.22 6.5E-2 -66.425 -2012.46 6.75E-2 -71.5565 -2093.16 7.E-2 -76.8933 -2176.82 7.25E-2 -82.4435 -2263.97 7.5E-2 -88.2165 -2355.15 7.75E-2 -94.223 -2450.91 8.E-2 -100.475 -2551.81 8.25E-2 -106.987 -2658.42 8.5E-2 -113.773 -2771.26 8.75E-2 -120.849 -2890.84 9.E-2 -128.233 -3017.6 9.25E-2 -135.943 -3151. 89 9.5E-2 -143.999 -3293.95 9.75E-2 -152.419 -3443.88 0.1 -161.224 -3601.57 0.1025 -170.433 -3766.69 0.105 -180.064 -3938.63 0.1075 -190.131 -4116.48 0.11 -200.65 -4290.97 0.1125 -211.629 -4484.47 0.115 -223.073 -4671. 02 0.1175 -234.983 -4856.33 0.12 -247.352 -5037.9 0.1225 -260.167 -5213.14 0.125 -273.41 -5379.5 0.1275 -287.237 -5727.32 0.13 -302.239 -6311. 97 0.1325 -310.972 -7106.05 0.135 -337.904 -8060.29 0.1375 -359.331 -9081. 69 0,14 -383.249 -1. 00253E4 0.1425 -409.242 -1.07111E4 0.145 -436.452 -1. 09774E4 0.1475 -463.715 -1.07507E4 0.15 -489.836 -1.00833E4 0.1525 -513.888 -9124.79 0.155 -535.367 -8053.48 0.1575 -554.187 -7018.17 0.16 -570.57 -6115.91 0.1625 -584.92 -5397.5 0.165 -597.728 -4883.26 0.1675 -609.511 -4577.78 0.17 -620.79 -4479.89 0.1725 -632.082 --4588.06 0.175 -643.901 -4901. 76 0.1775 -656.761 -5419.53 0.18 -671.163 -6133.15 0.1825 -687.57 -7017.47

PAGE 96

Simnon -an Intcractivc Simulation Program for Nonlinear Systems. Version 1 .00 Copyright (c) Department of Automatic Control Lund Institute of Technology, Lund, SWEDEN 1986 All Rights Reserved ---861020--011371---CONTINUOUS SYSTEM HCP Xl X2 DER DX1 DX2 TIMgT DX1=x2 DX2=P-B+SIN(x1+PI/180+D+180/PI) D:1.38 B:7358 P :-151.2 PI: 3 .1415926 END MACHO IICPM SYST HCP STORE xl x2 INIT x1:7 INIT x2:0 SIMU 0 .2/Cl AXES H 250 75 V -11000 10000 SHOW x2(xl)/Cl STORE xl x2 INIT x1:-190 INIT x2:0 PAR P : 111. 69 SHIU 0 .2/C2 SIIOW x2(x1)/C2 STROE xl x2 INIT x1:50 INIT x2:0 PAR P:-19.7 SIMU 0 .25/C3 SIIOW x2(xll/C3 xl x2 INIT xl:70 INIT x2:0 PAR P :-151.2 SIMU 0 .25/C4 SIIOW x2(xll/C4 tex t 'PHASE PLANE PLOT OF X2 Vs. Xl mark a 15 1.2 mark "State Variable Xl In Degree mark a 1 2 15 mnrk v"State VAR X2 DIID/S END 86

PAGE 97

87 For Pm12 -0.16 p.u For Prn12 c -0.12 p.u TIME Xl X2 TIME Xl X2 0.1 -190. O. 0.1 7. O. 0.102 -189.199 801.708 0.102 6.20184 -798.015 0.104 -186.788 1611. 39 0.104 3.8093 -1593.93 0.106 -182.744 2435.18 0.106 -0.169865 -2303.71 0.108 -177.036 3275.42 0.108 -5.71633 -3159.52 0.11 -169.634 4128.3 0.11 -12.79 -3908.15 0:112 -160.522 4981.47 0.112 -21.3178 -4'609.96 0.114 -149.723 5811.86 0.114 -31.1013 -5238.95
PAGE 98

88 For Pm12 -0.18 p.u For Pm12 -0.2 p.u TIME Xl X2 TIME Xl X2 0.1 50. O. 0.1 70. O. 0.102 49.3421 -656.912 0.103 66.9611 -682.49 0.104 47.3564 -1329.6 0.106 65.8668 -1403.26 0.106 44.0075 -2023.77 0.109 60.4855 -2199.91 0.108 39.2397 -2749.96 0.1l2 52.5569 -3106.68 0.11 32.9822 -3514.09 0.115 41.7098 -4148.32 0.112 25.1576 -4316.38 0.116 27.5286 -5326.92 0.114 15.6963 -5148.96 0.121 9.6511 -6600.34 0.116 4.55403 -5993. 0.124 -12.0591 -7056.06 0.118 -8.26209 -6816.24 0.127 -37.2844 -8906.23 0.12 -22.6665 -7572.59 0.13 -65.0603 -9521.16 0.122 -38.4701 -8205.36 0.133 -93.8087 -9539.25 0.121 -55.366 -8655.61 0.136 -121.694 -8961. 39 0.126 -72.9377 -8874.57 0.139 -147.144 -7952.33 0.128 -90.6923 -8636.49 0.142 -169.209 -6742.78 0.13 -108.115 -6546.2 0.145 -187.597 -5527.22 0.132 -124.731 -8037.67 0.148 -202.485 -4422.08 0.134 -140.156 -7364.29 0.151 -214.289 -3474.67 0.136 -154.118 -6585.62 0.154 -223.496 -2688.91 0.138 -166.464 -5755.89 0.157 -230.567 -20n.43 0.14 -177 .136 -4917.26 0.16 -235.899 -1525.0S 0.142 -186.147 -4097.86 0.163 -239.81 -1095.7 0.144 -193.551 -3313.03 0.166 -242.542 -735.061 .0.146 -199.425 -2566.02 0.169 -244.266 -421.266 0.148 -203.848 -1660.82 0.172 -245.097 -134.561 0 .15 -206.869 -1181.61 0.175 -245.004 143.161 0 .152 -208.601 -529.656 0.178 -214.227 430.743 0.154 -209.015 115.431 0.181 -242.4'12 745.715 0.156 -208.138 762.687 0.184 -239.705 1108.19 0.158 -205.955 1424.01 0.187 -235.753 1540.12 0.16 -202.425 2110.33 0.19 -230.371 2065.88 0.162 -197.491 2830.51 0.193 -223.238 2711.59 0.164 -191.077 3589.88 0.196 -213.956 3502.29 0 .166 -183.105 4388.16 0.199 -202.061 4454.95 0.168 -173.504 5216.84 0.202 -167.066 5564.64 0.17 -162.23 6056.24 0.205 -168.5. 63 6782.23 0.172 -149.293 6673.23 0.206 -146.362 7988.77 0.174 -134.783 7620.87 0.211 -120.837 6967.62 0.1.76 -116.894 8241. 96 0.214 -92.097 9540.36 0.178 -101.939 8677.66 0.217 -64.1512 9510.28 0.1.8 -84.3395 8880.08 0.22 -36.4349 6878.35 0.182 -66.5911 8624.91 0.223 -11.3105 7818.86 0.184 -49.2079 8518.64 0.226 10.2799 6559.25 0.186 -32.6608 7996.53 0.229 26.0358 5287.55 0.108 -17.3297 7312.56 0.232 42.1047 4112.96 0.19 -3.47866 6526.14 0.235 52.8524 3075.81 0.192 8.742.62 5690.83 0.238 60.6946 2173.01 0.194 19 28 4847.93 0.241 66. 1379.35 0.196 28.1477 4024.67 0.244 69.0455 660.371 0.198 35.4017 3235.64 0.247 69.9993 -21.5226 0.2 41.1164 2485.47 0.25 60.9152 -704.638 0.202 45.368 1771.7 0.253 65.732 -1427.27 0.204 48.2227 1087.10 0.256 60.2745 -2226.09 0.206 49.7296 421.927 0.259 52.2592 -3137.64 0.208 49.9156 -235.59 0.262 41.3126 -4183.77 0.21 48.7843 -897.445 0.265 27.0100 -5366.34 0.212 46.3151 -1575.41 0.268 9.01984 -6641.36 0.214 42.4649 -2260.02 0.271 -12.8103 -7094.6 0.216 37.1717 -3019.4 0.274 -38.1348 -6933.62 0.218 30.3613 -3797.54 0.277 -65.9684 -9531. 41 0.22 21.957 -4612.11 0.28 -94.7175 -9529.55 0.222 11.6959 -5451.63 0.283 -122.547 -8934. 8 0.224 0.149153 -6292.66 0.286 -147.9 -7915.0 0.226 -13.2512 -7097.89 0.289 -169.85 -6703.45 0.228 -28.1851 -7816.8 0.292 -100.122 -5490.04 0 .23 -44.4218 -8390.68 0.295 -202.905 -4389.48 0.232 -61.613 -8762.66 0.298 -214.619 -3447.31 0.234 -79.3094 -8890.91 0.301 -223.751 -2666.47 0.236 -97.0033 -8760.2 0.304 -230.761 -2029.17 0.238 -114.166 -6386.47 0.307 -236.044 -1510.15 0.24 -130.415 -7812.09 0.31 -239.914 -1083.35 0.242 -145.34 -7094.19 0.313 -242.612 -724.523

PAGE 99

APPENDIX C THE ENERGY METRIC ALGORITHM This topic is from the article "The generation of Liapunov function in control theory by an energy metric algorithm" Edward T. Wall, Professor. Department of Electrical Engineering, University of Colorado at Denver. The basic procedure involved consists of the following six steps: Step 1: Describe the system as a set of first-order nonlinear differential equations. X.= F.(X), 1 < i:5 n I I --(1) Where X is the vector (Xl' X2 ... X n). Step 2: Create a set of differential equations of integral curves by forming quotients of first-order nonlinear differential equations and eliminating dt in each case. dX. F. (X) __ 1 I -,(j>i), and (j;ei) (2) dX. F.(X) J J -There wIII be n(n-1) h t sue equa Ions. 2

PAGE 100

90 Step 3: The differential equations of integral curves are converted to the same number of differential one-forms by multiplying and clearing the denominator terms F.(X) dX. F.(X) dX. = 0, J-1 1-1 (3) Step 4: By addition and substitution these differential one-forms are reduced to a single one-form. m = m (X)dX + m (X)dX + 1-1 2-2 + m (X)dX n-n (4) where the m's are determined by the additions and substitutions. Step 5: Perform a line integration of the one-forms The path of integration is chosen along the state coordinates to form and "energy-like" scalar function directly from the system equation. x O ... OldT,+J 'm,(x,. T,. O ... OldT,+ ... o x ... + J nmn(X'.x' .... xn_'.TnldTn o The result of the line integration is the V function. Step 6: The total derivative of V with respect to t is taken to obtain V. The proper theorems are then considered to determine stability or a region of stability. (5)

PAGE 101

APPENDIX D NOTATION Symbol Quantity A B state vectors C constant of integration C constant of Liapunov function D damping constant Efd generator field voltage E ex exci ter vo ltage -Eq voltage proportional to field flux linkage -Eq flux decay Ei generator ith, terminal voltage Ej generator jth, terminal voltage Fi set of nonlinear functions Gii the sum of the admittance connected to node i H H constant of the machine Hessian matrix l< maximum value of V(X) on the stability boundary M inertia constant of the machine MW megawatts MVA megavoltampere

PAGE 102

MVAR max. min. MAG. n NO. P Pa Pc Pe Pm Pm Pmax. p.u rad r :R Q,S sec. t t cr tic V megavars maximum minimum magintude state vector number real power accelerating power dissipating power electrical power governor action mechanical power maximum power transfer per unit radian state vector region liner operator second time critical clearing time fault clearing time Liapunov function 92

PAGE 103

v Vq Vq w lU X XT Xl X2 X-d XL Yij 0:, f3 time derivative of V quadratic term of Liapunov function time derivative of Vq machine speed in rad./sec. scalar function "energy-like" state vector transpose of vector X state variable in deg. state variable in deg./sec. direct-axis transient reactance transformer reactance transmission line reactance admittance between node i and node j load impedance completing the square constant 93 cP the sum of the steady-state torque angle and the real component of transfer admittance S system torque angle S angular speed deviation pi=3.141592654 rad. Oij angle between node i and node j A engenvalue

PAGE 104

94 Subscripts 1 denotes to generator one 2 denotes to generator two 12 denotes the system equivalent power value o denotes pre-fault system value (degree or initial) ss denotes the steady-state post-fault system value