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.