7

SYSTEM STABILITY

7.1 Introduction
Traditionally, the concept of stability has been used to indicate that the synchronous
rotating components of a power system remain in synchronism with each other during
load changing conditions and following system disturbances. This interpretation is thus
related to relative rotor angle variations.
With the increasing pressure to use the transmission system further, and thus post-
pone transmission reinforcements, voltage instability is becoming a greater challenge.
Voltage stability is often assessed with the help of load flow models, complemented
by the use of P - V and V-Q curves. In general, however, the stability, both of the
voltage and of the rotor angle types, following some predetermined operating condition
is a dynamic problem and requires more elaborate plant component models than those
discussed in previous chapters. It is normally assumed that prior to the dynamic analysis,
the system is operating in the steady state and that a load flow solution is available.
In this respect, the mechanisms of voltage and rotor angle stability are difficult to
separate in other than extreme situations; conceptually, the former relates mostly to
the stability of loads, whereas the latter is a problem of the synchronous generators. In
a large interconnected system, voltage collapse of a load area is possible without loss
of synchronism of the generators.
The time frame of the phenomena is the determining factor for the analytical assess-
ment. The region between zero and, say, 10 s subsequent to a disturbance is referred
to as transient stability and the solution is obtained in the time domain. The solution
then provides important information of first and multiple generator swings as well as
transient voltage stability.
Longer-term stability can be assessed in either the frequency or the time domain,
although the latter is more often used in practice. In this chapter, longer-term stability
is treated as an extension of transient stability and is thus solved in the time domain.
Such extension normally requires modification of some plant component models and
often the introduction of new models, but because of the smaller perturbations and
longer study duration, the small time constant effects can be ignored.
The transient stability program described in sections 7.1-7.8 is sufficient for many
basic stability studies. It is more than adequate when first swing stability is being
evaluated and the machine detail and controllers will allow second and subsequent
swing stability to be examined also.
However, if synchronous machine saturation or compound thermal turbines have
to be modelled, it will be necessary to incorporate section 7.9.1 and 7.9.2 into the
230 7 SYSTEM STABILITY

program. The structure of this program should allow changes of this sort to be made
quite easily. Similarly, if other system components are to be included, this can be done
without difficulty.
A transient stability program need not necessarily contain all the models described
in order to describe completely a power system. Conversely, a program containing
all these refinements is not necessarily adequate for a particular system. It must be
anticipated that transient stability programs will be continuously refined as tighter
operating constraints coupled with new control strategies are introduced.

7.1.1 The form of the equations

To a greater or lesser extent, all system variables require time to respond to any
change in operating conditions and a large set of differential equations can be written
to determine this response. This is impractical, however, and many assumptions must
be made to simplify the system model. The assumptions made depend on the problem
being investigated and no clear definitive model exists.
A major problem with a time domain solution is the stiffness of the system
(Appendix I), That is, the time constants associated with the system variables vary enor-
mously. When only synchronous machines are being considered, rotor swing stability
is the principal concern. The main time constants associated with the rotor are of the
order of 1-10 s. The form of the solution is dominated by time constants of this order
and smaller or greater time constants have less significance.
The whole of the a.c. transmission network responds rapidly to configurational
changes as well as loading changes. The time constants associated with the network
variables are extremely small and can be considered to be zero without significant loss
of accuracy. Similarly, the synchronous machine stator time constants may be taken
as zero. The relevant differential equations for these rapidly changing variables are
transformed into algebraic equations.
When the time constant is large, or the disturbance is such that the variable will
not change greatly, the time constant may be regarded as infinite, i.e. the variable
becomes a constant. Excitation voltage, or mechanical power to the synchronous
machine, may often be treated as constant in short duration studies without appre-
ciable loss of accuracy. Depending on how the computer program is written, variables
that become constant may be treated by either:
(1) retaining the differential equation but assigning a very large value to the relevant
time constant;
(2) removing the differential equation.
For flexibility, both methods are usually incorporated in a program.
A system which after these initial simplifications contains (Y differential variables,
contained in the vector Y,, and 3, algebraic variables, contained in the vector Xg,may
be described by the matrix equations:

where p denotes the differential operator d/dt

 7.2 SYNCHRONOUS MACHINES-BASIC MODELS 231

7.1.2 Frames of reference

The choice of axes, or frame of reference, in which the system equations are formulated,
is of great importance as it influences the analysis.
For synchronous machines, the most appropriate frame of reference is one that is
attached to the rotor, i.e. it rotates at the same speed as the rotor. The main advantage
of this choice is that the coefficients of the equations developed for the synchronous
machine are not time dependent. The major axis of this frame of reference is taken as
the rotor pole or direct axis. The second axis lies 90 (electrical) from each pole and
is referred to as the quadrature axis.
In the dynamic state, each synchronous machine is rotating independently from the
others and transforming between synchronous machine frames through the network
is difficult. This is overcome by choosing an independent frame of reference for the
network and transforming between this frame and the synchronous machine frames at
the machine terminals. The most obvious choice for the network is a frame of reference
that rotates at synchronousspeed. The two axes are obtained from the initial steady
state load flow slack busbar. Although the network frame is rotating synchronously,
this does not stop each nodal voltage or branch current from having an independent
frequency during the dynamic analysis.

7.2 Synchronous Machines-Basic Models

7.2.1 Mechanical equations
The mechanical equations of a synchronous machine are very well established [ 1,2]
and need be only briefly outlined. Three basic assumptions are made in deriving the
equations:

(1) Machine rotor speed does not vary greatly from synchronous speed (1.O PA).
(2) Machine rotational power losses due to windage and friction are ignored.
(3) Mechanical shaft power is smooth, i.e, the shaft power is constant except for the
results of speed governor action.

Assumption (1) allows per unit power to be equated with per unit torque. From
Assumption (2), the accelerating power of the machine (Pa) is the difference between
the shaft power (Pm)as supplied by the prime mover or absorbed by the load and the
electrical power (Pe).The acceleration ( a ) is thus:

where Mg is the angular momentum.

The acceleration is independent of any constant speed frame of reference and it
is convenient to choose a synchronously rotating frame to define the rotor angle (6).
Thus:
(7.4)
232 7 SYSTEM STABILITY

The angular momentum may be further defined by the inertia constant Hg (measured
in MWSMVA), which is relatively constant regardless of the size of the machine, i.e.:

where f o is the system base frequency.

Eddy currents induced in the rotor iron or in the damping windings produce torques
which oppose the motion of the rotor relative to the synchronous speed. A deceleration
power can be introduced into the mechanical equations to account for this damping,
giving
d2S 1
- = - ( P m - Pe - Do"). dr
dt2 Mg
The damping coefficient (Da), measured in W.rad-' has been largely superseded
by a synchronous machine model that includes the subtransient effect of the damper
windings in the electrical equations, but it is still used in some programs.
Two single-order ordinary differential equations may now be written to describe the
mechanical motion of the synchronous machine, i.e.
1
po = -(Pm - Pe - Da(w - 2nfo)), (7.7)
Mg

7.2.2 Electrical equations

The derivation of equations to account for flux changes in a synchronous machine has
been given by Concordia [3] and Kimbark [4]. A brief outline only will be given in
this section, so that various electrical quantities may be defined and phasor diagrams
constructed. The approximations made in the derivation are as follows:

(1) The rotor speed is always sufficiently near 1.0 p.u. that it may be considered a
constant.
(2) All inductances defined in this section are independent of current. The effects due
to saturation of iron are considered in section 7.9.
(3) Machine winding inductances can be represented as constants plus sinusoidal
harmonics of rotor angle.
(4) Distributed windings may be represented as concentrated windings.
( 5 ) The machine may be represented by a voltage behind an impedance.
( 6 ) There are no hysteresis losses in the iron, and eddy currents are only accounted
for by equivalent windings on the rotor.
(7) Leakage reactance only exists in the stator.
Using these assumptions, classical theory permits the construction of a model for
the synchronous machine in the steady state, transient and subtransient states.
 7.2 SYNCHRONOUS MACHINES-BASIC MODELS 233

The per unit system adopted is normalized to eliminate factors of a,a, ?rand
turns ratio, although the term pruportiunal should be used instead of equal when
comparing quantities. Note that 1 p.u. field voltage produces 1.0 p.u. field current and
1.0 p.u. open-circuit terminal voltage at rated speed.

Steady state equations Figure 7.1 shows the flux and voltage phasor diagram for
a cylindrical rotor-synchronous machine in which all saturation effects are ignored.
The flux Ff is proportional to the field current If and the applied field voltage and it
acts in the direct axis of the machine. The stator open-circuit terminal voltage Ei is
proportional to Ff but lies on the quadrature axis. The voltage Ei is also proportional
to the applied field voltage and may be referred to as Ef.
When the synchronous machine is loaded, a flux F proportional to and in phase
with the stator current I is produced which, when added vectorially to the field flux
Ff,gives an effective flux Fe. The effective internal stator voltage El is due to Fe and
lags it by 90. The terminal voltage V is found from this voltage El by considering
the volt drops due to the leakage reactance X1 and armature resistance Ra. By similar
triangles, the difference between Ef and El is in phase with the IX1 volt drop and is
proportional to I. Therefore, the voltage difference may be treated as a volt drop across
an armature reactance Xu. The sum of Xa and X1 is termed the synchronous reactance.
For the salient pole synchronous machine, the phasor diagram is more complex.
Because the rotor is symmetrical about both the d- and q-axis, it is convenient to resolve
many phasor quantities into components in these axes. The stator current may be treated
in this manner. Although Fd will be proportional to l d and F, will be proportional to
I,, because the iron paths in the two axes are different, the total armature reaction
flux F will not be proportional to I nor necessarily be in phase with it. Retaining our
earlier normalizing assumptions, it may be assumed that the proportionality between
Id and Fd is unity but the proportionality between I, and F, is less than unity and is
a function of the saliency.

Direct axis
A
I

DFe
Figure 7.1 Phasor diagram of a cylindrical rotor synchronous machine in the steady state
234 7 SYSTEM STABILITY

Direct axis

Qa;sxi Quadiature

Figure 7.2 Phasor diagram of a salient pole synchronous machine in the steady state

In Figure 7.2, the phasor diagram of the salient pole synchronous machine is shown.
Note that the d and q axes armature reactances have been developed as in the cylindrical
rotor case. From these, direct and quadrature synchronous reactances ( X , and X , ) can
be established, i.e.
(7.9)
(7.10)
Ei - V , = Ral, -XdId, (7.11)
-vd = Raid + X,I,, (7.12)
where v d and V , are the axial components of the terminal voltage V .
In steady-state conditions it is quite acceptable to use, as the machine model, the
field voltage Efor the voltage equivalent to field current Ei behind the synchronous
reactances. In these circumstances, the rotor position (quadrature axis) with respect to
the synchronously rotating frame of reference is given by the angular position of Ef.
Only the salient pole machine will now be considered, as the cylindrical rotor model
may be regarded as a special case of a salient machine ( x d = X , ) .

Transient equations For faster changes in the conditions external to the synchronous
machine, the above model is no longer suitable. Due to the inertia of the flux linkages,
these changes cannot be reflected throughout the whole of the model immediately. It
is, therefore, necessary to create new fictitious voltages EL and Eb which represent
the flux linkages of the rotor windings. These transient voltages can be shown to exist
behind the transient reactances X& and X i .
Eb - v, = Ral, - x & I d , (7.13)
Eb - V d = R d d +XiI,. (7.14)
 7.2 SYNCHRONOUS MACHINES-BASIC MODELS 235

Figure 7.3 Phasor diagram of a synchronous machine in the transient state

The voltage Ei should now be considered as the sum of two voltages, Ed and E,,
and is the voltage behind synchronous reactance. In the previous section, where steady
state was considered, current flowed only in the field winding and, hence, in that case,
Ed = 0 and E , = Ei.
Where it is necessary to allow the rotor flux linkages to change with time, the
following ordinary differential equations are used:
pEb = (Ef - E ,)/ T'do = (Ef + (Xd - X&)ld - E&)/T&o, (7.15)

pE& = - E d / T b = ( - ( X q - Xi)', - E&)/Tbo. (7.16)

The phasor diagram of the machine operating in the transient state is shown in
Figure 7.3.
Subtransient equations Either deliberately, as in the case of damper windings, or
unavoidably, other circuits exist in the rotor. These circuits are taken into account if a
more exact model is required. The reactances and time constants involved are small and
can often be justifiably ignored. When required, the development of these equations is
identical to that for transients, and yields:
(7.17)
(7.18)
p E i = (Eb + (X&- x i ) I d - E t ) / T i o , (7.19)
pE&'= (E&+ (Xi- X:)I, - EZ)/T:o. (7.20)
The equations are developed assuming that the transient time constants are large
compared with the subtransient time constants. A phasor diagram of the synchronous
machine operating in the subtransient state is shown in Figure 7.4. It should be noted
that Equations (7.13) and (7.14) are now true only in the steady state mode of operation,
although once subtransient effects have decayed, the error will be small.
236 7 SYSTEM STABILITY

Quadrature
axis

Figure 7.4 Phasor diagram of a synchronous machine in the subtransient state

Machine models It is possible to extend the model beyond subtransient level but
this is seldom done in multi-machine programs. Investigations [ 5 ] using a generator
model with up to seven rotor windings have shown that using the standard machine
data the more complex models do not necessarily give more accurate results. However,
improved results can be obtained if the data, especially the time constants, are suitably
modified.
The most convenient method of treating synchronous machines of differing
complexity is to allow each machine the maximum possible number of equations
and then let the actual model used be determined automatically according to the data
presented.
Five models are thus possible for a four-winding rotor.

Model 1 -constant voltage magnitude behind d-axis transient reactance (X&)

requiring no differential equations. Only the algebraic Equations (7.13) and (7.14)
are used.
Model 2 -D-axis transient effects requiring one differential equation ( p E b ) .
Equations (7.13), (7.14) and (7.15) are used.
Model 3 -D- and q-axis transient effects requiring two differential equations
(PEL and pE&).Equations (7.13), (7.14), (7.15) and (7.16) are used.
Model 4 -D- and q-axis subtransient effects requiring three differential equations
( p E b , p E t and pE&).Equations (7.15). (7.17), (7.18), (7.19) and:
-(Xq - X t ) I q - E&
pE: = (7.21)
T$
are used. This last equation is merely Equation (7.16) with modified primes.
Whether it is a subtransient or transient equation is open to argument.
 7.3 SYNCHRONOUS MACHINE AUTOMATIC CONTROLLERS 237

Model 5 -D-and q-axis subtransient effects requiring four differential equations

(PEL, PEL, pE: and p E z ) . Equations (7.15), (7.16), (7.17), (7.18), (7.19) and
(7.20)are used.
The mechanical Equations (7.7) and (7.8) are also required to be solved for all these
models.
Groups of synchronous machines or parts of the system may be represented by
a single synchronous machine model. An infinite busbar, representing a large stiff
system, may be similarly modelled as a single machine represented by Model 1, with
the simplification that the mechanical Equations (7.7) and (7.8) are not required. This
sixth model is thus defined as:
Model 0 -Infinite machine-constant voltage (phase and magnitude) behind d-axis
transient reactance (Xb).Only Equations (7.13) and (7.14) are used.

7.3 Synchronous Machine Automatic Controllers

For dynamic power system simulations of 1 second or longer duration, it is necessary
to include the effects of the machine controllers, at least for the machine most affected
by the disturbance. Moreover, controller representation is becoming necessary, even
for first swing stability, with systems being operated at their limits with near critical
fault clearing times.
The two principal controllers of a turbine-generator set are the automatic-voltage
regulator (AVR) and the speed governor. The AVR model consists of voltage sensing
equipment, comparators and amplifiers controlling a synchronous machine which can
be generating or motoring. The speed governor may be considered to have similar
equipment but, in addition, it is necessary to take the turbine into account.

7.3.1 Automatic voltage regulators

Many different AVR models have been developed to represent the various types used
in a power system. The application of such models is difficult and a better approach is
to develop a single general purpose AVR model, on a similar basis to the synchronous
machine model. The model can then revert to any desired type by using the correct
data. The IEEE defined several AVR types [6], the main two of which (Type 1 and
Type 2) are shown in Figure 7.5.
A composite of these two AVR types can be constructed. This model may also
include a secondary signal, which can be taken from any source, but is usually either
machine rotor speed deviation from synchronous speed or rate of change of machine
output power. This model is shown in Figure 7.6 and has been found to be satisfactory
for all the systems studied so far. It is acknowledged that other AVR models may be
necessary for specific studies.
In many systems studied, the amount of data available for an AVR model is quite
small. The composite model can degenerate into a very simple model easily by
defaulting time constants to zero and gains to either zero, unity or an extremely large
value depending on their position.
238 7 SYSTEM STABILITY

n
Se= f(f) 4
Specified 0 vs
voltage Saturating
Regulator amplifier ,, exciter
Ef
t
Field
busbar voltage voltage

VO
0-
I kf.p ~

(1+ 7-r.p)
Other ,

f-$
(a) signals Feedback loop

(Ke+Te.P) Field
voltage

I kf.p 1,
Other (1+ rfp)(i+ w.p)
(b) signals Feedback loop

Figure 7.5 Block diagrams for two commonly used AVR models. (a) IEEE Type 1 AVR model;
(b) IEEE Type 2 AVR model. (0 1982 IEEE)

Specified vs
voltage

voltage

0-L

signals

Figure 7.6 Block diagram of a composite automatic voltage regulator model

The equations for the AVR model, shown in Figure 7.6, are as follows:
pVfl = (Vt - V>/Tr, (7.22)
PVU = ( K U ( 1 + Tb * p)Vh - Va)/Ta, (7.23)
subject to
IP V 4 i ~ m a x
 7.3 SYNCHRONOUS MACHINE AUTOMATIC CONTROLLERS 239

and

(7.24)
(7.25)

PVX = (KXVaux - V x ) / T x , (7.27)

~ V =O ((1 +Ty * p)Vx - VO)/TZ, (7.28)
Vh = V S -' V -'V + VO, (7.29)
Ve = SeEf, (7.30)
where Se = f (EB.
Vg = Ef [unless IEEE Type 2, when Vg = V a l , (7.31)
VaUx= a predefined signal
The IEEE [6] recommends that Se be specified at maximum field voltage (Semax)
and at 0.75 of maximum field voltage (Seo.75max). From this, Se may be determined
for any value of field voltage by either linear interpolation or by fitting a quadratic.
Where linear interpolation is used, Equation (7.30) may be transformed to:
V e = ( k ~. Ef - k2)Ef, (7.32)
where

or

A means of modelling lead-lag circuits, such as those in the regulator amplifier, the
stabilizing loop and the auxiliary signal circuits, is given at the end of this section.
Despite the advantages of one composite AVR model, if there are a great many
AVRs to be modelled, most of which have simple characteristics, then it is better to
make two models. One model, which contains only the commonly used parts of the
composite model, can then be dimensioned for all AVRs. The other model, which
contains only the less commonly used parts of the composite model, can be quite
small dimensionally. A connection vector is all that is necessary to interconnect the
two models whenever necessary.

7.3.2 Speed governors

For speed governors, as with AVRs, a composite model which can be reduced to any
desired level is the most satisfactory. The speed governor models recommended by the
IEEE [7] are shown in Figure. 7.7. It will be noticed that, if limits are not exceeded,
the two models are identical. The difference is due to the assumption that, in a hydro
240 7 SYSTEM STABILITY

Specified
power setting y Ps

(a) 2 n fo 4 Ref. speed I

Specified

Governor
speed power setting

Figure 7.7 Typical models of speed governors and valves. (a) Thermal governor and valve;
(b) Hydro governor and valve. 0 1982 IEEE

Boiler or
water power 7 Pb

2 n f, b Ref. speed R=
1
2 n foReg

Figure 7.8 Generalized model of a speed governor and valve

governor, gate servo and gate positions are the same. One model can be used for the
governors of both turbines provided that the limits are either internal or external to the
second transfer function block of Figure 7.8. Also, very little extra effort is required to
divorce the governor from the actual turbine power and keep it instead as a function
of valve position.
The equations of the speed governor shown in Figure 7.8 are:

pG1 = [R(1 + T2.p)(2nfo - - G l ] / T l ,

W) (7.33)
pG2 = (G1 - G2)/T3, (7.34)
GW = G2 + Gs. (7.35)

The valve/gate position setting (Gw) is subject to opening and closing rate limits
(Omaxand cmax,respectively) and to physical travel limits so that:

(7.36)
 7.3 SYNCHRONOUS MACHINE AUTOMATIC CONTROLLERS 241

The valve equation is:

Pgv = GV* Pb. (7.37)

For thermal turbines where a boiler is modelled in the steady state, Pb will be the
actual power delivered and Gs will be unity, i.e. the valve will be fully open. If a
boiler is not modelled, or a hydro turbine is being controlled, then in the steady state,
Pb will be the maximum output from the boiler or water system (i.e. maximum turbine
mechanical power output) and Gv, and hence Gs, will be such that Pgv is the actual
mechanical power output of the turbine.
This method of modelling a valve has the advantage that non-linearities between
valve position and power can be easily included and also the operation of the governor
and valve can be readily interpreted.
For a hydro governor where the limits are external, the model is as given in Equations
(7.33)-(7.37) but for a thermal governor, G2 is reset after the valve limits are applied
to be:
G2[1im)= Gv - Gs (thermal governor only). (7.38)

7.3.3 Hydro and thermal turbines

This section is restricted to the modelling of simple turbines only. Compound thermal
turbines may require a detailed model, as shown in section 7.9.2, but, for stability
studies of only 1 or 2 s duration, the effect of all but the high pressure (HP) turbine
can usually be ignored. The time constant associated with the steam entrained between
the HP turbine outlet and the IP or LP turbine inlet is usually very large (greater than
5 s) and the output from all turbines other than the HP turbine may be treated as
constant.
Simple linear models of hydro and thermal turbines are shown in Figure 7.9. The
hydro-turbine model includes the penstock which gives the characteristic lead-lag
response of this type of turbine. The model is generally sufficient for all hydro turbines
and, from Figure 7.9, the differential equation for the mechanical power output ( P m )
of the turbine is:
pPm = ((1 - T w . p)Pgv - P m ) / T 4 , (7.39)
with T 4 = 0.5Tw as a further close approximation.

Power delivered Mechanical power

(a) through gate to generator

Mechanical power
through valve to generator

PI IP and LP
turbine power

Figure 7.9 Simple linear models of turbines. (a) Hydro turbine; (b) Thermal turbine
242 I SYSTEM STABILITY

For the thermal turbine of Figure 7.9(b), this equation is:

pPm = ( K 1 . Pgv - P m ) / T 4 , (7.40)

with K 1 representing the fraction of power delivered by the HP turbine. For simple
turbines, K1 is thus unity. For compound turbines, the power ( P I ) from the IP and LP
turbines is obtained from:
PI = ( 1 - KI)Pmo, (7.41)

here, Pmo is the initial steady state mechanical power. Note that for this simple model,
the initial value of Pgv is PmIK1, even though all the steam passes through the valve.
Provided that the HP valve does not close fully then, rather than inject the power
from the IP and LP turbines as shown, it is easier to treat it as a simple turbine (PI = 0)
but with the speed regulation modified by:

(7.42)

7.3.4 Modelling lead-lag circuits

Lead-lag circuits may present a problem depending on the integration scheme adopted.
Where the differential equations are not used directly and the derivatives are not explic-
itly calculated, the following can be used to convert the model into a more acceptable
form.
For the circuit shown in block diagram form in Figure 7.10, the equation is:

(7.43)

This can be transformed to:

+
K * T 2 ( T l / T 2 ) T1 * p
vou, =-
T1 ( l+Tl.p
Vin 9

and then to:

(K;r2 K ( 1 - T2/T1)
vou, = -+ 1+Tl.p
Vin 3 (7.44)

which can be represented by the block diagram in Figure 7.11, and is a lag circuit in
parallel with a gain.
It is important to remember that the time constant T1 must be non-zero even if the
integration method can accommodate zero-time constants.

Figure 7.10 Typical lead-lag circuit block diagram

 1.4 LOADS 243

Figure 7.11 Modified block diagram of a lead-lag circuit

7.4 Loads
Early transient stability studies were concerned primarily with generator stability, and
little importance was attached to loads. In the two-machine problem, for example, the
remainder of the system, generators and loads were represented by an infinite busbar.
A great deal of attention has been given to load modelling since then.
Much of the domestic load and some industrial load consist of heating and lighting,
especially in the winter, and in early load models these were considered as constant
impedances. Rotating equipment was often modelled as a simple form of synchronous
machine and composite loads were simulated by a mixture of these two types of load.
A lot of work has gone into the development of more accurate load models. These
include some complex models of large induction motors, as described in section 7.9.3.
Most loads, however, consist of a large quantity of diverse equipment of varying levels
and composition and some equivalent model is necessary.
A general load characteristic [8] may be adopted, such that the MVA loading at a
particular busbar is a function of voltage ( V ) and frequency (f):
P = Kp(V)P" * ( f ) P f , (7.45)
Q = Kq(V)". ( f I q f , (7.46)
where K p and K q are constants which depend upon the nominal value of the variables
P and Q.
Static loads are relatively unaffected by frequency changes, i.e. pf = q f = 0, and
with constant impedance loads p v = qv = 2.
The importance of accurate load models has been demonstrated by Dandeno and
Kundur [8] when considering voltage sensitive loads. Figure 7.12 demonstrates the
power and current characteristics of constant power, constant current and constant
impedance loads. Berg [9] has identified the characteristic load parameters for various
homogeneous loads and these are given in Table 7.1. These characteristics may be
combined to give the overall load characteristic at a busbar. For example, a group of
n homogeneous loads, each with a characteristic of pvj and a nominal power of Pj
may be combined to give an overall characteristic of

(7.47)
j=l
The other three overall characteristics may be similarly determined.
244 7 SYSTEM STABILITY

Por Q

i'lw
Nominal IVI
voltage

pvor qv= 1

'KI /
pvor q v = 2 !

Nominal IVI
voltage

Figure 7.12 Characteristic of different load models. (a) Active and reactive power against
voltage; (b) Current against voltage

Table 7.1 Typical values of characteristic load parameters [9]

Load P" 4" Pf sf
Filament lamp 1.6 0 0 0
Fluorescent lamp 1.2 3.0 -1.0 -2.8
Heater 2.0 0 0 0
Induction motor half load 0.2 1.6 1.5 -0.3
Induction motor full load 0.1 0.6 2.8 1.8
Reduction furnace 1.9 2.1 -0.5 0
Aluminium plant 1.8 2.2 -0.3 0.6

7.4.1 Low-voltage problems

When the load parameters pv and qv are less than or equal to unity, a problem can
occur when the voltage drops to a low value. As the voltage magnitude decreases, the
current magnitude does not decrease. In the limiting case with zero-voltage magnitude,
a load current flows which is clearly irrational, given the non-dynamic nature of the
load model. From a purely practical point of view, then the load characteristics are only
 7.6 OVERALL SYSTEM REPRESENTATION 245

valid for a small voltage deviation from nominal. Further, if the voltage is small, small
errors in magnitude and phase produce large errors in current magnitude and phase.
This results in loss of accuracy and, with iterative solution methods, poor convergence
or divergence.
These effects can be overcome by using a constant impedance characteristic to
represent loads where the voltage is below some predefined value, for example 0.8 p.u.

7.5 The Transmission Network

It is usual to represent the static equipment which constitutes the transmission system
by lumped equivalent pi parameters independent of the changes occurring in the
generating and load equipment. This representation is used for multi-machine stability
programs because the inclusion of time varying parameters would cause enormous
computational problems. Moreover, frequency, which is the most obvious variable in
the network, usually varies by only a small amount and, thus, the errors involved are
small. Also, the rates of change of network variables are assumed to be infinite which
avoids the introduction of differential equations into the network solution.
The transmission network can thus be represented in the same manner as in the load
flow or short circuit programs, i.e. by a square complex admittance matrix.
The behaviour of the network is described by the matrix equation:
(7.48)
where []in,] is the vector of injected currents into the network due to generators and
loads and [V] is the vector of nodal voltages.
Any loads represented by constant impedances may be directly included in the
network admittance matrix with the injected currents due to these loads set to zero.
Their effect is thus accounted for directly by the network solution.

7.6 Overall System Representation

Two alternative solution methods are possible. The preferred method uses the nodal
matrix approach, while the alternative is the mesh matrix method.
Matrix reduction techniques can be used with both methods if specific network
information is not required, but this gives little advantage as the sparsity of the reduced
matrix is usually very much less.

7.6.1 Mesh matrix method

In this method, the system-loading components are treated as Thevenin equivalents
of voltages behind impedances. The network is increased in size to include these
impedances and the mesh impedance matrix of the increased network is created. This
is then inverted or the factorized form of the inverse determined.
The solution process is as follows:

(1) Calculate the Thevenin voltages of the system loading components by solving the
relevant differential and algebraic equations.
246 7 SYSTEM STABILITY

(2) Determine the network currents using the Y matrix or factors. As the network
current around a mesh containing the Thevenin voltage is the loading current,
this may affect the Thevenin voltage, in which case an iterative process will be
required.

7.6.2 Nodal matrix method

In this method, all network loading components are converted into Norton equivalents
of injected currents in parallel with admittance. The admittances can be included in
the network admittance matrix to form a modified admittance matrix which is then
inverted, or preferably factorized by some technique so that solution at each stage is
straightforward.
The following solution process applies:

(1) For each network loading component, determine the injected currents into the
modified admittance matrix by solving the relevant differential and algebraic
equations.
(2) Determine network voltages from the injected currents using the Z-matrix or
factors.
As the network voltages affect the loading components, an iterative process is often
required, although good approximations [8] can be used to avoid this.
With the nodal matrix method, the busbar voltages are available directly while the
branch currents can be calculated if necessary. With the mesh matrix method, the mesh
currents are available directly while the busbar voltages and branch currents can be
calculated, if necessary.
Although much work has been spent on the systematic construction of the mesh
impedance matrix, the nodal admittance matrix is easier to construct and has gained
wide acceptance in load flow and fault analysis. For this reason, the remainder of this
section will consider the nodal matrix method.

7.6.3 Synchronous machine representation in the network

The equations representing a synchronous machine, as defined in section 7.2, are given
in the form of Thevenin voltages behind impedances. This must be modified to a
current in parallel with an admittance by use of Nortons theorem. The admittance of
the machine thus formed may be added to the shunt admittance of the machine busbar
and treated as a network parameter. The vector [Ii,,] in Equation (7.48) thus contains
the Norton equivalent currents of the synchronous machines.
The synchronous machine equations are written in a frame of reference rotating with
its rotor. The real and imaginary components of the network equations, as illustrated
in Figure 7.13, are obtained from the following transformation:

(7.49)
 7.6 OVERALL SYSTEM REPRESENTATION 247

Network
Machine 1, imaginary axis
direct axis Machine
quodrature axis

> Network real axis

Figure 7.13 Synchronous machine and network frames of reference

This transformation is equally valid for currents as is the reverse transformation

(7.50)

When saliency exists, the values of Xz and Xi used in Equations (7.17) and (7.18)
and/or Xi and Xi used in Equations (7.13) and (7.14)are different. Therefore, the
Norton shunt admittance will have a different value on each axis, and when transformed
into the network frame of reference, will have time-varying components. However,
a constant admittance can be used, if the injected current is modified to retain the
accuracy of the Norton equivalent [lo]. This approach can be justified by comparing
the two circuits of Figure 7.14 in which y, is a time-varying admittance, whereas 70
is fixed.

Figure 7.14 Method of representing synchronous machines in the network. (a) Norton equiv-
alent of the synchronous machine. (b) Modified equivalent circuit
248 7 SYSTEM STABILITY

At any time t , the Norton equivalent of the machine is illustrated in Figure 7.14(a),
but the use of a fixed admittance results in the modified circuit of Figure 7.14(b).
The machine current is:
-
I = P,(E- v)= YO(E- v)+
(7.51)
where Tad, accounts for the fact that the apparent current source is not accurate in
this case.
The injected current into the network which includes 70is given by:
- -
linj = Iunadj + fadj, (7.52)
where - -
Inadj = YOE.
A suitable value for 70is found by using the mean of direct and quadrature admit-
tances, i.e.
- Ra - j X d q
Yo = (7.53)
Ra2 X j X : +
where
xdq = ?(xi
I + xy).
The unadjusted value of current injected into the busbar is:

1
(7.54)
Ra2 +X2X: -xdq EL
The adjusting current is not affected by rotor position in the machine frame of
reference but it is when considered in the network frame. From Equation (7.51) and
also Equations (7.17) and (7.18):

and transforming

(7.56)
Ra2 +X j X ; EL -V,
The total nodal injected current is, therefore:

(7.57)
 7.6 OVERALL SYSTEM REPRESENTATION 249

7.6.4 Load representation in the network

In the overall solution method, loads must be represented as injected currents into the
transmission network from which the terminal voltages can be calculated. A Norton
equivalent model of each load must, therefore, be created. In a similar way to that
adopted for synchronous machines, the Norton admittance may be included directly in
the network admittance matrix.
A constant impedance load is, therefore, totally included in the network admittance
matrix and its injected current is zero. This representation is extremely simple to imple-
ment, causes no computational problems and improves the accuracy of the network
solution by strengthening the diagonal elements in the admittance matrix.
Non-impedance loads may be treated similarly. In this case, the steady state values of
voltage and complex power obtained from the load flow are used to obtain a steady state
equivalent admittance (70) which is included in the network admittance matrix [ Y ] .
During the stability run, each load is solved sequentially along with the generators,
etc, to obtain a new admittance (Y), i.e.:

p = - s* (7.58)
IVI2
The current injected into the network thus represents the deviation of the load char-
acteristic from an impedance characteristic.
(7.59)
By converting the load characteristic to that of a constant impedance, when the
voltage drops below some predetermined value (Vlow), as described in section 7.4, the
injected current is kept relatively small. An example of a load characteristic and its
corresponding injected current is shown in Fig. 7.15.
In an alternative model, the low voltage impedance is added to the network and the
injected current compensates for the deviation from the actual characteristic. In this
case, there is a non-zero injected current in the initial steady state operating condition.

7.6.5 System faults and switching

In general, most power system disturbances to be studied will be caused by changes in
the network. These changes will normally be caused by faults and subsequent switching
action, but occasionally the effect of branch or machine switching will be considered.

Faults Although faults can occur anywhere in the system, it is much easier compu-
tationally to' apply a fault to a busbar. In this case, only the shunt admittance at the
busbar needs to be changed, i.e. a modification to the relevant self-admittance of the
Y matrix. Faults on branches require the construction of a dummy busbar at the fault
location and modification of the branch data, unless the distance between the fault
position and the nearest busbar is small enough to be ignored.
The worst case is a three-phase zero-impedance fault and this involves placing an
infinite admittance in parallel with the existing shunt admittance. In practice, a non-zero
but sufficiently low fault impedance is used, so that the busbar voltage is effectively
250 I SYSTEM STABILITY

III 4

Constant power

Figure 7.15 Load and injected currents for a constant power type load with low voltage adjust-
ment. (a) Load current; (b) Injected current

brought to zero. This is necessary to meet the requirements of the numerical solution
method.
The application or removal of a fault at an existing busbar does not affect the
topology of the network and, where the solution method is based on sparsity exploiting
ordered elimination, the ordering remains unchanged and only the factors required for
the forward and backward substitution need to be modified. Alternatively, the factors
can remain constant and diakoptical techniques [ 111 can be used to account for the
network change.

Branch switching Branch switching can easily be carried out by either modifying
the relevant mutual- and self-admittances of the Y matrix or using diakoptical tech-
niques. In either case, the topology of the network can remain unchanged, as an open
branch is merely one with zero admittance. While this does not fully exploit sparsity,
the gain in computation time by not reordering exceeds the loss by retaining zero
elements, in almost all cases.
The only exception is the case of a branch switched into a network where no inter-
connection existed prior to that event. In this case, either diakoptical or reordering
techniques become necessary. To avoid this problem, a dummy branch may be included
with the steady state data; it should be of sufficiently high impedance so that the power
 7.6 OVERALL SYSTEM REPRESENTATION 251

flow is negligible under all conditions, or alternatively, the branch resistance may be
set negative to represent an initial open circuit. A negative branch reactance should
not be used as this is a valid parameter where a branch contains series capacitors.
If a fault occurs on a branch but very close to a busbar, the non-unit protection at
the near busbar will normally operate before that at the remote end. Therefore, there
will be a period when the fault is still being supplied from the remote end. There are
two methods of accounting for this type of fault.
The simplest method only requires data manipulation. The fault is initially assumed
to exist at the local busbar rather than on the branch. When the specified time for the
protection and local circuit breaker to operate has elapsed, the fault is removed and
the branch on which the fault is assumed to exist is opened. Simultaneously, the fault
is applied at the remote busbar but, in this case, with the fault impedance increased by
the faulted branch impedance; similarly, the fault is maintained until the time specified
for the protection and remote circuit breaker to operate has elapsed.
The second method is generally more involved but it is better when protection
schemes are modelled. In this case, a dummy busbar is located at the fault position,
even though it is close to the local busbar and a branch with a very small impedance
is inserted between the dummy busbar and the local busbar. The faulted branch then
connects the dummy busbar to the remote busbar and the branch shunt susceptance
originally associated with the local busbar is transferred to the dummy busbar. This
may all be done computationally at the time when the fault is being specified. The
two branches can now be controlled independently by suitable protection systems. An
advantage of this scheme is that the fault duration need not be specified as part of
the input data. Opening both branches effectively isolates the fault, which can remain
permanently attached to the dummy busbar or, if auto-reclosing is required, it can be
removed automatically after a suitable deionization period.
The second method will give problems if the network is not being solved by a
direct method. During the iterative solution of the network, slight voltage errors will
cause large currents to flow through a branch with a very small impedance. This will
slow convergence and, in extreme cases, will cause divergence. With a direct method,
based on ordered elimination, an exact solution of the busbar voltage is obtained for the
injected currents specified at that particular iteration. Thus, provided that the impedance
is not so small that numerical problems occur when calculating the admittance, and the
subsequent factors for the forward and backward substitution, then convergence of the
overall solution between machines and network will be unaffected. The value of the
low impedance branch between the dummy and local busbars may be set at a fraction
of the total branch impedance, subject to a minimum value. If this fraction is under
MOO, the change in branch impedance is very small compared with the accuracy of
the network data input and it is unnecessary to modify the impedance of the branch
from the remote to the dummy busbar.

Machine switching Machine switching may be considered either as a network or

as a machine operation. It is a network operation if a dummy busbar is created to
which the machine is connected. The dummy busbar is then connected to the original
machine busbar by a low impedance branch.
Alternatively, it may be treated as a machine operation by retaining the original
network topology. In this case, when the machine is switched out, it is necessary to
252 7 SYSTEM STABILITY

remove its injected current from the network solution. Also, any shunt admittance
included in the network Y matrix, which is due to the machine, must be removed.
Although a disconnected machine can play no direct part in system stability, its
response should still be calculated as before, with the machine stator current set to zero.
Thus machine speed, terminal voltage, etc., can be observed even when disconnected
from the system and in the event of reconnection, sensible results are obtained.
Where an industrial system is being studied, many machines may be disconnected
and reconnected at different times as the voltage level changes. This process will require
many recalculations of the factors involved in the forward and backward substitution
solution method of the network. However, these can be avoided by using the method
adopted earlier to account for synchronous machine saliency. That is, an appropriate
current is injected at the relevant busbar, which cancels out the effect of the shunt
admittance.

7.7 Integration
Many integration methods have been applied to the power system transient stability
problem and the principal methods are discussed in Appendix 11. Of these, only three
are considered in this section. They are simple and easily applied methods which have
gained wide acceptance. The purpose of the third method is not to provide another
alternative but to clarify the differences between the other two methods.
Explicit Runge-Kutta methods have been used extensively in transient stability stu-
dies. They have the advantage that a 'packaged' integration method is usually available or
quite readily constructed and the differential equations are incorporated with the method
explicitly. Stability problems, arising from the introduction of more detailed system
component models with very small time constants, have caused interest in other methods.
Fourth-order methods ( p = 4) have probably been the most popular and among these
the Runge-Kutta-Gill method has the advantage that round-off error is minimized.
With reference to Equations (11.21)-(11.23), for this method the number of function
substitutions is four (v = 4) and:
W I = 1/6,

~2 = (2 - f i ) / 6 , (7.60)
w3 = (2 + J2))/6,
~4 = 1/6,

(7.61)
k3 = h f ( t , + h / 2 , yn + (JZ- l)k1/2 + (2 - &)k*/2),
k4 = h f ( t , + h, y, - f i k 2 / 2 + (2 + JZ)k3/2).
The characteristic root of this fourth-order method, when applied to Equation
(Kll), is:
~1 = 1 + + +
hh {h2A2 $h3A3 &h4A4, + (7.62)
and to ensure stability, the step length h must be sufficiently small that ZI is less than
unity.
 7.7 INTEGRATION 253

The basic trapezoidal method is very well known, having been established as a useful
method of integration before digital computers made hand calculation redundant.
More recently, an implicit trapezoidal integration method has been developed for
solving the multi-machine transient stability problem [lo], and has gained recognition
as being very powerful, having great advantages over the more traditional methods.
The method is derived from the general multistep equation given by Equation (11.10)
+
with k equal to unity and is thus a single-step method. The solution at the end of n 1
steps is given by:
hfl+l
Yfl+I = Yfl + --(pyn+l
2 + pyfl). (7.63)

It has second-order accuracy with the major term in the truncation error being - A h 3 .
The characteristic root when applied to Equation (11.11) is:

ZI = 1 - 26,+1, (7.64)

where
(7.65)

If R e @ ) c 0, then 0 Ibfl+l 5 1.0 and lzl I I1.0. The trapezoidal method is, there-
fore, A-stable, a property which is shown in Appendix I1 to be more important in
the solution process than accuracy. The trapezoidal method is linear and thus in a
multivariable problem, like power system stability, the method is C-stable.
It can be shown that an A-stable linear multistep method cannot have an order of
accuracy greater than two, and that the smallest truncation error is achieved by the
trapezoidal method. The trapezoidal method is thus the most accurate C-stable finite
difference method possible.
The method, as expressed by Equation (7.63). is implicit and requires an iterative
solution. However, the solution can be made direct by incorporating the differential
equations into Equation (7.63). Rearranging forms algebraic equations, as described in
Appendix 11.
For example, consider the trivial transfer function shown in Figure 7.16. The differ-
ential equation for this system is given by:

P Y ( 0 = (G * t ( t >- y ( t ) ) / T , (7.66)
with the input variable being denoted by z to indicate that it may be either integrable
or non-integrable.
The algebraic form of Equation (7.66) has a solution at the end of the (n 1)th+
step of
Yn+l = cfl+1 %+I+ * Zn+l* (7.67)

Figure 7.16 Simple transfer function

254 7 SYSTEM STABILITY

where

(7.70)
Provided that the step length h remains constant, it is unnecessary to re-evaluate b
or m at each step. That is:

(7.71)

There is little to be gained by this, however, as it is a simple process and it is often

desirable to change h during a study.
A comparison between the Runge-Kutta-Gill and the trapezoidal methods when
used to solve two power system transient stability problems is given in Tables 7.2 and
7.3. The comparison is made in terms of maximum error (based on results using very
small step lengths) and central processor unit (CPU) execution time.
The advantages of the C-stable trapezoidal method are apparent from both tables
but the results are sufficiently different to show that an absolute comparison between
methods cannot be made. The non-linearity of the equations in any system also affects

Table 7.2
Step length Runge-Kutta-Gill Trapezoidal Backward Euler
(ms) Max. error CPU time Max. error CPU time Max. error CPU time
(deg) (S) (deg) 6) (deg) 6)
100.0 - - 2.2 0.26 - -
50.0 - - 0.7 0.27 - -
25.0 21.0 0.43 0.1 0.29 5.7 0.4 1
10.0 13.0 0.72 - 0.49 2.4 0.47
5.0 7.8 1.18 - 0.69 1.3 0.67
2.0 3.7 2.57 - 1.34 .5 1.31
1.o 1.9 4.88 - 2.42 .2 2.35
0.5 1.o 9.52 - 4.60 - 4.2
0.2 0.4 24.19 - - - 10.58
0.1 0.2 47.95 - - - -

Table 7.3
~ _ _ _

Step length Runge-Kutta-Gill Trapezoidal Backward Euler

(ms) Max. error CPU time Max. error CPU time Max. error CPU time
(deg) 6) (deg) (s) (deg) (s)
10.0 8.6 1.67 0.5 2.37 - -
5.0 4.4 3.06 0.1 2.31 8.5 2.76
2.0 1.7 7.24 - 3.74 3.8 364
1.o 1.2 14.19 - 7.12 1.8 6.80
0.5 0.9 28.00 - 13.88 0.6 13.24
 7.7 INTEGRATION 255

the errors obtained. CPU time using the Runge-Kutta-Gill method is a function of
the step length but this is not so with the trapezoidal method. For very small step
lengths, only one iteration per step is needed using the trapezoidal method, but as the
step length increases so does the number of iterations. The relationship between step
length and iterations is non-linear, with the result that there is an optimum step length
in which the iterations per step are small but greater than one.
For comparison, the backward Euler method is also included. This is a first-order
method with the solution given by:

Y ~ +=
I Yn +h PY~+I 9 (7.72)
and the characteristic root when applied to Equation (11.11) is:
Z] = 1/(1 - hA), (7.73)
Despite the three orders of accuracy difference between it and the Runge-Kutta-Gill
method, the backward Euler method compares well.
The results for the trapezoidal and backward Euler methods were obtained using
linear extrapolation of the non-integrable variables at the beginning of each step. This
required the storing of machine terminal voltages and currents together with other
non-integrable variables obtained at the end of the previous step.

7.7.1 Problems with the trapezoidal method

Although the trapezoidal method is %stable and the step length is not constrained
by the largest negative eigenvalue, the accuracy of the solution corresponding to the
largest negative eigenvalues will be poor if a reasonable step length is not chosen.
With the backward Euler method, the larger the step length the smaller the charac-
teristic root
ZI(BE) +. 0 as hA -+ -00, (7.74)
whereas, for the trapezoidal method

ZI(TRAP) + -1 as hA -+ -00, (7.75)

For small step lengths, the characteristic roots of both methods tend towards, but
never exceed, unity (positive).

ZI(BE) and ZI(TRAP) + +1 as hA + 0. (7.76)

The effect of too large a step length can be shown in a trivial but extreme example.
The system in Figure 7.16 and Equation (7.66) with a zero time constant T, and unity
gain G, is such an example.
If the input z ( t ) is a unit step function from an initial value of zero, then with a zero-
time constant, the output y ( t ) should follow the input exactly, i.e. a constant output of
unity. In fact, the output oscillates with yl = 2, y2 = 0, y3 - 2, etc.
Table 7.4 shows the effect of different step lengths on this simple system with a
non-zero time constant T. This table shows that oscillations occur when hX is smaller
than -2, i.e. when the characteristic root tl is negative. The oscillations decay with a
256 7 SYSTEM STABILITY

Table 7.4 The effect of different step lengths on the solution of a simple system (Figure 7.16)
by the trapezoidal method

Step hh = -0.5 hh = -2.0 hh =z -8.0 hh = -32.0

No' Trap Exact Trap Exact Trap Exact Trap Exact
method soh method soh. method soln. method soln.
0 0 0 0 0 0 0 0 0
1 0.4000 0.3935 1.0000 0.8664 1.6000 0.9997 1.8824 1.0000
2
3
4
5
0.6400
0.7840
0.8704
0.9222
0.6321
0.7769
0.8647
0.9179
I
1.0000
0.9817
0.9975
0.9997
1.0000
0.6400
1.2160
0.8704
1.0778
1.Oi)OO

1.0000
0.2215
1.6870
0.3938
1.5349
I
1.0000

rate dependent on hA, i.e. the rate is dependent on the magnitude of t i , It can also be
seen that accuracy is good provided that hA is greater than or equal to -0.5.
Oscillations are only initiated at a discontinuity. Provided that there is no step func-
tion input, the output of a transfer function with zero-time constant duplicates the
input.
The example given is an extreme case and, for the power system stability problem,
this usually only occurs in the input circuit of the AVR.
For the mechanical equation of the synchronous machine, the speed is given by:
1
pw = - ( P a ) , (7.77)
Mg
where Pa is the accelerating power given by Pa = P,,, - Pe, and the damping factor
Da is zero. Therefore, in this case

(7.79)

(7.80)

and oscillations do not occur. Da,when it does exist, is usually very small and any
oscillations will similarly be very small.
For the electrical equations of the synchronous machine, only the current can change
instantaneously, and the effect is not as pronounced as for a unit step function.
Techniques 1121 are available to remove the oscillations but they require a lot of
storage and it is simpler to reduce the step length.

7.7.2 Programming the trapezoidal method

There is no means of estimating the value of errors in the trapezoidal method but the
number of iterations required to converge at each step may be used as a very good
indication of the errors. As previously mentioned, the number of iterations increases
 7.7 INTEGRATION 257

more rapidly than the step length and thus the number of iterations is a good reference
for the control of the step length. It is suggested [13] to double the step length if the
number of iterations per step is less than 3 and to halve it if the number of iterations
per step exceeds 12. The resulting bandwidth (3-12) is necessary to stop constant
changes in the step length.
To avoid problems of step length chattering, a factor of about 1.5 (instead of 2)
may be used. Unfortunately, it is difficult to maintain a regular printout interval if a
non-integer factor is used.
Even using a step length changing factor of 2, it is difficult to maintain a regular
printout interval. Step halving can be carried out at any time but indiscriminate step
doubling may mean that there is no solution at the desired printout time. Doubling the
step length thus should only be done immediately after a printout and the step length
should not be allowed to exceed the printout interval.
On rare occasions, it is possible that the number of iterations at a particular step
greatly exceeds the upper desired limit. It can be shown that the convergence pattern is
geometric and usually oscillatory [13] after the first five or six iterations. Even when
diverging, the geometric and oscillatory pattern can be observed. Schemes can thus
be devised which estimate the correct solution. However, these schemes are relatively
costly to implement in terms of programming, storage and execution and a more
practical method is to stop iterating after a fixed number of iterations and start again
with a half-step. It is not necessary to store all the information obtained at the end of
the previous step, in anticipation of a restart, as this information is already available for
the non-integrable variables if an extrapolation method is being used at the beginning
of each step. Further, much information is available in the C and M constants of the
algebraic form of the integration method.
For example, with the two-variable problem given by Equation (7.66), if ~ ( t is) a
non-integrable variable, then its value at the end of the nth step Z, is stored. The value
of the integrable variable y ( t ) at the end of the nth step, y l r ,can be re-evaluated from
Equations (7.68) and (7.69) to be:

(7.81)

In only a few cases, where the differential equation is complex, need the value of
y, be stored at the beginning of each step. While the method requires programming
effort, it is very economical on storage and the few instances where it is used do not
affect the overall execution time appreciably.
Linear extrapolation of non-integrable variables at the beginning of each step is a
very worthwhile addition to the trapezoidal method. Although not essential, the number
of iterations per step is reduced and the storage is not prohibitive. Higher orders of
extrapolation give very little extra improvement and, as they are not effective until
some steps after a discontinuity, their value is further reduced.
It is only at the first step after a discontinuity that linear extrapolation cannot be
used. As this often coincides with a large rate of change of integrable variables, the
number of iterations to convergence can be excessive. This is overcome by automatic
step length reduction after a discontinuity. Two half-step lengths, before returning to
the normal step length, have been found to be satisfactory in almost all cases [13].
258 7 SYSTEM STABILITY

7.7.3 Application of the trapezoidal method

The differential equations developed in this chapter have all been associated with the
synchronous machine and its controllers. These equations can be transformed into alge-
braic form of the trapezoidal method given by Equation (7.67). While these algebraic
equations can be combined to make a matrix equation, this has little merit and makes
discontinuities such as regulator limits more difficult to apply.
In order to simplify the following equations, the subscripts on the variables have
been removed. It is rarely necessary to retain old values of variables and, where this is
necessary, it is noted. The variable values are thus overwritten by new information as
soon as they are available. The constants C and M associated with the algebraic form
are evaluated at the beginning of a new integration step and, hence, use the information
obtained at the end of the previous step.

Synchronous machine The two mechanical differential equations are (7.7) and (7.8)
and the algebraic form is:
w = C,,, + M,(Pm - Pe), (7.82)
where
C,,, = (1 - 2M,Du)w + M,,,(Pm - Pe + 4 nf o Da)
M , = h/(2Mg +~ D u ) ,
and also
8 = Cs + Ms(w), (7.83)
where
cs = s + Ms(w - 4n f o ) ,
Ma = 0.5h.
It would be possible to combine these equations to form a single simultaneous
solution of the form:
+
S = C;i M i ( P m - P e ) , (7.84)
where

but machine speed w is a useful piece of information and would still require evaluation
in most problems.
It is also more convenient to retain the electrical power ( P e ) as a variable rather
than attempt to reduce it to its constituent parts:

Thus, Pe is extrapolated after C, and M , have been evaluated.

 7.7 INTEGRATION 259

The mechanical power Pm is an integrable variable which, in the absence of a speed

governor model for the machine, is constant.
There are four electrical equations associated with the change in flux in the
synchronous machine, i.e. Equations (7.15), (7.16) (7.19) and (7.20). In algebraic form,
these equations are:
+
Eb = c9 M q ( E f + ( x d - & ) I d ) , (7.86)

where

(7.87)

where

(7.88)

where

(7.89)

Cdd = (1 - 2Mdd)Ei + Mdd(E& - (xb - xt)19),

Mdd = h / ( 2 T i 0 +h).
Synchronous machine controller limits There are usually limits associated with
AVRs and speed governors and these require special consideration when applying
the algebraic form of the trapezoidal rule. It is best to ignore the limits at first and
develop the whole set of 'limitless' equations. Rather than confuse this discussion, it
is easier to consider a simple AVR system, as shown in Figure 7.17, for which can be
written:
pvout = (Gl(Vin - Vfb) - Vout)/Tl, (7.90)

subject to

(7.91)
260 7 SYSTEM STABILITY

Figure 7.17 Block diagram of a simple controller

Figure 7.18 Modified block diagram of a simple controller of Figure 7.17

The feedback loop can be rearranged to avoid the derivative of Vou, being explicitly
required, as described in section 7.3, and this is shown in Figure 7.18. Equation (7.91)
is now replaced by:
G2
vjb = -Vou, - Va, (7.92)
T2
(
p V a = -VOut
G2
T2
- Va) /T2. (7.93)

Equations (7.90) and (7.93) can be transformed into the algebraic form:

Vout = C1 + Mi (Vin - Vfb), (7.94)

V a = c2 + M2(Vout), (7.95)

where C I , C2, M I and M2 may be determined in the usual way.

A simultaneous solution for the whole system is now possible by combining Equa-
tions (7.92). (7.94)and (7.95) to give:

Vout = C3 + M3(Vin)v (7.96)

 7.7 INTEGRATION 261

where

and

After a solution of VOutis obtained, it may be subjected to the limits of Equation

(7.90). If it is necessary, Vfb can now be evaluated from Equations (7.92) and (7.95)
using the limited value of VOut.
Where this simple controller model represents an AVR, the input Vi, may well be
the (negated) deviation of terminal voltage V , from its specified value (Vs),
It is simple
to treat V , as an extra non-integrable variable rather than incorporate

in the model.
The usual models of speed governors do not have feedback loops associated with
them but the input to the governor (machine speed) is related to the turbine output
(mechanical power) by differential equations. It is, therefore, necessary to solve a set
of simultaneous equations in a similar manner to the example above. The simultaneous
solution should be first made at a point at which limits are applied (i.e. at the valve)
and then, after ensuring the result conforms to the limits, all the other variables around
the loop (including machine speed and rotor angle) can be evaluated.

Solution for saturating AVR exciter Another problem occurs when a non-linear
function is encountered. Equations (7.24) and (7.32) may be combined to form a
single differential equation but this, then, involves a term in Ef which complicates the
evaluation of Ef. As the saturation function is approximate, it may be further simplified
to give:
V e = (k;Ef - k:)Ef, (7.97)
where Ef is the value of Ef at the previous iteration and kT and k; are determined
from E s . The equation describing the saturating exciter is thus:
pEf = (Vam - ( K e - k;)Ef - k?EfEf)/Te, (7.98)
and applying the trapezoidal rule, the algebraic form of solution is:
Ef = [Cef + Mef(Vam)l/El + Mef(k;Ef - k; - Kef)I> (7.99)
where
Cef = (1 - 2(Ke + K e f ) M e f)Ef + M e f Varn,
M , f = h/(2Te + h(Ke + K e f ) ) ,
K e f = klEf - k2.
C c f , M e f and K e f are evaluated once at the beginning of the step and, hence, only
contain information obtained at the end of the previous step.
 Next Page

262 7 SYSTEM STABILITY

Section 1
7Read in steady-state system data

Read in switching data

1
1
Section 2

1. Calculate machine initial conditions

2. Calculate sparse factored inverse of
unfaulted nodal network matrix
I (bifactorization)
&
I
Section 3

1 2. Make switching changes (if necessary)

3. Bifactorization (if necessary)

1
Either Solve for machines and network iteratively
or recalculate nonintegrablevariables
(initially and when switching has occurred)
Section 4
1
1. Print-out results (if necessary)
2. Make power balance check of initial
- conditions (if necessary) J

Section 5

Figure 7.19 Overview of a transient stability program structure