Arrillaga
Arrillaga
and
B. J . HARKER
New Zealand Electricity
A Wiley-Interscience Publication
1. INTRODUCTION . 1
1.1. General background . . 1
1.2. Transmission system development . 2
1.3. Theoretical models and computer programs . 3
4. LOAD FLOW . 58
4.1. Introduction . 58
4.2. Basic nodal method . . 59
4.3. Conditioning of Y matrix . 61
4.4. The case where one voltage is known . 62
4.5. Analytical definition of the problem . 64
4.6. Newton-Raphson method of solving load flows . . 66
Equations relating to power-system load flow
4.7. Techniques which make the Newton-Raphson
method competitive in load flow . 72
Sparsity programming; Triangular factorization; Optimal
ordering; Aids to convergence
4.8. Characteristics of the Newton-Raphson
load flow . . 75
4.9. Decoupled Newton load flow . 76
4.10. Fast decoupled load flow . . 78
4.11. Convergence criteria and tests . 82
4.12. Numerical example . 83
4.13. References . . 88
This chapter deals with the derivation of equivalent circuits of transmission plant
components and with the formation of the system admittance matrix relating the
current and voltage at every node of the transmission system. This constitutes the
basic framework for the algorithms developed in following chapters.
or in matrix form
(v) The actual network admittance matrix which relates the nodal currents to
the voltages by,
p r 1[
circuit tests.
This yields the following matrix equation:
- Y s c + Yoclz
where
ys, is the short-circuit or leakage admittance,
yo, is the open-circuit or magnetizing admittance.
The use of a three-terminal network is restricted to the single-phase
representation and cannot be used as a building block for modelling three-phase
transformer banks.
The magnetizing admittances are usually removed from the transformer model
and added later as small shunt-connected admittances at the transformer
terminals. In the per unit system the model of the single-phase transformer can
then be reduced to a lumped leakage admittance between the primary and
secondary busbars.
A simple equivalent 7c circuit can be deduced from equations (2.3.7) and (2.3.8) the
elements of which can be incorporated into the admittance matrix. This circuit is
illustrated in Fig. 2.4(b)
(a (b)
Fig. 2.4. Transformer with off-nominal tap setting
The equivalent circuit of Fig. 2.4(b) has to be used with care in banks
containing delta-connected windings. In a star-delta bank of single-phase
transformer units, for example, with nominal turns ratio, a value of 1.0 per unit
voltage on each leg of the star winding produces under balanced conditions 1.732
per unit voltage on each leg of the delta winding (rated line to neutral voltage as
base). The structure of the bank requires in the per unit representation an effective
tapping at fi nominal turns ratio on the delta side, i.e. a = 1.732.
For a delta-delta or star-delta transformer with taps on the star winding, the
equivalent circuit of Fig. 2.4(b) would have to be modified to allow for effective
taps to be represented on each side. The equivalent-circuit model of the single-
phase unit can be derived by considering a delta-delta transformer as comprising
a delta-star transformer connected in series (back to back) via a zero-impedance
link to a star-delta transformer, i.e. star windings in series. Both neutrals are
solidly earthed. The leakage impedance of each transformer would be half the
impedance of the equivalent delta-delta transformer. An equivalent per unit
representation of this coupling is shown in Fig. 2.5. Solving this circuit for
terminal currents :
I' (V' - V " ) ,
Ip=-=
a a
Fig. 2.5. Basic equivalent circuit in p.u. for coupling between primary and
secondary coils with both primary and secondary off-nominal tap ratios of a and B
or in matrix form:
These admittance parameters form the primitive network for the coupling
between a primary and secondary coil.
Phase-shifting transformers
To cope with phase-shifting, the transformer of Fig. 2.5 has to be provided with a
complex turns ratio. Moreover, the invariance of the product VI* across the ideal
transformer requires a distinction to be made between the turns ratios for current
and voltage, i.e.
Thus the circuit of Fig. 2.5 has two different turns ratios, i.e.
a, = a +jb for the voltages
and
a,= a -jb for the currents
Solving the modified circuit for terminal currents:
and
and a set of invariant matrices [TI exist. Transformation (2.4.6) will then yield a
diagonal matrix [y],, ,.
In this case the mutually coupled three-phase system has been replaced by
three uncoupled symmetrical systems. In addition, if the generation and loading
are balanced, or may be assumed balanced, then only one system, the positive
sequence system, has any current flow and the other two sequences may be
ignored. This is essentially the situation with the single-phase load flow.
If the original phase admittance matrix [yak] is in its natural unbalanced state
then the transformed admittance matrix [yo, ,] is full. Therefore, current flow of
one sequence will give rise to voltages of all sequences, i.e. the equivalent circuits
for the sequence networks are mutually coupled. In this case the problem of
analysis is no simpler in sequence components than in the original phase
components and symmetrical components should not be used.
From the above considerations it is clear that the asymmetry inherent in all
power systems cannot be studied with any simplification by using the symmetri-
cal component frame of reference. Data in the symmetrical component frame
should only be used when the network element is balanced, for example,
synchronous generators.
In general, however, such an assumption is not valid. Unsymmetrical
interphase coupling exists in transmission lines and to a lesser extent in
transformers and this results in coupling between the sequence networks.
Furthermore, the phase shift introduced by transformer connections is difficult to
represent in sequence component models.
With the use of phase coordinates the following advantages become apparent:
- Any system element maintains its identity.
- Features such as asymmetric impedances, mutual couplings between phases
and between different system elements, and line transpositions are all readily
considered.
- Transformer phase shifts present no problem.
The primitive admittance matrix relates the nodal injected currents to the
branch voltages as follows:
Graphically we represent this partitioning as grouping the six coils into two
compound coils (a) and (b), each composed of three individual admittances. This
is illustrated in Fig. 2.8.
if, and only if yik = yki for i = 1 to 3 and k = 4 to 6. That is, if and only if the
couplings between the two groups of admittances are bilateral.
In this case equation (2.4.9) may be written
(ii)
(iv)
Fig. 2.1 1. Primitive networks and corresponding admittance matrices.
(i) Primitive network using single admittances; (ii) Primitive admittance matrix ;
(iii) Primitive network using compound admittances; (iv) Primitive admittance
matrix
Network subdivision
To enable the transmission system to be modelled in a systematic, logical and
convenient manner the system must be subdivided into more manageable units.
These units, called subsystems, are defined as follows: A SUBSYSTEM is the unit
into which any part of the system may be divided such that no subsystem has any
mutual couplings between its constituent branches and those of the rest of the
system. This definition ensures that the subsystems may be combined in an
extremely straightforward manner.
The system is first subdivided into the most convenient subsystems consistent
with the definition above.
The most convenient unit for a subsystem is a single network element. In
previous sections the nodal admittance matrix representation of all common
elements has been derived.
The subsystem unit is retained for input data organization. The data for any
subsystem is input as a complete unit, the subsystem admittance matrix is
formulated and stored and then all subsystems are combined to form the total
system admittance matrix.
and substituting
1, = I,+ 1, + 1, + I,
Va - Vi = Ia(Ra +joL,) + I,joL,, + I,joL,,
+jwL,;I, -joL,,(I, + I, + I, + I,) + V,
Regrouping and substituting for V,, i.e.
19
+ I,(jwL,, - jwL,, + R, +jwL, -jwL,,)
+ IC(jojL, -jwL,, + R, +jwL, -jwL,)
+ I,O'wL, -jwL,, + R, +jwL, -jwL,,)
+
AVa = Ia(Ra jwL, - 2jwLa, + R, +joL,)
+ Ib(jwLab-jwL,, -jwL,, + R, +jwL,)
+ Ic(jwL, -jwLc, - jwL,, + R, +jwL,)
+ Ig(jwLag-jwL,, -jwL,, + R, +jwL,)
or
A v a = zaa- +
,Ia + Zab-,Ib ZaC-,Ic Zag-,I#+ (2.5.1)
and writting similar equations for the other phases the following matrix equation
results:
From (2.5.3).
From equations (2.5.3) and (2.5.4) and assuming that the ground wire is at zero
potential :
SHUNT ADMI'ITANCE
With reference to Fig. 2.12(b) the potentials of the line conductors are related to
the conductor charges by the matrix e q ~ a t i o n ' ~ ) :
where PA, is a 3 x 3 matrix which includes the effects of the ground wire. The
capacitance matrix of the transmission line of Fig. 2.12 is given by
The series impedance and shunt admittance lumped-.rr model representation of
the three-phase line is shown in Fig. 2.13(a)and its matrix equivalent is illustrated
in Fig. 2.13(b).These two matrices can be represented by compound admittances,
(Fig. 2.13(c)) as described earlier.
t-blf shunt
admittance
Fig. 2.13. Lumped-n model of a short three-phase line series impedance. (a) Full
circuit representation; (b) Matrix equivalent; (c) Using three-phase compound
admittances
Following the rules developed for the formation ofthe admittance matrix using
the compound concept, the nodal injected currents of Fig. 2.13(c) can be related to
the nodal voltages by the equation,
This forms the element admittance matrix representation for the short line
between busbars i and k in terms of 3 x 3 matrix quantities.
This representation may not be accurate enough for electrically long lines. The
physical length at which a line is no longer electrically short depends on the
wavelength, therefore if harmonic frequencies are being considered, this physical
length may be quite small. Using transmission line and wave propagation theory
more exact models may be However, for normal mains frequency
qnalysis, it is considered sufficient to model a long line as a series of two or three
nominal-n sections.
12 x 1 12 x 12 12x 1
It is assumed here that the mutual coupling is bilateral. Therefore.Y2, = Y T2
etc.
The subsystem may be redrawn as Fig. 2.15. The pairs of coupled 3 x 3
compound admittances are now represented as a 6 x 6 compound admittance.
The matrix representation is also shown. Following this representation and the
labelling of the admittance blocks in the figure, the admittance matrix may be
written in terms of the 6 x 6 compound coils as,
(ii)
Fig. 2.15. 6 x 6 compound admittance representation of two coupled three-
phase lines. (i) 6 x 6 Matrix representation; (ii) 6 x 6 Compound admittance
representation
This is clearly identical to equation (2.5.10) with the appropriate matrix
partitioning.
The representation of Fig. 2.15 is more concise and the formation of equation
(2.5.11 ) from this representation is straightforward, being exactly similar to that
which results from the use of 3 x 3 compound admittances for the normal single
three-phase line.
The data which must be available, to enable coupled lines to be treated in a
similar manner to single lines, is the series impedance and shunt admittance
matrices. These matrices are of order 3 x 3 for a single line, 6 x 6 for two coupled
lines, 9 x 9 for three and 12 x 12 for four coupled lines.
Once the matrices [Z,] and [Y,] are available, the admittance matrix for the
subsystem is formed by application of equation (2.5.11).
When all the busbars of the coupled lines are distinct, the sub-system may be
combined directly into the system admittance matrix. However, if the busbars are
not distinct then the admittance matrix as derived from equation (2.5.11)must be
modified. This is considered in the following section.
The admittance matrix derived previously related the currents and voltages at
the four busbar A l , A2, B1 and B2. This relationship is given by:
25
The nodal injected current at busbar A,(',), is given by,
A
' = 'A1 + 'A2
similarly
B
' =IBl +' ~ 2
Also from inspection of Fig. 2.16,
VA= VAl = VA2
VB= VB1 = VB1
The required matrix equation relates the nodal injected currents, 1, and I,, to
the voltages at these busbar. This is readily derived from equation (2.5.12)and the
conditions specified above. This is simply a matter of adding appropriate rows
and columns and yields,
This matrix [YAB] is the required nodal admittance matrix for the subsystem.
It should be noted that the matrix in equation (2.5.12) must be retained as it is
required in the calculation of the individual line power flows.
Shunt elements t
Shunt reactors and capacitors are used in a power system for reactive power
control. The data for these elements are usually given in terms of their rated MVA
and rated kV; the equivalent phase admittance in p.u. is calculated from these
data.
Consider, as an example, a three-phase capacitor bank shown in Fig. 2.17. A
similar triple representation as that for a line section is illustrated. The final two
forms are the most compact and will be used exclusively from this point on.
Series elements
Any element connected directly between two buses may be considered a series
element. Series elements are often taken as being a section in a line sectional-
ization which is described later in the chapter.
A typical example is the series capacitor bank which is usually taken as
uncoupled, i.e. the admittance matrix is diagonal.
This can be represented graphically as in Fig. 2.18.
The admittance matrix for the subsystem can be written by inspection as:
1- S E I I
Fig. 2.1 8. Graphic representation of series capacitor bank between nodes i and k
where
yk is the mutual admittance between primary coils;
yk is the mutual admittance between primary and secondary coils on
different cores ;
yz is the mutual admittance between secondary coils.
For three separate single-phase units all the primed values are effectively zero.
In three-phase units the primed values, representing parasitic interphase
coupling, do have a noticeable effect. This effect can be interpreted through the
symmetrical component equivalent circuits.
Let us now consider the Wye G-Delta connection illustrated in Fig. 2.23.
Moreover, if the primitive admittances are expressed in per unit, with both the
primary and secondary voltages being one per unit, the Wye-Delta transformer
. model must include an effective turns ratio offi. The upper right and lower left
quadrants of matrix (2.6.10)must be divided by fi and the lower right quadrant
by 3.
In the particular case of three-single phase transformer units connected in Wye
G-Delta all the y' and y" terms will disappear. Ignoring off nominal taps (but
keeping in mind the effective fi turns ratio in per unit) the nodal admittance
matrix equation relating the nodal currents to the nodal voltages is
The models for the other common connections can be derived following a
similar procedure.
In general, any two-winding three-phase transformer may be represented using
two coupled compound coils. The network and admittance matrix for this
representation is illustrated in Fig. 2.25.
Table 2.1. Characteristic submatrices used in forming the transformer admittance matrices
where a,, a, and a, are the off-nominal taps on windings 1,2, and 3 respectively.
In addition any windings connected in delta will, because of the per unit system,
have an effective tap offi.
The nodal admittance matrix for the transformer windings is:
where
36
Therefore
Secondary side
The delta connection on the secondary side introduces an effective ,,hturns ratio
and the sequence components admittance matrix is
Mutual terms
The mutual admittance submatrix of equation (2.6.10), modified for effective
turns ratio is transformed as follows:
37
Recombining the sequence components submatrices yields
'sc
Wye G w P
J'-
Wye G Delta
W Y ~ W Y ~
w e blta
Delta Delta
Line sectionalization
A line may be divided into sections to account for features such as the following:
- Transposition of line conductors.
- Change of type of supporting towers.
- Variation of soil permitivity.
39
- Improvement of line representation. (Series of two or more equivalent n
networks.)
- Series capacitors for line compensation.
- Lumping of series elements not central to a particular study.
An example of a line divided into a number of sections is shown in Fig. 2.29.
Section 4 Section6
Bus A PI P2
Section I Section Section
1 : 2 j 3
r
I I
I
1
0
11 II 4
I\
b
1 II ti
I
Transposition /
Change of
Series &aci+ors
obc obc
Phases configuration
FEATURES OF INTEREST
2 2 0 kV
22O/66 kV
Bus B
Bus A
2 2 0 kV
I I
I
Section No1 ! ~ e c t i o n ~ o 2Section
! No 3 1
( ii 1
Fig. 2.31. Sample system to illustrate line sectionalization.
(i) System single line diagram;(ii) system redrawn to illustrate line
sectionalization
Table 2 . 3 . ABCD Parameter matrices for the common section types
Shunt element
I [.I [OI I
Series element
i ;; I-;;J
In Table 2.3 [u] is the unit matrix, [O] is a matrix of zeros, and all other matrices
have been defined in their respective sections.
Note, all the above matrices have dimensions corresponding to the number of
coupled three-phase elements in the section.
Once the resultant ABCD parameters have been found the equivalent nodal
admittance matrix for the subsystem can be calculated from the following
equation.
CYl =
I [Dl CBI - '
CBI - '
CCI - [Dl CBI - ' CAI
'
- CBI - CAI
2.8. REFERENCES
1. E. Clarke, 1943. Circuit Analysis of A.C. Power Systems, Vol. I. Wiley and Sons Ltd.,
New York.
2. G. Kron, republished in 1965. Tensor Analysis of Networks. MacDonald, London.
3. M. S. Chen and W. E. Dillon, July 1974. 'Power-system modelling', Proc. IEEE, 62 (7),
901.
4. M. A. Laughton, 1968. 'Analysis of unbalanced polyphase networks by the method of
phase coordinates, Part I. System representation in phase frame of reference', Proc.
IEE, 115 (8), 1 163- 1 172.
5. K. J. Bowman and J. M. McNamee, 1964. 'Development of equivalent pi and T matrix
circuits for long untransposed transmission lines', IEEE, PAS-84, 625-632.
6. L. M. Wedepohl and R. G. Wasley, 1966. 'Wave propagation in multiconductor
overhead lines', Proc. IEE, 113 (4), 627-632.
7. W. E. Dillon and M. S. Chen, 1972. 'Transformer modelling in unbalanced three-phase .
networks.' IEEE Summer Power Meeting, Vancouver.
Modelling of Static A C - D . C .
Conversion Plant
3.1. INTRODUCTION
Although the power electronic device is basically a switch, it is only explicitly
represented as such in dynamic studies (Chapters 11 and 12). The periodicity of
switching sequences can be used in steady-state studies to model the active and
reactive power loading conditions of a.c.-d.c. converters at the relevant busbars.
Such modelling is discussed here with reference to the most common con-
figuration used in power systems, i.e. the three-phase bridge rectifier shown in
Fig. 3.1.
For large power ratings static converter units generally consist of a number of
series and/or parallel connected bridges, some or all bridges being phase-shifted
relative to the others. With these configurations twelve-pulse and higher pulse
numbers can be achieved to reduce the distortion of the supply current with
limited or no filtering. A multiple bridge rectifier can therefore be modelled as a
single equivalent bridge with a sinusoidal supply voltage at the terminals.
The following basic assumptions are normally made in the development of the
steady-state m ~ d e l ( ' ) . ( ~ ) . ( ~ ) . ( ~ ) :
(i) The forward voltage drop in a conducting valve is neglected so that the
valve may be considered as a switch. This is justified by the fact that the
voltage drop is very small in comparison with the normal operating
\
voltage. It is, further, quite independent of the current and should,
therefore, play an insignificant part in the commutation process since all
valves commutating on the same side of the bridge suffer similar drops.
Such a voltage drop is taken into account by adding it to the d.c. line
resistance. The transformer windings resistance is also ignored in the
development of the equations, though it should also be included to
calculate the power loss.
(ii) The converter transformer leakage reactances as viewed from the secon-
dary terminals are identical for the three phases, and variations of leakage
reactance caused by on-load tap-changing are ignored.
(iii) The direct current ripple is ignored, i.e. sufficient smoothing inductance is
assumed on the d.c. side.
3.2. RECTIFICATION
Rectifier loads can use diode and thyristor elements in full or half-bridge
configurations. In some cases the diode bridges are complemented by transfor-
mer on-load tap-changer and saturable reactor control. Saturable reactors
produce the same effect as thyristor control over a limited range of delay angles
v = a v sin (3.2.1)
43
0
With delay angle control the average rectified voltage (shown in Fig. 3.3) is thus
In practice the voltage waveform is that of Fig. 3.4, where a voltage area (6A)is
lost due to the reactance (X,) of the a.c. system (as seen from the converter),
Fig. 3.4. Effect of commutation reactance. (a) Alternating cur-
rent; (b) D.c. voltage wavefoms
where e,, e, are the instantaneous voltages of phases a and b respectively, and i, is
the incoming valve (commutating) current.
It must be emphasized that the commutating voltage (Vterm) is the a.c. voltage
at the closest point to the converter bridge where sinusoidal waveforms can be
assumed. The commutation reactance (X,) is the reactance between the point at
which Vtermexists and the bridge. Where filters are installed the filter busbar
voltage can be used as Vterm.In the absence of filters, Vtermmust be established at
some remote point and X, must be modified to include the system impedance
from the remote point to the converter.
With perfect filtering, only the fundamental component of the current
waveform will appear in the a.c. system. This component is obtained from the
Fourier analysis of the current waveform in Fig. 3.4, and requires information
of i, and u.
Taking as a reference the instant when the line voltage (e, - e,) is zero, equation
(3.2.3) can be written as
X, di
$a. Vtermsin o t = 2 -2
o dt
and integrating with respect to (ot):
1
- -a. Vtermcos ot + K = X,i,
fi
From the initial condition: i, = 0 at ot = a the following expressions for K and
i, are obtained.
where
,/ [cos 2a - cos 2(a + u)I2 + [2u + sin 2a - sin 2(a + u)]
k=
~ [ C O Sa - cos(a + u)]
When using per unit values based on a common power and voltage base on both
sides of the converter, the direct current base has to be fi
times larger than the
a.c. current base (as explained in Section 6.3) and equation (3.2.10) becomes
Using the fundamental components of voltage and current and assuming perfect
filtering at the converter terminals the power factor angle at the converter
terminals is 4 (the displacement between fundamental voltage and current
waveforms) and we may write
and
3.3. INVERSION
Owing to the unidirectional nature of current flow through the converter valves,
power reversal (i.e. power flow from the d.c. to the a.c. side) requires direct voltage
polarity reversal. This is achieved by delay angle control, which, in the absence of
commutation overlap produces rectification between 0" < a < 90" and inversion
between 90" < a < 180". In the presence of overlap, the value of 'a' at which
inversion begins is always less than 90". Moreover, unlike with rectification, full
inversion (i.e. a = 180")can not be achieved in practice. This is due to the existence
of a certain deionization angle y at the end of the conducting period, before the
voltage across the commutating valve reverses, i.e.
If the above condition is not met (yo being the minimum required extinction
angle) a commutation failure occurs; this event would upset the normal
conducting sequence and preclude the use of the steady-state model derived in
this chapter.
The inverter voltage, although of opposite polarity with respect to the rectifier,
is usually expressed as positive when considered in isolation.
Typical inverter voltage and current waveforms are illustrated in Fig. 3.5. By
similarity with the waveforms of Fig. 3.4, the following expression can be written
for the inverter voltage in terms of the extinction angle:
,~, ,Q,~~=
,~i f ,~, ,4,~~=
sin ,~i f sin (n - 4 ) (3.3.3)
Fig. 3.6. P and Q vector diagram
Equations (3.3.2) and (3.3.3) indicate that the inversion process still requires
reactive power supply from the a.c. side. The vector diagram of Fig. 3.6 illustrates
the sign of P and Q for rectification and inversion.
Fig. 3.7. 'n' bridges connected in series on the d.c. side and in
parallel on the a.c. side
where j represents any of the n bridges. For bridges connected in parallel on the
d.c. side the equivalent bridge commutation reactance is:
It should be noted that with perfect filtering or when many bridges are used
with different transformer phase shifts the voltage on the a.c. side of the converter
transformers may be assumed to be sinusoidal and hence Xs, has no influence on
the commutation.
Moreover, the presence of local plant components at the converter terminals
may affect the commutation reactance. By way of example, let us consider the two
ends of the New Zealand h.v.d.c. link (with reference to Fig. 3.8). It must be noted
that h.v.d.c. schemes are normally designed for twelve-pulse operation and that
filters are always provided (i.e. the system impedance can be ignored).
AC Benrnore Haywards
220 kV
To South
I
1 Pole I
t
To North
Sea + Islandsystem
Synchronous
cornpensotors
where
X,-transformer secondary leakage reactance
X,-transformer primary leakage reactance
X,-transformer tertiary leakage reactance
X&' is the subtransient reactance of the synchroncus condenser unit
Fig. 3.9. Equivalent circuit for the calculation of the com-
mutation reactance at Haywards end
1
A. C.
t+ 4 system
I
where
X-two-winding transformer leakage reactance
X,-interconnecting transformer secondary leakage reactance (note
filters connected to tertiary winding)
X:-generator subtransient reactance
n- number of generators connected
where Rd is the resistance of the link and includes the loop transmission resistance
(if any), the resistance of the smoothing reactors and the converter valves.
Convertor I Convertor I1
lnvertor c.c., I
- C
In cases where the power rating of the d.c. link is comparable with the rating of
either the sending or receiving a.c. systems interconnected by the link, the
frequency of the smaller a.c. system is often controlled to a large extent by the d.c.
link. With power-frequency (P.F.) control if the frequency goes out of pre-
specified limits, the output power is made proportional to the deviation of
frequency from its nominal value. Frequency control is analogous to the current
control described earlier, i.e. the converter with lower voltage determines the
direct voltage of the line and the one with higher voltage determines the
frequency. Again, current limits have to be imposed, which override the frequency
error signal.
C.P.1E.A. and C.C.1E.A. controls were evolved principally for bulk point-to-
point power transmission over long distances or submarine crossings and are still
the main control modes in present use.
Multiterminal d.c. schemes are also being considered, based on the basic
controls already described. Two alternatives are possible, i.e. constant voltage
parallel"' and constant current series(@schemes.
3.6. REFERENCES
1. C. Adamson and N. G. Hingorani, 1960.H.VD.C. Power Transmission. Garraway Ltd.,
London.
2. E. W. Kimbark, 1971. Direct Current Transmission.Vol. I . Wiley Interscience, London.
3. B. J. Cory (ed.), 1965. High-Voltage Direct-Current Converters and Systems. Macdonald,
London.
4. E. Uhlmann, 1975. Power Transmission by Direct Current. Springer-Verlag,
Berlin/Heidelberg.
5. U. Lamm and E. Uhlmann and P. Danfors, 1963. 'Some aspects of tapping of h.v.d.c.
transmission systems', Direct Current, 8, 124-1 29.
6. C. Adamson and J. Arrillaga, 1968. 'behaviour of multiterminal a.c.4.c. intercon-
nections with series-connected stations', Proc. LEE, 115 (1 l), 1685-1692.
Load Flow
4.1. INTRODUCTION
Under normal conditions electrical transmission systems operate in their steady-
state mode and the basic calculation required to determine the characteristics of
this state is termed load flow (or power flow).
The object of load-flow calculations is to determine the steady-state operating
characteristics of the power generation/transmission system for a given set of
busbar loads. Active power generation is normally specified according to
economic-dispatching practice and the generator voltage magnitude is normally
maintained at a specified level by the automatic voltage regulator acting on the
machine excitation. Loads are normally specified by their constant active and
reactive power requirement, assumed unaffected by the small variations of
voltage and frequency expected during normal steady-state operation.
The solution is expected to provide information of voltage magnitudes and
angles, active and reactive power flows in the individual transmission units, losses
and the reactive power generated or absorbed at voltage-controlled buses.
The load flow problem is formulated in its basic analytical form in this chapter
with the network represented by linear, bilateral and balanced lumped para-
meters. However the power and voltage constraints make the problem non-linear
and the numerical solution must therefore be iterative in nature.
The currents problem faced in the development of load flow are: an ever
increasing size of systems to be solved, on-line applications for automatic control,
and system optimization. Hundreds of contributions have been offered in the
literature to overcome these problems.(')
Five main properties are required of a load-flow solution method.
(i) High computational speed. This is especially important when dealing with
large systems, real time applications (on-line), multiple case load-flow
such as in system security assessment, and also in inter-active
applications.
(ii) Low computer storage. This is important for large systems and in the use
of computers with small core storage availability, e.g. mini-computers for
on-line application.
(iii) Reliability of solution. It is necessary that a solution be obtained for
illconditioned problems, in outage studies and for real time applications.
(iv) Versatility. An ability on the part of load-flow to handle conventional and
special features (e.g. the adjustment of tap ratios on transformers; different
representations of power-system apparatus), and its suitability for
incorporation into more complicated processes.
(v) Simplicity. The ease of coding a computer program of the load-flow
algorithm.
The type of solution required from a load-flow also determines the method
used :
accurate or approximate
unadjusted or adjusted
off-line or on-line
single case or multiple cases
The first column are requirements needed for considering optimal load-flow and
stability studies, and the second column those needed for assessing security of a
system. Obviously, solutions may have a mixture of the properties from either
column.
The first practical digital solution methods for load flow were the Y matrix-
iterative methods.(2)These were suitable because of the low storage requirements,
but had the disadvantage of converging slowly or not at all. Z matrix methodd3)
were developed which overcame the reliability problem but a sacrifice was made
of storage and speed with large systems.
The Newton-Raphson m e t h ~ d ( ~ )was
. ( ~developed
) at this time and was found
to have very strong convergence. It was not, however, made competitive until
sparsity programming and optimally ordered Gaus~ian-elimination(~)~(~)~(*) were
introduced, which reduced both storage and solution time.
Nonlinear programming and hybrid methods have also been developed, but
these have created only academic interest and have not been accepted by
industrial users of load flow. The Newton-Raphson method and techniques
derived from this algorithm, satisfy the requirements of solution type and
programming properties better than previously used techniques and are
gradually replacing them.
.,
Fig. 4.1. Simple network showing nodal quantities
Figure 4.1 gives a simple network showing the nodal currents, voltages and
powers.
In the nodal method, it is convenient to use branch admittances rather than
impedances. Denoting the voltages of nodes k and i as Ekand Ei respectively, and
the admittance of the branch between them as yki,then the current flowing in this
branch from node k to node i is given by:
Let the nodes in the network be numbered 0, 1 . . .n, where 0 designates the
reference node (ground). By Kirchhoffs current law, the injected current I, must
be equal to the sum of the currents leaving node k, hence:
If this equation is written for all the nodes except the reference, i.e. for all busbar
in the case of a power-system network, then a complete set of equations defining
the network is obtained in matrix form as:
where
n
Suppose that the injected currents are known, and nodal voltages are
unknown. In this case no solution for the latter is possible. The Y matrix is
described as being singular, i.e. it has no inverse, and this is easily detected in this
example by noting that the sum of the elements in each row and column is zero,
which is a sufficient condition for singularity, mathematically speaking. Hence, if
it is not possible to express the voltages in the form E = Y -'. I, it is clearly
impossible to solve equation (4.3.1.)by any method, whether involving inversion
of Y or otherwise.
The reason for this is obvious: we are attempting to solve a network whose
reference node is disconnected from the remainder, i.e. there is no effective
reference node, and an infinite number of voltage solutions will satisfy the given
injected current values.
When, however, a shunt admittance from at least one of the busbar in the
network of Fig. 4.2. is present, the problem of insolubility immediately vanishes in
theory, but not necessarily in practice. Practical computation cannot be
performed with absolute accuracy, and during a sequence of arithmetic
operations, rounding errors due to working with a finite number of decimal
places accumulate. If the problem is well-conditioned and the numerical solution
technique is suitable, these errors remain small and do not mask the eventual
results. If the problem is ill-conditioned, and this usually depends upon the
properties of the system being analysed, any computational errors introduced are
likely to become large with respect to the true solution.
It is easy to see intuitively that if a network having zero shunt admittances
cannot be solved even when working with absolute computational accuracy, then
a network having very small shunt admittances may well present dficulties when
working with limited computational accuracy. This reasoning provides a key to
the practical problems of network, i.e. Y matrix, conditioning. A network with
shunt admittances which are small with respect to the other branch admittances
is likely to be ill-conditioned, and the conditioning tends to improve with the size
of the shunt admittances, i.e. with the electrical connection between the network
busbars and the reference node.
The first row of this set may now be eliminated, leaving (n - 1)equations in (n - 1)
unknowns, E 2 . .. En. In matrix form, this becomes:
The new matrix Y is obtained from the full admittance matrix Y merely by
removing the row and column corresponding to the fixed-voltage busbar, both in
the present case where it is numbered 1 or in general.
In summation notation, the new equations are
n
1,- y , , ~ =
, 1 YkiEi for k = 2 ...n (4.4.5)
i=2
the two unknown variables at each node in a system. A second set of variable
equations, which are linear, are derived from the first set, and an iteration method
is applied to this second set.
The basic algorithm which load-flow programs use is depicted in Fig. 4.3.
System data, such as busbar power conditions, network connections and
impedance, are read in and the admittance matrix formed. Initial voltages are
+
specified to all buses; for base case load flows, P, Q buses are set to 1 j O while
P, V busbars are set to V + jO.
The iteration cycle is terminated when the busbar voltages and angles are such
that the specified conditions of load and generation are satisfied. This condition is
accepted when power mismatches for all buses are less than a small tolerance, q,,
or voltage increments less than q2. Typical figures for q, and q2 are 0.01 p.u. and
0.001 p.u. respectively. The sum of the square of the absolute values of power
mismatches is a further criterion sometimes used.
When a solution has been reached, complete terminal c lditions for all buses are
computed. Line power flows and losses and system totals can then be calculated.
66
4.6. NEWTON-RAPHSON METHOD OF SOLVING LOAD mows
The generalized Newton-Raphson method is an iterative algorithm for solving a
set of simultaneous nonlinear equations in an equal number of unknowns.
f,(x,)=O for k = 1 + N (4.6.1)
m=l+N
At each iteration of the N - R method, the nonlinear problem is approximated
by the linear matrix equation. The linearizing approximation can best be
visualized in the case of a single-variable problem.
In Fig. 4.4. xPis an approximation to the solution, with error AxP at iteration p.
Then
If the initial estimates of the variable xP is near the solution value, AxP will be
relatively small and all terms of higher powers can be neglected. Hence,
+
f ( x p ) AxP f ' ( x P =
) 0 (4.6.4)
and represents the slopes of the tangent hyperplanes which approximate the
functions fk(xm)at each iteration point.
The Newton-Raphson algorithm will converge quadratically if the functions
have continuous first derivates in the neighbourhood of the solution, the
Jacobian matrix is nonsingular, and the initial approximations of x are close to
the actual solutions. However the method is sensitive to the behaviours of the
functions fk(xm)and hence to their formulation. The more linear they are, the
more rapidly and reliably Newton's method converges. Nonsmoothness, i.e.
humps, in any one of the functions in the region of interest, can cause convergence
delays, total failure or misdirection to a nonuseful solution.
where I, is the current injected into a bus k. The power at a bus is then given by
In polar coordinates the real and imaginary parts of equation (4.6.10) are
Pk= x
me k
VkVm(Gkm
cos 8 , +B, sin 8,,) (4.6.11)
Linear relationships are obtained for small variations in the variables 0 and V by
forming the total differentials, the resulting equations being:
For a PQ busbar:
APk = c ap
~
m s k aOm
0 + ,LA
ap
1 vm
m s k avm
and
For a PV busbar:
Only equation (4.6.13)is used, since Qk is not specified.
For a slack busbar:
No equations.
The voltage magnitudes appearing in equations (4.6.13)and (4.6.14)for PV and
slack busbars are not variables, but are fixed at their specified values. Similarly 0
at the slack busbar is fixed.
The complete set of defining equations is made up of two for each PQ busbar
and one for each PV busbar. The problem variables are V and 0 for each PQ
busbar and 0 for each PV busbar. The number of variables is therefore equal to
the number of equations. Algorithm (4.6.7) then becomes:
P mismatches 9 corrections
for all P Q for all PQ and
and PV busbars PV busbars
Q mismatches AVp V corrections
for all PQ busbars
The division of each AVP by Vf -'does not numerically affect the algorithm, but
simplifies some of the Jacobian-matrix terms. For busbars k and rn (not row k and
column rn in the matrix)
apk = VkVm(Gkm
Hkm= -- sin Ohm - B , cos 0,)
a0m
apk
cos ohm+ Bkmsin Ohm)
Nkm= Vm-= VkVm(Gkm
av m
aQk
Lkm= Vm- = VkVm(Gkm
sin 0 , - B, cos Ohm)
av m
and for m = k:
At a voltage controlled bus the voltage magnitude is fixed but not the phase angle.
Hence both ek and f k vary at each iteration. It is necessary to provide another
equation
V,"=e,t+f,Z
to be solved with the real power equation for these buses.
Linear relationships are obtained for small variations in e andf , by forming the
total differentials,
apk
AP,= C ---At?,+ C -apk
Afm
m s k aem msk ' f m
= C
msk
U k m Aem +C
msk
Wkm A f m
70
for all nonvoltage controlled buses;
and the values of the partial differentials, which are the Jacobian elements, are
given by
Wkk = a, - G k k e k - B k k f k
Tkk = bk - B k k e k +Gkkfk
Ukk = - bk - B k k e k +Gkkfk
E E , = 2ek
FFk = 2 f k
For voltage controlled buses, V is specified, but not the real and imaginary
components of voltage, e and f . Approximations can be made, for example, by
ignoring the off-diagonal elements in the Jacobian matrix, as the diagonal
elements are the largest. Alternatively for the calculation of the elements the
voltages can be considered as E = 1 +jO. The off-diagonal elements then become
constant.
The polar coordinate representation appears to have computational advan-
tages over rectangular coordinates. Real power mismatch equations are present
for all buses except the slack bus, while reactive power mismatch equations are
needed for nonvoltage controlled buses only.
The Jacobian matrix has the sparsity of the admittance matrix [Y] and has
positional but not numerical symmetry. To gain in computation, the form of
[A& AV/ VIT is normally used for the variable voltage vector. Both increments
are dimensionless and the Jacobian coefficients are now symmetric in structure
though not in value. The values of [J] are all functions of the voltage variables V
and 8 and must be recalculated each iteration.
As an example, the Jacobian matrix equation for the four-busbar system of
Fig. 4.5 is given as equation (4.6.17)
QSlack
1,
I
Fig. 4.5. Sample system
results
Sparsity programming
In conventional matrix programming, double subscript arrays are used for the
location of elements. With sparsity programming,(6)only the nonzero elements
are stored, in one or more vectors, plus integer vectors for identification.
For the admittance matrix of order n the conventional storage requirements
are nZwords, but by sparsity programming 66 + 3n words are required, where b is
the number of branches in the system. Typically b = lSn, and the total storage is
12n words. For a large system (say 500 buses) the ratio of storage requirements of
conventional and sparse techniques is about 40: 1.
Triangular factorization
To solve the Jacobian matrix equation (4.6.15), represented here as
for increments in voltage, the direct method is to find the inverse of [JI and solve
for [AEJ from
CAE3 = [A - [Asl' (4.7.1)
In power systems [J] is usually sparse but [JI-' is a full matrix.
The method of triangular factorization solves for the vector [AEJ by
eliminating [J] to an upper triangular matrix with a leading diagonal, and then
back-substituting for [ A q , i.e. eliminate to
CAS'l =C u l CAE3
and back-substitute
[U] - [AS'] = [AEI
The triangulation of the Jacobian is best done by rows. Those rows below the
one being operated on need not be entered until required. This means that the
maximum storage is that of the resultant upper triangle and diagonal. The lower
triangle can then be used to record operations.
The number of multiplications and additions to triangulate a full matrix is $N3,
compared to N3 to find the inverse. With sparsity programming the number of
operations varies as a factor of N. If rows are normalized N further operations are
saved.
Optimal ordering
In power-system load flow, the Jacobian matrix is usually diagonally dominant
which implies small round-off errors in computation. When a sparse matrix is
triangulated, nonzero terms are added in the upper triangle. The number added
is affected by the order of the row eliminations, and total computation time
increases with more terms.
The pivot element is selected to minimize the accumulation of nonzero terms,
and hence conserve sparsity, rather than minimizing round off error. The
diagonals are used as pivots.
Optimal ordering of row eliminations to con ,sparsity is a practical
impossibility due to the complexity of programming m d time involved. However,
semioptimal schemes are used and these can be divided into two sections.
(a) Pre~rdering.'~'
Nodes are renumbered before triangulation. No complicated
programming or storage is required to keep track of row and column
interchanges.
(i) Nodes are numbered in sequence of increasing number of connected lines.
(ii) Diagonal banding-nonzero elements are arranged about either the major
or minor diagonals of the matrix.
(b) Dynamic ordering.") Ordering is effected at each row during the elimination.
(i) At each step in the elimination, the next row to be operated on is that with
the fewest nonzero terms.
(ii) At each step in the elimination, the next row to be operated on is that
which introduces the fewest new nonzero terms, one step ahead.
(iii) At each step in the elimination, the next row to be operated on is that
which introduces the fewest new nonzero terms, two steps ahead. This may
be extended to the fully optimal case of looking at the effect in the final
step.
(iv) With cluster ordering, the network is subdivided into groups which are
then optimally ordered. This is most efficient if the groups have a
minimum of physical intertie. The matrix is then anchor banded.
The best method arises from a trade-off between a processing sequence which
requires the least number of operations, and time and memory requirements.
The dynamic ordering scheme of choosing the next row to be eliminated as that
with the fewest nonzero terms, appears to be better than all other schemes in
sparsity conservation, number of arithmetic operations required, ordering times
and total solution time.
However, there are conditions under which other ordering would be prefer-
able, e.g. with system changes affecting only a few rows these rows should be
numbered last; when the subnetworks have relatively few interconnections it is
better to use cluster ordering.
Aids to convergence
The N-R method can diverge very rapidly or converge to the wrong solution if
the equations are not well behaved or if the starting voltages are badly chosen.
Such problems can often be overcome by a variety of techniques. The simplest
device is to impose a limit on the size of each A8 and AV correction at each
iteration. Figure 4.7. illustrates a case which would diverge without this device.
Another more complicated method is to calculate good starting values for the
8's and V's, which also reduces the number of iterations required.
In power-system load flow, setting voltage controlled buses to V+jO and
nonvoltage controlled buses to 1 + j O may give a poor starting point for the
Newton-Raphson method.
If previously stored solutions for a network are available these should be used.
One or two iterations of a Y matrix iterative method(2)can be applied before
Fig. 4.7. Example of diverging solution
commencing the Newton method. This shows fast initial convergence unless the
problem is ill-conditioned, in which case divergence occurs.
A more reliable method is the use of one iteration of a d.c. load flow (i.e.
neglecting losses and reactive power conditions) to provide estimates of voltage
angles; followed by one iteration of a similar type of direct solution to obtain
voltage magnitudes. The total computing time for both sets of equations is about
50 percent of one Newton-Raphson iteration and the extra storage required is
only in the programming statements. The resulting combined algorithm is faster
and more reliable than the formal Newton method and can be used to monitor
diverging or difficult cases, before commencing the Newton-Raphson algorithm.
where for the reference node 80 = 0 and Vk= Vo. The values of Sk and 22,
represent real and reactive power quantities respectively and [7'l and [U]are
given by
where Zkmand X,, are the branch impedance and reactance respectively. [U] is
constant valued and needs be triangulated once only for a solution. [7-'J is
recalculated and triangulated each iteration.
The two equations (4.9.1)and (4.9.2)are solved alternately until a solution is
obtained. These equations can be solved using Newton's method, by expressing
the Jacobian equations as
where
CAP] = CAP]
If the submatrices N and J are neglected, since they represented the weak
coupling between 'P-8' and 'Q-V', the following decoupled equations result
CAP] = [HI [A81 (4.9.11)
It has been found that the latter equation is relatively unstable at some distance
from the exact solution due to the nonlinear defining functions. An improvement
in convergence is obtained by replacing this with the polar current-mismatch
formulation(7)
CAI] = [Dl [A vl (4.9.13)
Alternatively the right-hand side of both equations (4.9.11) and (4.9.12)is divided
by voltage magnitude V
[APIvI = [ A ] [A81 (4.9.14)
The equations are solved successively using the most up-to-date values of V and 8
available. [ A ] and [ C ] are sparse, nonsymmetric in value and are both functions
of V and 8.
They must be calculated and triangulated each iteration.
Further approximations that can be made are to assume that E , = 1.0 p.u., for
all buses, and G,, < B,, in calculating the Jacobian elements. The off-diagonal
terms then become symmetric about the leading diagonal.
The decoupled Newton method compares very favourably with the formal
Newton method. While reliability is just as high for ill-conditioned problems, the
decoupled method is simple and computationally efficient. Storage of the
Jacobian and matrix triangulation is saved by a factor of 4, or an overall saving of
30-40 percent on the formal Newton load flow. Computation time per iteration is
also less than the Newton method.
However, the convergence characteristics of the decoupled method are linear,
the quadratic characteristics of the formal Newton being sacrified. Thus, for high
accuracies, more iterations are required. This is offset for practical accuracies by
the fast initial convergence of the method. Typically, voltage magnitudes
converge to within 0.3 percent of the final solution on the first iteration and may
be used as a check for instability. Phase angles converge more slowly than voltage
magnitudes but the overall solution is reached in 2 + 5 iterations. Adjusted
solutions (the inclusion of transformer taps, phase shifters, interarea power
transfers, Q and V limits) take many more iterations.
The Newton methods can be expressed as follows(12):
where
E= 1 for the full Newton-Raphson method
E =0 for the decoupled Newton algorithm
A Taylor series expansion of the Jacobian about E = 0 results in a first-order
approximation of the Newton-Raphson method whereas the decoupled method
is a zero order approximation. The method has quadratic convergence properties
because of the coupling, but retains storage requirements similar to that of the
decoupled method.
Qk = vk 1
msk
vm(Gkm sin O k m - Bkm cOs e k m ) (4.10.2)
where Okm = 8, - 8,.
A decoupled method which directly relates powers and voltages is derived
using the series approximations for the trigonometric terms in equations (4.10.1)
and (4.10.2).
The equations, over all buses, can be expressed in their simplified matrix form
where P and Q are terms of real and reactive power respectively and
'kk= 1
mok
'kmBkm
and Bkmare the imaginary parts of the admittance matrix. To simplify still further,
line resistances may be neglected in the calculation of elements of [B].
An improvement over equations (4.10.5) and (4.10.6) is based on the decoupled
equations (4.9.14) and (4.9.15) which have less nonlinear defining functions.
Applying the same assumptions listed previously, we obtain the equations
[APIV] = [B*] [A61 (4.10.7)
The equations are solved alternatively using the most recent values of V and 0
available as shown in Fig. 4.8.(*)
, converged
NO
converged
OUTPUT RESULTS
The matrices Bf and B" are real and are of order ( N - 1) and ( N - M )
respectively, where N is the number of busbars and M is the number of PV
busbars. B" is symmetric in value and so is B' if phase shifters are ignored; it is
found that the performance of the algorithm is not adversely affected. The
elements of the matrices are constant and need to be evaluated and triangulated
only once for a network.
Convergence is geometric, two to five iterations are required for practical
accuracies, and more reliable than the formal Newton's method. This is because
the elements of B' and B" are fixed approximations to the tangents of the defining
functions APIV and AQ/V, and are not susceptible to any 'humps' in the defining
functions.
If AP/V and AQ/V are computed efficiently, then the speed for iterations of the
fast decoupled method is about five times that of the formal Newton-Raphson or
about two-thirds that of the Gauss-Seidel method. Storage requirements are
about 60 percent of the formal Newton, but slightly more than the decoupled
Newton method.
Changes in system configurations are easily effected, and while adjusted
solutions take many more iterations these are short in time and the overall
solution time is still low.
The fast decoupled Newton load flow can be used in optimization studies for
a network and is particularly useful for accurate information of both real and
reactive power for multiple load-flow studies, as in contingency evaluation for
system security assessment.
00000000000000000
OOOOOGOOOOOOOOOOO
C h UC
¶
I-
I
0
w h: r 0 w
(roo' 0 (D (rcn
00000000C000Ul000r I C
e................ <
00000000000000000 b tn
00000000000000030
C W N 0' G
m cno 9
00000000000000000 b
.................
00000000000000000
ccccccccCc~cc
ccccc Bq
00000000000000000
-I
zu
CO
x
z z z
C C
x
c m m
C
X
~ m m
I) a
x
.................
OGOOOOOOOOOOOGOOO
00000000000000000 X X
00000000000000000 C
xD-
CZ
Pm
...*..
CCCCLC
4
rl
OCOOOO 0
n
0
a
'I-
1 Q
a 00000000000000000000 4
Z NNNNNNNNNNNNNNNNNNNN 3
WNNYNNWNNNNNNNNNNNNW n
>>ZZ>>><Jc<Z>'-'->~~dT4 Q.
zza~ZZzOOOW~>>~3L9u3 0
uus z u m a a a m J q a J t + I - r n c J 0
3
1 0
I : .
0
eus DATA
GENERATION LOAO SHLNT
US N A M ~ U~JLTS WGLE MY HVAR rn MVAR MVPR
.....................................
3 MAN220
3
5
MAN220
TIC220
I
49'339
* r:g*%:$
39:647
2 Xb:%
MISCATCh
49'339
51:092
-0.014
3 9 47
-49:997
-0.004
3 MAN220
3 MAN220
1 INV220
I 60.722
I-49.244
60.722
-49.933
-49.933
-42.559
1 INV220 -49 244 -42.659
MISCATCh -0:067 -0.016
i! ROX220
MISVATCh
.....................................
7 BEN220 0.006 223.602
MISVATCh -0.006 0.000
.........................................................................
12 LIV220 2 1 370 92.023
7 BEhZ20 88'957 0.732
7 BEN220 88:957 0.732
10 A V I O l l -199.261 -93.686
MISV ATCh -0.023 -0.000
.........................................................................
4 AVI220 199.993 115.661
MISV A T C r 0.007 0.000
17 ~ ~ 1 2 2 0 350.004 113.378
MISCATCh -0.004 0.000
--I---------------------------------
13 W228
15 T E K O l l
ii3:118 -i8:0#
-149.320 12. 90
MISVATCh -0.007 0.001
--i
ISL 220
OHA TWI 2 2 0
I
LIV 2 2 0
MAN 220
TIW 2 2 0
Fig. 4.9. Reduced primary a.c. system for the South Island of
New Zeaiand
4.13. REFERENCES
1. B. Stott, 1974. 'Review of load-flow calculation methods', Proc. I EEE, 62 (7),9 16-929.
2. J. B. Ward and H. W. Hale, 1956. 'Digital computer solution of power-flow problems',
Trans. A1EE., PAS-75, 398-404.
3. H. E. Brown, G. K. Carter, H. H. Happ and C. E. Person, 1963. 'Power-flow solution
by impedance matrix iterative method', IEEE Trans. PAS-82, 1-10.
4. J. E. Van Ness and J. H. Griflin, 1961. 'Elimination methods for load-flow studies',
Pans. A1EE., PAS-80, 299-304.
5. W. F. Tinney, C. E. Hart, 'Power flow solution by Newton's method', IEEE Trans.,
PAS86 (1 l), 1449-1460.
6. E. C. Ogbuobiri, W. F. Tinney and J. W. Walker, 1970. 'Sparsity-directed decom-
position for Gaussian elimination on matrices', IEEE, Trans., PAS-89 (I), 141-150.
7. B. Stott and E. Hobson, 1971. 'Solution of large power-system networks by ordered
elimination: a comparison of ordering schemes', Proc. IEE, 118 (I), 125-134.
8. W. F. Tinney and J. W. Walker, 1967. 'Direct solutions of sparse network equations
by optimally ordered triangular factorization', Proc. IE EE, 55 ( 1I), 1801-1 809.
9. S. T. Despotovic, 1974.'A new decoupled load-flow method', IEEE Trans., PAS93 (3),
884-89 1.
10. B. Stott, 1972. 'Decoupled Newton load flow', IEEE Trans., PAS-91, 1955-1959.
11. B. Stott and 0.Alsac, 1974. 'Fast decoupled load flow', IEEE Trans., PAS93 (3), 859-
869.
12. J. Medanic and B. Avramovic, 1975.'Solution of load-flow problems in power systems
by &-couplingmethod', Proc IEE, 122 (8), 801-805.
13. B. Stott, 1972. 'Power-system load flow' (M.Sc. lecture notes). University of
Manchester Institute of Science and Technology.
Three-Phase Load Flow
5.1. INTRODUCTION
For most purposes in the steady-state analysis of power systems, the system
unbalance can be ignored and the single-phase analysis described in Chapter 4 is
adequate. However, in practice it is uneconomical to balance the load completely
or to achieve perfectly balanced transmission system impedances, as a result of
untransposed high-voltage lines and lines sharing the same right of way for
considerable lengths.
Among the effects of power-system unbalance are: negative sequence currents
causing machine rotor overheating, zero sequence currents causing relay
maloperations and increased losses due to parallel untransposed lines.
The use of long-distance transmission motivated the development of analytical
techniques for the assessment of power-system unbalance. Early techniques(').(2)
were restricted to the case of isolated unbalanced lines operating from known
terminal conditions. However, a realistic assessment of the unbalanced operation
of an interconnected system, including the influence of any significant load
unbalance, requires the use of three-phase load-flow a l g o r i t h m ~ . ( ~ ) . (The
~)?(~)
object of the three-phase load flow is to find the state of the three-phase power
system under the specified conditions of load, generation and system
configuration.
The basic three-phase models of system plant and the rules for their
combination into overall network admittance matrices, discussed in Chapter 2,
are used as the framework for the three-phase load flow described in this chapter.
The storage and computational requirements of a three-phase load-flow
program are much greater than those of the corresponding single-phase case. The
need for efficient algorithms is therefore significant even though, in contrast to
single-phase analysis, the three-phase load flow is likely to remain a planning,
rather than an operational exercise.
The basic characteristics of the Fast Decoupled Newton-Raphson algorithms
described in Chapter 4, have recently been shown(@to apply equally to the three-
phase load-flow problem. Consequently, this algorithm is now used as a basis
for the development of an efficient three-phase load-flow program. When the
program is used for post-operational studies of important unbalanced situations
. ... , I
1
5.2. NOTATION
A clear and unambiguous identification of the three-phase vector and matrix
elements requires a suitable symbolic notation using superscripts and subscripts.
The a.c. system is considered to have a total of n busbars where:
+
n = nb ng
n b is the number of actual system busbars
ng is the number of synchronous machines
Subscripts i, j, etc. refer to system busbars as shown in the following examples:
i = 1, nb identifies all actual system busbars, i.e. all load busbars plus all
generator terminal busbars.
i = n b + 1, nb + ng - 1 identifies all generator internal busbars with the
exception of the slack machlne.
i = n b + ng identifies the internal busbar at the slack machine.
The following subscripts are also used for clarity:
reg-refers to a voltage regulator
int-refers to an internal busbar at a generator
gen-refers to a generator
Superscripts p, m identify the three phases at a particular busbar.
where
and 'a' is the complex operator ej2"I3.The phase component impedance matrix is
thus
The phase component model of the generator is illustrated in Fig. 5.l(a) The
machine excitation acts symmetrically on the three phases and the voltages at the
internal or excitation busbar form a balanced three-phase set, i.e.
and
For three-phase load flow the voltage regulator must be accurately modelled as
it influences the machine operation under unbalanced conditions. The voltage
regulator monitors the terminal voltages of the machine and controls the
excitation voltage according to some predetermined function of the terminal
voltages. Often the positive sequence is extracted from the three-phase voltage
measurement using a sequence filter.
Before proceeding further it is instructive to consider the generator modelling
from a symmetrical component frame of reference. The sequence network model
of the generator is illustrated in Fig. 5.l(b). As the machine excitation acts
symmetrically on the three-phases, positive sequence voltages only are present at
the internal busbar.
The influence of the generator upon the unbalanced system is known if the
voltages at the terminal busbar are known. In terms of sequence voltages, the
positive sequence voltage may be obtained from the excitation and the positive
sequence voltage drop caused by the flow of positive sequence currents through
the positive sequence reactance. The negative and zero sequence voltages are
derived from the flow of their respective currents through their respective
impedances. It is important to note that the negative and zero sequence voltages
are not influenced by the excitation or positive sequence impedance.
There are infinite combinations of machine excitation and machine positive
sequence reactance which will satisfy the conditions at the machine terminals and
give the correct positive sequence voltage. Whenever the machine excitation must
be known (as in fault studies) the actual positive sequence impedance must be
used. For load flow however, the excitation is not of any particular interest and
the positive sequence impedance may be arbitrarily assigned to any value.@)The
positive sequence impedance is usually set to zero for these studies.
I V term,
I f v e sequence
-ve sequence
zero sequence
Thus the practice with regard to three-phase load flow in phase coordinates, is
to set the positive sequence reactance to a small value in order to reduce the
excitation voltage to the same order as the usual system voltages with a
corresponding reduction in the angle between the internal busbar and the
terminal busbar. Both these features are important when a fast decoupled
algorithm is used.
Therefore, in forming the phase component generator. model using equation
(5.3.4), an arbitrary value may be used for Z , but the actual values are used for 2,
and 2,. There is no loss of relevant information as the influence of the generator
upon the unbalanced system is accurately modelled.
The nodal admittance matrix, relating the injected currents at the generator
busbars to their nodal voltages, is given by the inverse of the series impedance
matrix derived from equation (5.3.4).
Specified variables
The following variables form a minimum and sufficient set to define the three-
phase system under steady-state operation:
- The slack generator internal busbar voltage magnitude Vin,j where j =
+
nb ng. (The angle Ointj is taken as a reference.)
- The internal busbar voltage magnitude Vintj and angles Ointj at all other
+
generators, i.e. j = n b 1, n b + ng - 1 .
- The three voltage magnitudes (Vf') and angles (Of)at every generator terminal
busbar and every load busbar in the system, i.e. i = 1, n b and p = 1,3.
Only two variables are associated with each generator internal busbar as the
three-phase voltages are balanced and there is no need for retaining the
redundant voltages and angles as variables. However, these variables are retained
to facilitate the calculation of the real and reactive power mismatches. The
equations necessary to solve for the above set of variables are derived from the
specified operating conditions, i.e.
- The individual phase real and reactive power loading at every system busbar.
- The voltage regulator specification for every synchronous machine.
- The total real power generation of each synchronous machine, with the
Derivation of equations
The three-phase system behaviour is described by the equation
= (Qf)""- Vf
n
k= 1 m= 1
x
3
V r [Gim sin Of7- B f j cos 0im] (5.3.9)
= - z xz
p=l
3
Vint
n 3
k=l m=l
V r [Gjql"cos 04r + Bjql"sin 0 4 3
(5.3.11)
where, although the summation for k is over all system busbars, the
mutual terms Gj, and Bjkare nonzero only when k is the terminal busbar of
the jth generator.
It should be noted that the real power specified for the generator is the total real
power at the internal or excitation busbar whereas in actual practice the specified
quantity is the power leaving the terminal busbar. This in effect means that the
generator's real power loss is ignored.
The generator losses have no significant influence on the system operation and
may be calculated from the sequence impedances at the end of the load-flow
solution, when all generator sequence currents have been found. Any other
method would require the real power mismatch to be written at busbars remote
from the variable in question, that is, the angle at the internal busbar. In addition,
inspection of equations (5.3.8) and (5.3.1 1) will show that the equations are
identical except for the summation over the three phases at the generator internal
busbar.
That is, the sum of the powers leaving the generator may be calculated in
exactly the same way and by the same subroutines as the power mismatches at
other system busbars. This is possible because the generator internal busbar is
not connected to any other element in the system. Inspection of the Jacobian
submatrices derived later will show that this feature is retained throughout the
study. In terms of programming the generators present no additional complexity.
Equations (5.3.8) to (5.3.11) form the mathematical formulation of the three-
phase load flow as a set of independent algebraic equations in terms of the system
variables.
The solution to the load-flow problem is the set of 'variables which, upon
substitution, make the left-hand side mismatches in equations (5.3.8) to (5.3.1 1)
equal to zero.
5.4. FAST DECOUPLED THREE-PHASE ALGORITHM
The standard Newton-Raphson algorithm may be used to solve equations (5.3.8)
to (5.3.11). This involves an iterative solution of the matrix equation
for the right-hand side vector of variable updates. The right-hand side matrix in
equation (5.4.1) is the Jacobian matrix of first order partial derivatives.
Following decoupled single-phase load-flow practice, the effects of A 8 on
reactive power flows and A V on real power flows are ignored. Equation (5.4.1)
may therefore be simplified by assigning
and
[C]=[GI = O
In addition, the voltage regulator specification is assumed to be in terms of the
terminal voltage magnitudes only and therefore,
[Dl=[HI = 0
Equation (5.4.1) may then be written in decoupled form as:
where [Fill = 0 for all j # 1because the jth generator has no connection with the
lth generator's internal busbar, and
p= 1
3 3
+ m2= 2 (Vin,J2[Gf;'"sin O
1 p=l
r - BEmcos Of;"]
m f P
where
K$" = V r Vf [Gfp sin Of; - Bfp cos O$"]
except
K g " = - B;:'(v;)2 + Q;:
let [L;] = V;[L;]' where k is the terminal busbar of the jth generator and
L$ = 0 otherwise.
Jacobian approximations
Approximations similar to those applied to the single-phase load flow are
applicable to the Jacobian elements as follows:
(i) at all nodes (i.e. all phases of all busbars)
(iii) Moreover the phase-angle unbalance at any busbar will be small and
hence an additional approximation applies to the three-phase system, i.e.
z f 120" for p # m
(iv) Finally, as a result of (ii) and (iii) the angle between different phases of
connected busbars will be approximately 120", i.e.
or
cos @" z - 0.5
and
sin Of? z 0.866
These values are modified for the f 30" phase-shift inherent in the
star-delta connection of three-phase transformers.
The final approximation (iv), necessary if the Jacobians are to be kept constant,
is the least valid, as the cosine and sine values change rapidly with small angle
variations around 120 degrees. This accounts for the slower convergence of the
phase unbalance at busbars as compared with that of the voltage magnitudes and
angles.
It should be emphasised that these approximations apply to the Jacobian
elements only, i.e. they do not prejudice the accuracy of the solution nor do they
restrict the type of problem which may be attempted.
Applying approximations (i) to (iv) to the Jacobians and substituting into
equations (5.4.2) and (5.4.3) yields,
and
where
with
Recalling that [LEI' = [aA Vreg,/a V,"], as Vreg is normally a simple linear
function of the terminal voltages, [L'] will be a constant matrix.
Therefore, the Jacobian matrices [B'] and [B"] in equations (5.4.6) and (5.4.7)
have been approximated to constants.
Zero diagonal elements in equation (5.4.7) may result from the ordering of the
equations and variables. This feature causes no problems if these diagonals are
not used as pivots until the rest of the matrix has been factorized (by which time,
fill-in terms will have appeared on the diagonal). This causes a minor loss of
efficiency as it inhibits optimal ordering for the complete matrix. Although this
could be avoided by reordering the equations, the extra program complexity is
not justified.
Based on the reasoning of Stott and Alsac,") which proved successful in the
single-phase load flow the [B'] matrix in equation (5.4.6) is further modified by
omitting the representation of those elements that predominantly affect MVAR
flows.
The capacitance matrix and its physical significance is illustrated in Fig. 5.2, for
a single three-phase line. With n capacitively-coupled parallel lines the matrix will
be 3n x 3n.
In single-phase load flows the shunt capacitance is the positive sequence
capacitance which is determined from both the phase-to-phase and the phase-to-
earth capacitances of the line. It therefore appears that the entire shunt
capacitance matrix predominantly affects MVAR flows only. Thus, following
single-phase fast decoupling practice the representation of the entire shunt
capacitance matrix is omitted in the formulation of [B']. This increases
dramatically the rate of real power convergence.
With capacitively coupled three-phase lines the interline capacitance in-
fluences the positive sequence shunt capacitance. However, as the values of
interline capacitances are small in comparison with the self capacitance of the
phases, their inclhsion makes no noticeable difference. The effective tap of fi
introduced by the star-delta transformer connection is interpreted as a nominal
tap and is therefore included when forming the [B'] matrix.
( ii )
==F=T
Calculate [AO/V],[AV,, ] Output
-
Fig. 5.3. Iteration sequence for three-phase a.c. load flow
A further difficulty arises from the modelling of the star-gldelta transformer
connection. The equivalent circuit, illustrated in Section 2.6 shows that large
shunt admittances are effectively introduced into the system. When these are
excluded from [B'], as for a normal shunt element, divergence results. The entire
transformer model, must therefore be included in both [B'] and [B"].
With the modifications described above the two final algorithmic equations
may be concisely written, i.e.
Data input
The input data routine implements the system modelling techniques described in
Chapter 2 and forms the system admittance model from the raw data for each
system component. Examples of the raw data are given in Section 5.7. with
reference to a particular test system.
The structure and content of the constant Jacobians B' and B" are based upon
the system admittance matrix and are thus formed simultaneously with this
matrix.
Both the system admittance matrix and the Jacobian matrices are stored and
processed using sparsity techniques which are structured in 3 x 3 matrix blocks
to take full advantage of the inherent block structure of the three-phase system
matrices.
Data Input
Form and store system
admittance matrix, [ B ' ]
and [ B " ]Jacobian
I matrices. (1645)
I
1
I Z Z i z e [ B ' ]and [ B * ]
I
I Initialize three-phase
system voltages (200) I
procedure (see Fig. 5.3)
(450)
Starting values
Starting values are assigned as follows:
- The nonvoltage controlled busbars are assigned 1 p.u. on all phases.
- At generator terminal busbars all voltages are assigned values according to
the voltage regulator specifications.
- All system busbar angles are assigned 0, - 120°, + 120" for the three phases
respectively.
- The generator internal voltages and angles are calculated from the specified
real power and, in the absence of better estimates, by assuming zero reactive
power. For the slack machine the real power is estimated as the difference
between total load and total generation plus a small percentage(say 8 percent)
of the total load to allow for losses.
For cases where convergence is excessively slow or difficult it is recommended
to use the results of a single-phase load flow to establish the starting values. The
values will, under normal steady-state unbalance, provide excellent estimates for
all voltages and angles including generator internal conditions which are
calculated from the single-phase real and reactive power conditions.
Moreover, as a three-phase iteration is more costly than a single-phase
iteration, this practice can be generally recommended to provide more efficient
overall convergence and to enable the more obvious data errors to be detected at
an early stage.
For the purpose of investigating the load-flow performance, flat voltage and
angle values are used in the examples that follow.
Iterative solution
The iterative solution process (Fig. 5.3) yields the values of the system voltages
which satisfy the specified system conditions of load, generation and system
configuration.
Output results
The three-phase busbar voltages, the line power flows and the total system losses
are calculated and printed out. An example is given in Table 5.6. In addition the
sequence components of busbar voltages are also calculated as these provide a
more direct measure of the unbalance present in the system under study.
ltemtion lterition
Fig. 5.5. Power convergence patterns for three-phase and single-phase load flow
VoltaJe
(plr.1 A
1.04 Three-phase busbar
voltages
3
2
2
1 1'
I I t
2 3 4 5 Iteration
(i 1
Voltage
( pu.)
Generator
Synchronous
condenser
Genemtor
excitation voltopes
-r
&
MAN GN
MAN M4
Excitation
voltages
Secondory ( S 1
Primary ( P
MAN 220
0
Fig. 5.9. Test system 3 x 3 matrix representation
5.7. TEST SYSTEM AND RESULTS
A single-line diagram of the test system under consideration is illustrated in
Fig. 5.7. Some features of interest are:
- An example of a line sectionalization is included. One section contains four
mutually coupled three-phase power lines. The other section contains two
sets of two mutually coupled three-phase lines.
- All parallel lines are represented in their unbalanced mutually coupled state.
- Both transformers are star-delta connected with the star neutrals solidly
earthed. Tap ratios are present on both primary and secondary sides.
The system is redrawn in Fig. 5.8 using 3 x 3 compound coil notation and
substituting for the generator and line models. Following this, Fig. 5.9 illustrates
the system graphically in terms of 3 x 3, 6 x 6 and 12 x 12 matrix blocks,
The effect of subsystem 3 (the synchronous condenser) is not included in the numerical example
Busbar names
Leakage Tap ratio
primary secondary reactance prim.
Both these matrices are in p.u. for the total line length.
Subsystem 7
This subsystem consists of a pair of parallel, mutually coupled three-phase lines.
These lines are represented in their-natural coupled unbalanced state.
Terminal busbars
Line 1 INV220 - TIW220
Line 2 INV220 - TIW220
The series impedance matrix for the total length (2,) is:
Line 1 Line 2
b b
0.0012
Line 1 b
+jO.008
Line 2 b
+jO.0061
The shunt admittance matrix for the total line length is:
Line 1 Line 2
a b c a b c
Line 1 b
Line 2 b
The lower diagonal half only is shown as all line matrices are symmetrical.
Subsystem 8
Sectionalized mutually coupled lines: Section 1 consists of four mutually coupled
three-phase lines and has 12 x 12 characteristic matrices, [Z,,] and [ Y , , ] , as
indicated in the system diagrams. These are given in Figs. 5.1 1 and 5.12 in per unit
length of line and Section 1 is taken as having a length of 0.75 units.
Section 2 consists of two sets of two mutually coupled three-phase lines. To
ensure consistent dimensionality with Section 1, the second section is considered
as being composed of four mutually coupled three-phase lines, the elements
representing the coupling between the two separate double-circuit lines being set
to zero. The characteristic matrices for Section 2 become,
Section 2 4
1 2 3 4
c Y*I
12 x 12
Section 2
where 101 is a matrix of zeros. The submatrices are labelled as those in Fig. 5.8.
These 12 x 12 matrices are given in Fig. 5.13 and Fig. 5.14. in per unit length of
line. Section 2 is taken as having a length of 0.25 units.
Once the overall admittance matrix for the combined sections has been found
it must be stored in full. This is to enable calculation of the power flows in the four
individual lines. The matrix is modified, as described in Section 2.5. for the
terminal connections. This modified matrix is the subsystem admittance matrix
to be combined into the overall system admittance matrix.
Line 1 Line 4
Line I b
Line 2 b
Line 3 b
Line 4 b
Case 10.0 1O
. 0.1
i
11
...
111
iv
v
vi
vii
viii
ix
X
Table 5.3. Sequence components of busbar voltages
Case (vii)
Case (viii)
- An apparent gain in active power flow in any one phase. This power flows
through the mutual coupling terms between phases. The overall active power
shows a net loss as expected for a realistic system.
INV220
ROX220
MAN220
MAN014
TIW220
ROXO 11
MAN.GN
ROX. GN
Table 5.6 Computed print-out power flows
where
P is a vector of the voltage magnitudes at all a.c. system busbars.
0 is a vector of the angles at all a.c. system busbars (except the reference bus
which is assigned 8 = 0).
2 is a vector of d.c. variables.
Chapter 4 has described the use of P and 8 as a.c. system variables and the
selection of d.c. variables 2 is discussed in Section 6.3.
The development of a Newton-Raphson based algorithm requires the
formulation of n independent equations in terms of the n variables.
The equations which relate to the a.c. system variables are derived from the
specified a.c. system operating conditions. The only modification required to the
usual real and reactive power mismatches occurs for those equations which relate
to the converter terminal busbars. These equations become:
where
Pterm(ac)is the injected power at the terminal busbar as a function of the a.c.
system variables
Pterm(dc)is the injected power at the terminal busbar as a function of the d.c.
system variables
c g m is the usual a.c. system load at the busbar.
and similarly for Qterm(dc)and Qterm(ac)
The injected powers Qterm(dc)and Pterm(dc)are functions of the convertor a.c.
terminal busbar voltage and of the d.c. system variables, i.e.
The equations derived from the specified a.c. system conditions may therefore be
summarized as :
where the mismatches at the converter terminal busbars are indicated separately.
A further set of independent equations are derived from the d.c. system
conditions. These are designated,
where the subscript 'term' refers to the converter a.c. terminal busbar.
Converter variables
Under balanced conditions similar converter bridges attached to the same a.c.
terminal busbar will operate identically regardless of the transformer connection.
They may therefore be replaced by an equivalent single bridge for the purpose of
single-phase load flow analysis. With reference to Fig. 6.1 the set of variables
illustrated, representing fundamental frequency or d.c. quantities permits a full
description of the converter system operation.
An equivalent circuit for the converter is shown in Fig. 6.2 which includes the
modification explained in Section 6.2 as regards the position of angle reference.
The variables, defined with reference to Fig. 6.2 are:
Vterm& converter terminal busbar nodal voltage (phase angle referred to
converter reference)
-
Is le"
I-
I
IP
--t
- -
I
Fig. 6.1. Basic d.c. convertor (angles refer to a.c. system reference)
Fig. 6.2. Single-phase equivalent circuit for basic converter. (Angles
referred to d.c. reference)
I,(p.u.) f i J~.I,(~.U.)
= -.
n
and if commutation overlap is taken into account, as described in Chapter 3, this
equation becomes
where k is very close to unity. In load-flow studies, equation (6.3.2) can be made
sufficiently accurate in most cases by letting:
k = 0.995
Derivation of equations
The following relationships are derived for the variables defined in Fig. 6.2. The
equations are in per unit.
(i) The fundamental current magnitude on the converter side is related to the
direct current by the equation
(iii) The d.c. voltage may be expressed in terms of the a.c. source commutating
voltage referred to the transformer secondary, i.e.
(v) The assumptions listed at the beginning of this section prevent any real
power of harmonic frequencies at the primary and secondary busbars.
Therefore, the real power equation relates the d.c. power to the
transformer secondary power in terms of fundamental components only,
i.e.
VdeId= E.I;cos $ (6.3.7)
(vi) As the transformer is lossless, the primary real power may also be equated
to the d.c. power, i.e.
These control equations are simple and are easily incorporated into the solution
algorithm. In addition to the usual control modes, nonstandard modes such as
specified a.c. terminal voltage may also be included as converter control
equations (see Section 6.5).
During the iterative solution procedure the uncontrolled converter variables
may go outside prespecified limits. When this occurs, the offending variable is
usually held to its limit value and an appropriate control variable is freedj4)
Inverter operation
All the equations presented so far are equally applicable to inverter operation.
However, during inversion it is the extinction advance angle (y) which is the
subject of control action and not the firing angle (a). For convenience therefore,
equation R(2) of (6.3.11) may be rewritten as
Unified solution
The unified method gives recognition to the interdependence of a.c. and d.c.
system equations and simultaneously solves the complete system. Referring to
equation (6.2.7) the standard Newton-Raphson algorithm involves repeat
solutions of the matrix equation:
and
Applying the a.c. fast decoupled assumptions to all Jacobian elements related
to the a.c. system equations, yields:
where all matrix elements are zero unless otherwise indicated. The matrices [B']
and [B"] are the usual single-phase fast decoupled Jacobians and are constant in
value. The other matrices indicated vary at each iteration in the solution process.
A modification is required for the element indicated as B:i in equation (6.4.6).
This element is a function of the system variables and therefore varies at each
iteration.
The use of an independent angle reference for the d.c. equations results in
i.e. the diagonal Jacobian element for the real power mismatch at the converter
terminal busbar depends on the a.c. equations only and is therefore the usual fast
decoupled B' element.
In addition,
aR/a&erm =0
therefore: DD = 0
Similarly
1
CAA'I = - - [ a ~ ~ ~ ~ ~ ~ l a x l
Vterm
1
[AA"] = -[aAQterm/a2]
Vterm
1
B!! = -CaAQtermla Vterml
" Vterm
In the above formulation the d.c. variables 2 are coupled to both the real and
reactive power a.c. mismatches. However, equation (6.4.6) may be separated to
enable a block successive iteration scheme to be used.
134
The d.c. mismatches and variables can be appended to the two fast decoupled
a.c. equations in which case the following two equations result
Bii AA"
.
I ~ a l c u l a t ePkm (dc) and d r residuals ( 5 1 I
I
1
*
I Calculate A d ( a . c . svstem onlv
+
I Calculate qerm
(d.c.1 and d c residuals ( #1
I [sropb
Converge
*
I Form d.c. Jacobian Matrix
w
Fig. 6.3. Flow chart for unified single-phase a.c.-d.c. load flow
The Jacobian elements related to the d.c. variables are non-constant and must
be re-evaluated at each iteration. It is therefore necessary to separate the constant
and nonconstant parts of the equations for the solution routine.
Initially, the a.c. fast decoupled equations are formed with the d.c. link ignored
(except for the minor addition of the filter reactance at the appropriate a.c.
busbar). The reactive power mismatch equation for the a.c. system is:
136
where
AQIerm = Q:Lm - Qterm(ac).
is the mismatch calculated in the absence of
the d.c. converter.
and
B" is the usual constant a.c. fast decoupled Jacobian.
After triangulation down to, but excluding the busbars to which d.c. converters
are attached, equation (6.4.11) becomes
where (AQ/ V)" and (AQterm/ Vterm)"signify that the left-hand side vector has been
processed and matrix B"' is the new matrix B" after triangulation.
This triangulation (performed before the iterative process) may be achieved
simply by inhibiting the terminal busbars being used as pivots during the optimal
ordering process.
The processing of AQ indicated in the equation is actually performed by the
standard forward reduction process used at each iteration.
The d.c. converter equations may then be combined with equation (6.4.12)as
follows :
AV
A2
where
Sequential method
The sequential method results from a further simplification ofthe unified method,
i.e. the a.c. system equations are solved with the d.c. system modelled simply as a
real and reactive power injection at the appropriate terminal busbar. For a d.c.
is written as one of the two control equations. This would lead to a zero row in
equation (6.4.7) and therefore during the solution of equation (6.4.7) some other
variable (e.g. tap ratio) must be specified instead. Although d.c. convergence is
marginally slower for the PDC, QDC iteration, the d.c. system is overconverged in
this iteration scheme and the overall convergence rate is practically unaffected.
With the sequential method equation (6.5.1)cannot be written. The terminal
busbar is specified as a P - V busbar and the control equation
is used, where QsLrn(dc) is taken as the reactive power required to maintain the
voltage constant. The specified reactive power thus varies at each iteration and
this discontinuity slows the overall convergence.
%+
82-Bus
t
North Island
system
This example indicates the ease of extension to the multiple converter case.
Sequential methods
Specified d.c. Unified methods
link constraints (5 variables) (5 variables) (4 variables)
m-rectifier end
n-inverter end ~ P D C QDC
. 2 ~QDC
. ~P,Q.DC~P.DC,Q.DC ~P.Q.DC~P.DC,Q.DC
Converter I Converter 2
Various control strategies have been applied to the link and the convergence
results for the various algorithms are given in Table 6.1. The number of iterations
(i,j) should be interpreted as follows:
- i is the number of reactive power-voltage updates required.
-j is the number of real power-angle updates.
Although the number of d.c. iterations varies for the different sequences, this is of
secondary importance and may if required be assessed in each case from the
number 0fa.c. iterations. In this respect, a unified QDC iteration is equivalent to a
Q iteration and a DC iteration executed separately.
The d.c. link data and specified controls for Case 1 are given in Table 6.2 and
the corresponding d.c. link operation is illustrated in Fig. 6.6. The specified
conditions for all cases are derived from the results of Case 1. Under those
conditions, the a.c. system in isolation, (with each converter terminal modelled as
an equivalent a.c. load) requires (4, 3) iterations. The d.c. system in isolation
(operating from fixed terminal voltages) requires two iterations under all control
strategies.
Bus 5 Bus 4
All angles are in degrees. D.C.voltages and current are in kVand Amp respectively.
D.C.resistance is in ohms. AC.powers (P,Q)are in MW and MVARs.
Fig. 6.6. D.c. link operation for Case I
Unified cases
The results in Table 6.1 show that the unified methods provide fast and reliable
convergence in all cases.
For the unified methods 1 and 2 the number of iterations did not exceed the
number required for the a.c. system alone.
Sequential cases
The sequential method (P,Q,DC) produces fast and reliable convergence
although the reactive power convergence is slower than for the a.c. system alone.
With the removal of the variable 4, Q,,,, (dc) converges faster but the
convergence pattern is more oscillatory and an overall deterioration of a.c.
voltage convergence results.
With the second sequential method, (P, DC, Q, DC) convergence is good in all
cases except 2 and 3, i.e. the cases where the transformer tap and d.c. voltage are
specified at the inverter end. However, this set of specifications is not likely to
occur in practice.
Q Filters
Sequential Sequential
P.Q,DC P,Q,DC
m-rectifier Unified Unified
n-inverter P,QDC (i) (ii) P,QDC (1) (ii)
(i) using the five variable formulation; (ii) using the four variable formulation.
The reactive power compensation of the filters was adjusted to give similar d.c.
operating conditions as previously.
The number ofiterations to convergence for the most promising algorithms are
shown in Table 6.3 for the control specifications corresponding to Cases 1 to 4 in
the previous results.
The different nature of the sequential and unified algorithms is clearly
demonstrated. The effect of the type of converter control is also shown. For Case
11, both the d.c. real power and the d.c. reactive powers are well constrained by
the converter control strategy. Convergence is rapid and reliable for all methods.
In all other cases, where the control angle at one or both converters is free, an
oscillatory relationship between converter a.c. terminal voltage and the reactive
power of the converter is possible. This leads to poor performance of the
sequential algorithms.
To illustrate the nature of the iteration, the convergence pattern of the
converter reactive power demand and the a.c. system terminal voltage of the
rectifier is plotted in Fig. 6.8.
A measure of the strength of a system in a load-flow sense is the short circuit to
converter power ratio (SCR) calculated with all machine reactances set to zero.
This short-circuit ratio is invariably much higher than the usual value.
In practice, converter operation has been considered down to a SCR of 3. A
survey ofexisting schemes shows that, almost invariably, with systems of very low
SCR, some form of voltage control, often synchronous condensers, is an integral
part of the converter installation. These schemes are therefore often very strong in
a load-flow sense.
It may therefore be concluded that the sequential integration should converge
in all practical situations although the convergence may become slow if the
system is weak in a load-flow sense.
Fig. 6.8. Convergence pattern for a.c.-d.c. load flow with weak a.c.
system. (i) Sequential method (P,Q, DC five variable); (ii) Unified
method (P, QDC)
Discussion of convergence properties
The overall convergence rate of the a.c.-d.c. algorithms depends on the successful
interaction of the two distinct parts. The a.c. system equations are solved using
the well-behaved constant tangent fast decoupled algorithm, whereas the d.c.
system equations are solved using the more powerful, but somewhat more erratic,
full Newton-Raphson approach.
The powerful convergence of the Newton-Raphson process for the d.c.
equations can cause overall convergence difficulties. If the first d.c. iteration
occurs before the reactive power-voltage update then the d.c. variables are
converged to be compatible with the incorrect terminal voltage. This introduces
an unnecessary discontinuity which may lead to convergence difficulties in the
sequential method. In the unified approach the powerful convergence of the d.c.
equations is dampened by the reflection of the a.c. mismatches onto the changes
in d.c. variables. This gives faster and better behaved convergence. The solution
time of the d.c. equations is normally small compared to the solution time of the
a.c. equations. The relative efficiencies of the alternative algorithms may there-
fore be assessed by comparing the total numbers of voltage and angle updates.
In general, those schemes which acknowledge the fact that the d.c. variables are
strongly related to the terminal voltage give the fastest and most reliable
performance. The unified methods are the more reliable and the P, QDC solution
is the most efficient. Of the sequential methods, the (P,Q, DC) solution is only
marginally inferior to the unified method.
When the busbar to which the converters are attached is voltage controlled (as
is often the case) the two approaches become virtually identical as the interaction
between a.c. and d.c. systems is much smaller.
The only computational difference between a unified and a comparable
sequential iteration is that the d.c. Jacobian equation for the unified method
(equation 6.4.14) is slightly larger. The difference is one additional row and
column for each converter present. In terms of computational cost per iteration
the corresponding unified and sequential algorithms are virtually identical.
(i) In cases where the a.c. system is strong, both the unified and sequential
algorithms may be programmed to give fast and reliable convergence.
(ii) If the a.c. system is weak the sequential algorithm is susceptible to
convergence problems. Thus, in general, the unified. method is recom-
mended due to its greater reliability.
-
CI
1 0 0 0 0
0 . 0 0 0 0
C l o o o o
u o o o o
-
C(
1 0 0 0 0
m o o 0 0
C l o o o o
m
.....
u o o
o o
o
0
o
0
-
rn
C(
1 0 0 0 0
P \ o o o o
a 0 0 0 0
.....
u o o o o
-.
C(
1 0 0 0 0
9 0 0 0 0
G O O O O
.....
u o o o o
C( I
1 0 0 0 0
Y l o o o o
a 0 0 0 0
~ 0 0 0 0
-
C(
1 0 0 0 0
- 0 0 0 0
G O O O O
.....
0
& O O O O
000000000m000000000
L 000000000ln000000000 - 0 0 0 0
a
i- ...................
ooooooooooooooooooo
"dddd-ld-l"""-ld""""d"
G
..
1 0 0 .
"
0 0
m o - l o o
C l o o o o
u
.....o o
P Y 0 0 6 0
o o
n
C(
1 0 0 0 0
r.4
0 0 0 0
U O O O O
.....
C L O O O
r 4 0 0 0 0
O
a
SOLUTION C O N V E R G E D IN 7 P - D AND 6 n-v ITERATIONS
CONVERTER AC VOLTAGE TRANSFORMER TAP CONTROL ANGLE COHMUTATION ANGLE DC CURRENT DC VOLTAGE
(K-VOLTS) (PER CENT) (DEGS) (DEGS) (K-AMPS) (K-VOLTS)
97.08 -5.46 8.00 19.93 0.81 -245.67
POWER TRANSFERS
L I N K TERMINAL POWER = -200.00 HW 110.89 HVAR
FROM TRANSFORMER TO CONVERTER = - 2 0 0 . 0 0 HW 74.63 HVAR
REACTIVE POWER OF F I L T E R S = 120.83 HVAR
COWIFTATION R E A C T L N E (P.U. )
. F I R I N G ANGLE:MINIWUM(D€G)
MAXICUM(DEG)
TRANSFORMER TAP:MINIMUM(P.C.)
MAXIMUM ( P o c o )
INCRIiMENT
F l L T l i R REACTANCE(P.U.)
NUMBER OF BRIDGES I N SERIES*
15 1
6.10 REFERENCES
1. B. Scott and 0.Alsac, 1974.'Fast decoupled load flow', Trans. IEEE, PAS-93 (3), 859-
869.
2. H. Sato and J. Arrillaga, 1969. 'Improved load-flow techniques for integrated a.c.-d.c.
systems', Proc. IEE, 116 (4), 525-532.
3. J. Reeve, G. Fahmy, and B. Stott, 1976. 'Versatile load-flow method for multi-terminal
h.v.d.c. systems', Paper F76-354-1. IEEE PES Summer Meeting, Portland,
4. B. Stott, 'Load flow for a.c. and integrated a.c.-d.c. systems.' Ph.D. Thesis, University of
Manchester, 1971.
5. D. A. Braunayal, L. A. h a f t , and J. L. Whysong, 1976. 'Inclusion of the converter and
transmission equations directly in a Newton power flow', IEEE Trans. PAS-75 (I), 76-
88.
6. J. Arrillaga, B. J. Harker, K. S. and Turner, 1980. 'Clarifying an ambiguity in recent
a.c.-d.c. load-flow formulations', Proc. IEE, 127, Pt. C (5), 324-325.
7. P. Bodger, 1977. 'Fast decoupled a.c. and a.c.-d.c. load flows.' Ph.D. Thesis, University
of Canterbury, New Zealand.
Three-Phase A.C. -D.C. Load
Flow
7.1. INTRODUCTION
Any converter which is operatin g from an unbala.nced a.c. system will itself
operate with unbalanced power flows and unsymmetric valve conduction
periods. In addition any unbalance present in the converter control equipment or
any asymmetry in the converter transformer will introduce additional unbalance.
Considerable interaction exists between the unbalanced operation of the a.c.
and d.c. systems. The exact nature of this interaction depends on features such as
the converter transformer connection and the converter firing controller.
High-power converters often operate in systems of relatively low short-circuit
ratios where unbalance effects are more likely to be significant and require
additional consideration. The steady-state unbalance and its effect in converter
harmonic current generation may also influence the need for transmission line
transpositions and the means of reactive power compensation.
The converter model for unbalanced analysis is considerably more complex
than those developed for the balanced case in Chapters 3 and 6. The additional
complexity arises from the need to include the effect of the three-phase converter
transformer connection and of the different converter firing control modes. Early
h.v.d.c. control schemes were based on phase angle control, where the firing of
each valve is timed individually with respect to the appropriate crossing of the
phase voltages. This control scheme has proved susceptible to harmonic stability
problems when operating from weak a.c. systems. An alternative control, based
on equidistant firings on the steady state, is generally accepted to provide more
stable operation, ( 1 ) . ( 2 ) ~ ( 3 ) Under normal steady-state and perfectly-balanced
operating conditions, there is no difference between these two basic control
strategies. However, their effect on the a.c. system and d.c. voltage and current
waveshapes during normal, but not balanced, operation, is quite different. A
three-phase converter model must be capable of representing the alternative
control strategies.
Similar techniques are available for the integration of the three-phase
converter model into the load-flow analysis as were discussed with respect to the
balanced single-phase analysis. Based upon the extensive investigations into the
behaviour of single-phase a.c.-d.c. load flow described in Chapter 6, the
sequential approach is considered the most appropriate for the integration of the
d.c. model into the three-phase load flow. The complexity of the unified approach
is not considered justified in the three-phase case because, in cases of difficult
convergence such as those involving very weak a.c. systems, it is possible to use
starting values derived from single-phase analysis.
This chapter describes the development of a model for the unbalanced
converter and its sequential integration with the three-phase fast decoupled load
flow described in Chapter 5.
where :
are vectors of the balanced internal voltages at the generator
internal busbars;
pH are vectors of the three-phase voltages at every generator
terminal busbar and every load busbar;
-
x is a vector of the d.c. variables (as yet, undefined).
The significance of the three-phase a.c. variables was discussed in Chapter 5
and the selection of d.c. variables iis discussed in Section 7.3.
To enable a Newton-Raphson based technique to be used, it is necessary to
formulate a set of n independent equations in terms of the n variables describing
the system. As explained in Chapter 5, the equations which relate to the a.c.
system variables are derived from the specified a.c. system operating conditions.
The only modification to these equations, which results from the presence of the
d.c. system, occurs at the converter terminal busbars. These equations become:
where Pferm(dc)
and Q,Perm(dc)
are function of the a.c. terminal conditions and the
converter variables, i.e.
The equations for the a.c. system may therefore be summarized as,
where the mismatches at the converter terminal busbars are indicated separately.
Further equations are derived from the d.c. system conditions.
That is, for each converter k a set of equations
Primary Secondary
equivalent single bridge. The dimensions of the three-phase d.c. model, will,
therefore, be much greater than the balanced d.c. model.
All converters, whether rectifying or inverting, are represented by the same
model (Fig. 7.2) and their equations are of the same form.
Basic assumptions
To enable the formulation of equation (7.2.6) and to simplify the selection of
variables 2 the following assumptions are made:
(i) The three a.c. phase voltages at the terminal busbar are sinusoidal.
(ii) The direct voltage and direct current are smooth.
(iii) The converter transformer is lossless and the magnetizing admittance is
ignored.
Assumptions (ii) and (iii) are equally valid for unbalanced three-phase analysis
as for single-phase analysis. Assumption (i) is commonly used in unbalanced
converter studies (4)*(5) and appears to be backed from the experience of existing
schemes. However, a general justification will require more critical examination
of the problem.
Under balanced operation only characteristic harmonics are produced and, as
filtering is normally provided at these frequencies, the level of harmonic voltages
will be small. However, under even small amounts of unbalance, significant
noncharacteristic harmonics may be produced and the voltage harmonic
distortion at the terminal busbars will increase. The possible influence of
harmonic voltages at the converter terminal, on the fundamental frequency
power flows to the converter, is considered in Section 7.7.
Phase1
Phase3 Phase 2
7 Phase 1
1
Phase 2
D.C. voltage
The d.c. voltage, found by integration of the waveforms in Fig. 7.3(ii), may be
expressed in the form:
Id, - Id3 = 0
Id, - Id, =0
For a six-pulse unit, the interval between firing pulses in specified as 60". This
provides two more equations.
In the equation above, phase (i) is selected during the solution procedure such
that the other two phases will have, in the unbalanced case, firing angles greater
than amin.
With conventional phase angle control, the firing angle on each phase is
specified as being equal to a,,,, i.e.
a, - amin= O (7.3.6)
a, - a,,, =0 (7.3.7)
a j - a,,, = 0 (7.3.8)
The remaining three-control equations required are derived from the operat-
ing conditions. Usually, the off-nominal taps are specified as being equal, i.e.
The final equation will normally relate to the constant current or constant
power controller, e.g.
Similar equations apply to the other two phases with a cyclic change of suffices.
"ode = mi
The three-phase converter transformer is represented by its nodal admittance
By subtracting #, in the above equation, the terminal busbar angles are related
to the converter angle reference.
Separating this equation into real and imaginary components, the following six
equations result:
3
Iicosoi = 1 [bit E, sin 4, + bit V:,,
k= 1
sin (Ofer, - O:,,)] (7.3.16)
Iisin oi= 1 [ - b:
k= 1
E, cos 4, - b;t Vf,,, cos (Of,,, - O:,,)] (7.3.17)
Three further equations are derived from approximate expressions for the
fundamental r.m.s. components of the line current waveforms as shown in
Fig. 7.3, i.e.
4.1,
li= 0.995 -sin (Ti/2) (7.3.18)
Jr
where Ti is the assumed conduction period for phase i.
The sum of the real powers on the three phases of the transformer secondary
may be equated to the total d.c. power, i.e.
3
C E,.I,-cos(q5i - a,)- Vd.Id= 0 (7.3.19)
i=l
The derivation of the last two equations is influenced by the position of the
fundamental frequency voltage reference for the secondary of the converter
transformer.
The voltage reference for the a.c. system is earth, while in d.c. transmission the
actual earth is placed on one of the converter d.c. terminals. This point is used as a
reference to define the d.c. transmission voltages and the insulation levels of the
converter transformer secondary windings.
In load-flow analysis, it is possible to use arbitrary references for each
converter unit to simplify the mathematical model. The actual voltages to earth, if
required, can then be obtained from knowledge of the particular configuration
and earthing arrangements.
With a star-connected secondary winding an obvious reference is the star point
itself. If the nodal admittance matrix is formed for a star-g-star-g connection then
this reference is implicitly present through the admittance model of the
transformer. In this case, however, the converter transformer does not restrict
the flow of zero-sequence currents and the following two equations may be
written :
3
C I,lo,=O (7.3.20)
i= 1
These two equations (real and imaginary parts) complete the set of twelve
independent equations in terms of the twelve additional variables.
However, the above considerations do not apply to delta-connected secondary
windings.
To obtain a reference which may be applied to all transformer secondary
windings, an artificial reference node is created corresponding to the position of
the zero-sequence secondary voltage. This choice of reference results in the
following two equations:
3
E , c o s ~= O~ (7.3.21)
i= 1
3
R(2) = Eisin 4i= 0
i= 1
4 I d sin ( T 2 / 2 )
R(5) = I, - -.---
nfi
3
R(9) = 13.cosw 3 - 1 [b: Eksin 4,+ b,3,k.VFermsin (O:,,, - 8:e,m)]
k= 1
3
R(11) = I,.sin w 2 x 1 [b:: E kcos 4,+ b:;. v:,, cos (qerm
- O;,,)]
k= 1
where [B'] and [ B ] are the three-phase fast-decoupled a.c. Jacobian matrices as
developed in Chapter 5, and [J] is the d.c. Jacobian of first order partial
derivatives.
Equations (7.4.1) and (7.4.2) are the three-phase fast-decoupled algorithmic
equations from Chapter 5. For the solution of the equations (7.4.1) and (7.4.2), the
d.c. variables 2 are treated as constants and, in effect, the d.c. system is modelled
simply by the appropriate real and reactive power injections at the converter
terminal busbar.
These power injections are calculated from the latest solution of the d.c. system
equations and are used to form the corresponding real and reactive power
mismatches. For the d.c. iteration, the a.c. variables at the terminal busbars are
considered to be constant.
The iteration sequence for the solution of equations (7.4.1), (7.4.2) and (7.4.3) is
illustrated in Fig. 7.5. It is based on the P, Q, DC sequence described in Chapter 6
which proved the most successful sequential technique in the single-phase case.
This sequence acknowledges the fact that the converter operation is strongly
related to the magnitude of the terminal voltages and more weakly dependent on
their phase angles. Therefore, the converter solution follows the update of the a.c.
1Evaluate real power mismatches]
1
"' Evaluate d.c. residuals
I
t
Solve equation ( 7.4.3 )
and update i
1 I
kP=O
terminal voltages. It should be noted, however, that in the three-phase case, final
convergence is comparatively slow because the d.c. system behaviour is
dependent on the phase-angle unbalance as much as on the voltage unbalance.
-initialize d.c.
variables (410)
Iterative solution
procedure. (Fig. 7.5)
(1325)
Fig. 7.6.
By using row pivoting only during the solution procedure, the block sparsity of
Fig. 7.7 is preserved. Each block containing nonzero elements is stored in full, but
only nonzero elements are processed.
This routine requires less storage than a normal sparsity program for
nonsymmetrical matrices and the solution efficiency is improved.
Test system
The performance of the algorithm is discussed with reference to the test system
illustrated in Fig. 7.8. The system consists of two a.c. systems interconnected by a
600 kV, 600 MW h.v.d.c. link.
The twenty-bus system is a representation of the 220 kVa.c. network of the
South Island of New Zealand. It includes mutually-coupled parallel lines,
Twenty-bus system Five - bu; system
(a)
BUS 01
I
BUS 05
Fig. 7.8. Three-phase a.c.-d.c. test system. (a) H.V.D.C. interconnection;(b) Five-bus a.c.
system
Generator data
Sequence
reactances Voltage
Power regulator
Name Xo x I *, (MW) Va
Busbar loadings
1 2 3 4lteration(b) 1 2 3 4 Iteration
(vi) As for case Mi); with large unbalanced load at BUS03 8,7 7,6
(vii) As for case b(ii); with large unbalanced load at BUS03
(viii) As for case Hi); with loss of 1 line BUS01 to BUS03
(ix) Symmetrical firing; ai = a,,,,a, = - 10%,a, = O,a, = + 10% 7,6
(XI Phase-angle control; a , =a, =a, = d p ,a , = a, = a,, p,, = P z 8,7 7,6
slowed are (viii) and (xi) where the system is weakened by the loss of one
transmission line. This is to be expected from the discussion of single-phase
sequential algorithms given in Section 6.8.
To examine the effect of a weak system in the three-phase case, the convergence
patterns for the terminal powers and voltages are shown for case (xi) in Fig. 7.10.
The reactive power and voltage unbalance vary considerably over the first few
iterations but this initial variation does not cause any convergence problems.
With weaker systems, the unbalance increases and the convergence patterns
become more oscillatory. The corresponding convergence pattern of the single-
phase load flow for case (xi) is shown in Fig. 7.1 l(b) where a similar oscillatory
pattern is observable. Moreover the sum of the three-phase reactive powers and
x-X-X-X %rm
I (tve sequence)
the average phase voltage of Fig. 7.10, plotted in Fig. 7.1 l(a), shows an even
closer similarity between the three-phase case and the single-phase behaviour.
Sample results
.The operating states of the two converters connected to BUS03 are listed for the
most typical cases in Table 7.3. The corresponding a.c. system voltage profiles
and generation results are given for cases a(i), b(i) and b(ii) in Table 7.4. The
following discussion is with reference to these results.
The results of the realistic three-phase converter model (case b(i)), although
Table 7.3(a) Converter 1 results
Converter 1 (Star-Star)
Converter 2 (Star-G-Delta)
Commu-
Firing tation
angle angle Real Reactive Voltage Current
Case Phase a, (deg) ui (deg) Pi (MW) Qi (MVAr) V d , (kV) I d ,
7.00 29.80
7.00 29.60
7.00 29.32
8.03 28.97
7.00 29.57
8.55 28.08
7.00 30.63
7.00 28.92
7.00 28.90
7.00 30.48
14.95 23.25
13.41 24.25
8.08 25.42
8.38 27.30
7.00 26.96
Table 7.4. Bus voltages and generation results
Case a(i)
Case Mi)
Case Mii)
distinguishable from those of the balanced model a(i), are not significantly
different as regards the a.c. system operation. They are definitely significant
however, as regards converter operation, particularly when consideration is
given to the harmonic content.
A comparison of cases b(i) and b(ii) shows an increase in reactive power
consumption in case b(ii) due to two phases having greater than minimum firing
angles.
The results also show that the transformer connection modifies the converter
source voltages and the phase distribution of power flows. Under balanced
conditions, a zero-sequence voltage may appear at system busbars. As the
converter has no zero-sequence path, zero sequence current will only flow when
the converter transformer provides a path, as in the case of the star-G/delta
transformer. A typical example is illustrated in Fig. 7.12 where the zero-sequence
voltages and currents are shown for case b(i). Accurate converter transformer
models must therefore be included in the converter modelling.
a
(degrees) Firing angle errors due to shift of zero crossings
Errors in calculated phase current magnitudes
Appreciable errors occur with the case of phase-angle control; errors with
symmetrical firing are limited to the effects of commutation angle unbalance.
Theeffect of up to two degrees modulation has been investigated on the basis of
the current waveform of a converter with commutation angle of ten degrees and
balanced voltages. Ignoring the shifts in commutation angle which will occur
with an alteration of firing angle, the effect of alteration in the period of
conduction has been investigated using a Fourier Transform algorithm. The
results are shown in Table 7.6. In addition to the fundamental, the percentages of
the harmonics are also given. Note that no even harmonics are present as the
waveform was assumed to be symmetrical.
It is important to note that the percentage errors presented in the Table are the
maximum that can occur; in all practical cases there will be an alteration in the
commutation angle which will inhibit the change in the waveform. This effect will
be most noticeable at small firing angles. In general therefore, the errors will be
less than those indicated in the table.
The change in fundamental is not, therefore, expected to exceed 1 percent of the
value calculated by the load flow.
7.8. REFERENCES
1. J. D. Ainsworth, 1967. 'Harmonic instability between controlled static converters and
a.c. networks', Proc. IEE, 114 (7), 949-957.
2. J. D. Ainsworth, 1968. 'The phase-locked oscillator-A new control system for
controlled static converters', I EE E, PAS, PAS-87 (3), 859-865.
3. J. Kauferle, R. Mey and Y. Rogowsky, 1970. 'H.V.D.C. stations connected to weak a.c.
systems', I EEE, PAS, PAS49 (7), 1610- 1617.
4. A. G. Phadke and J. H. Harlow, 1966. 'Unbalanced Converter Operation', IEEE
PAS, P A S S , 233-239.
5. J. Arrillaga and A. E. Efthymiadis, 1968. 'Simulation of converter performance under
unbalanced conditions', Proc. IE E, 115 (12), 1809-1 818.
Faulted System Studies
8.1. INTRODUCTION
The main object of Fault Analysis is the calculation of fault currents and voltages
for the determination of circuit-breaker capacity and protective relays
performance.
Early methods used in the calculation of fault levels involved the following
approximations:
- All voltage sources assumed a one per unit magnitude and zero relative phase,
which is equivalent to neglecting the prefault load current contribution.
- Transmission plant components included only inductive parameters.
- Transmission line shunt capacitance and transformer magnetizing imped-
ance were ignored.
Based on the above assumptions, simple equivalent sequence impedance
networks were calculated and these were interconnected according to the fault
specification. Conventional circuit analysis was then used to calculate the
sequence voltage and currents and with them, by means of the inverse sequence
component transformation, the phase components.
Although the basic procedure of the computer solution is still the same, the
need for the various approximations has disappeared.
The three-phase models of transmission plant developed in Chapter 2, which
included interphase and parallel line mutual effects, could be easily combined to
produce the faulted system matrix admittance or matrix impedance and hence
provide an accurate model for the analysis of a.c. system faults.
However, the main reasons given for the use of the phase frame of reference in
load flows are less relevant here. Extra losses and harmonic content are less of a
problem in the short period of time prior to fault clearance. Fault studies are
normally performed on systems reasonably well balanced either at the oper-
ational or planning stage; in the latter case only after prospective system
configurations have been proved acceptable through load-flow studies.
Moreover, unbalanced faults involving large converter plant, while in need of a
three-phase representation, cannot be assessed by steady-state models as
explained in Chapter 11.
Finally, faulted system studies constitute an integral part of multi-machine
transient stability programs, the complexity of which will not normally permit
the three-phase approach.
A single-phase representation, achieved with the help of the symmetrical
components transformation(') is used in this chapter as a basis for the
development of a fault-study p r ~ g r a m . ( ~ ) . ( ~ ) . ( ~ ) ~ ( ~ )
When analysing the first two or three cycles following the fault, the
subtransient admittance of the machine is normally used, whilst for longer times,
it is more appropriate to use the transient admittance. The machine model,
illustrated in Fig. 8.l(a), is then converted to a nodal equivalent by means of
Norton's Theorem which changes the voltage source into a current source
injected at the bus j as shown in Fig. 8.l(b). This is most effective as otherwise a
further node at j' is necessary to define the machine admittance yM.
The injected nodal current is given by
where
so that
17 is the current required at the voltage V, to produce the machine power
P7 +j Q 7 , so
(I?)* Vj = P y +j - Q r (8.2.4)
Thus from the load-flow data of P M ,Q~ and V Mwe may calculate the injected
nodal current I ; , as
The following equations may then be written for the network of Fig. 8.4.
or in matrix form after grouping together the terms common to each voltage:
Fig. 8.3. Model substitution
where [a and [V] are the current and voltage vectors and [Y] is the nodal
admittance matrix of the system of Fig. 8.2.
It can be seen from equations (8.2.6) to (8.2.10) that nonzero elements only
occur where branches exist between nodes. Since each node or busbar is normally
connected to less than four other nodes, there are usually quite a number of zero
elements in any system with more than ten busbars. Such sparsity is exploited by
only storing and processing the nonzero elements. Moreover, the symmetry of the
matrix (Yij = Yji) permits using only the upper right hand terms in the
calculations.
[v] = [Y]-~.[I]
=P I . [I1
This equation uses the bus nodal impedance matrix [Z] and permits using the
Thevenin Equivalent circuit as illustrated in Fig. 8.5, which, as will be shown
later, provides a direct solution of the fault conditions at any node. However, the
use of conventional matrix inversion techniques results in an impedance matrix
with nonzero terms in every position Zij.
The sparsity of the [Y] matrix may be retained by using an efficient inversion
and the nodal impedance matrix can then be calculated directly
from the factorized admittance matrix.
Fault calculations
From the initial machine data, the values of [aare first calculated from equation
(8.2.5) using one per unit voltages. These may now be used to obtain a better
estimate of [V], the prefault voltage at every node from equation (8.2.13). If the
initial data is supplied from a load flow, this calculation will not make any
difference.
The program now has sufficient information to calculate the voltages and
currents during a fault.
The voltage at the fault bus k is from Fig. 8.6.
I/{ = z * . I*
where
k is the bus to be faulted,
Z* is the fault impedance and
I* is the fault current.
From equations (8.2.19) and (8.2.22) the fault voltages at every bus in the
system may be calculated, each calculation requiring only one column of the
impedance matrix. Using the bifactorization routines the kth column can be
obtained by multiplying the impedance matrix by a vector which has a '1' in the
kth row and '0"s elsewhere, i.e.
where
Admittance matrices
The data specifying each element of the system is then used to form the following
three nodal equations:
OIi = O V i O y i i + ( O V i + O V l ) O y .i .i .++ ( O V i - O V n ) O ~ n i (8.3.4)
' I i = ' V i ' ~ i i + ( ' V i - ' V l ) ' ~ i i -+. . + ( ' V i - ' V n ) ' ~ n i (8.3.5)
2 1 i = 2 ~ i 2 y i i + ( 2 V 2i -~ l ) 2 y i i + . , . + ( 2 V 2Vn)2yni
i- (8.3.6)
where
' I i is the zero-sequence injected current at bus i.
' V i is the positive-sequence voltage at bus i.
and
2yniis the negative-sequence admittance between nodes n and i.
The above equations can be expressed as:
[OI] = [O Y ] [OV]
[ ' I ] = ['Y ] [' V ]
C21] = C2 Y ] C2 V ]
where
and
n
YYii = C Yyik for i = 1, n; y = 0,l, or 2
k= 1
Fault calculations
As already explained for the three-phase fault, the nodal impedance matrices may
now be calculated directly from the reduced admittance matrices and the
following sequence impedance matrix equations result.
[O V] = [OZ] [Ol]
Because the system is assumed to be balanced prior to the fault, the vectors of
negative and zero-sequencecurrents are zero, i.e. there are no prefault negative or
zero-sequence voltages.
The positive-sequence network then models the prefault network condition
and equation (8.3.11) is used to calculate the prefault voltages. If the original
voltages used in the machine models were obtained from a load-flow calculation,
then the use of equation (8.3.1 1)will make no difference to those results; however,
if the voltages were assumed at one p.u. with zero angle then this calculation will
provide more accurate prefault voltages.
The single-phaseequivalent circuit is then set up by linking the three sequence
networks together according to the type of fault to be analysed.(*)
Short-circuit faults
A convenient way of simulating the fault location F for the analysis of short-
circuit faults is illustrated in Fig. 8.8.
It includes three fault impedances a Z , b and~ 'Z and three injected currents
'I*, bl*, and 'I*
For each type of fault, it is possible to write 'boundary conditions' for the
currents and voltages at the fault location. For example, Fig. 8.9 shows the case
of a line to ground fault at bus k.
0
A
-
a
b .r
C
A
4
Ov bv cv
-
I
Fig. 8.8. The fault location
0
A -
b v
a
A
C a
A
c'., olf'r b;fJf; 0
"% *'.,
-
1
-
Fig. 8.9. Single line to ground fault
Also, the sequence voltages at the fault location may be described by the
equations
Ovf -
k -
-0 ~ ~ ~ ' O l f
(8.3.17)
lv,f =lvk- lzkk.llf (8.3.18)
2 v f= - 2 z k k . 2 1 f (8.3.19)
From equations (8.3.15) to (8.3.19), the following relationships are obtained.
Similar considerations yield the fault currents for other types of short-circuit
fault. The results for line-to-ground, line-to-line, line-to-line-to-ground, and line-
to-line-to-line faults are illustrated in Table 8.1.
Table 8.1. Fault currents for short-circuit faults
L-G 'If
- '1f
K - OZfi.' I f - 2Zfi.'If
G L-G
+
( 2 z ; i . 0 z ; ) / ( 2 z ; i OZ;) + 'ZIi 2Zii + OZIi 2zIi + OZIi
L-L-GG
where
+
'Zli = ' Z i i 0.5 Z*
'Z,',= ' Z i i + 0.5 Z*
OZIi = OZii + 0.5 Z f
These fault currents at the fault location are then added to the current vectors
['I], ['I] and ['I] to produce the fault current vectors [OI*], [ l I*] and ['I*]. For
a fault at bus k these are
0 for i # k
-'I* for i = k
'1-I =
{::: for i # k
- lI* for i = k
0 for i # k
- 'I* for i = k
The fault voltages are then obtained from equations (8.3.10) to (8.3.12) by
substituting the fault current vector for the prefault current vector, i.e.
[O V*] = [OZ] (8.3.24)
[l V * ] = [ l Z ] [ l l * ] (8.3.25)
[' V * ] = ['Z] ['I*] (8.3.26)
Open-circuit faults
The system is now represented by a two-port network across which the faulty line
is connected as shown in Fig. 8.10. In this case, the prefault voltages have to be
obtained from a load-flow study.
System
Oz * 01'
q, * dl'
A
A
A
Cz CIf
A
+ A
ck$ *k$ '5 dl(
-
-L
Fig. 8.10. Two-part network with faulty line where " Z ,bZ
or 'Z may be on open circuit
For an open-circuit fault on phases '6' and 'c' the boundary conditions are
Equations (8.3.27) and (8.3.28) define the connection of the Thevenin equiva-
lent sequence networks at the fault location to solve for the fault currents.
The equivalent Thevenin impedances are the sequence impedances of the
system between the two buses k and 1, i.e.
and the equivalent Thevenin voltage is given by the difference between the
voltage at buses 1 and k with the faulted line disconnected.
During the fault the sequence voltages OV,/,,'V{~and 'V,/, have the same
expressions as equations (8.3.17) (8.3.18) (8.3.19).Thus similar considerations, as
in the case of the line-to-ground short-circuit, lead to the following expression for
the fault currents.
The case of a single open-circuit fault can be analysed in a similar manner and
the final relevant equations are shown in Table 8.2.
The fault current vector is formed as follows:
0 for i = l , n i#korl
i=l
'Ii for i = l,n if korl
'Ik - f i =k
I ,I i= 1
ONEO-C
TWOO-C
where
'z = 'zkk+ 'zll- 'zk l - l Z l k
Oz = Oz,, + Oz,,- Oz,, - Oz,,
2~ = 'Zkk + 'Zll - 2zkl- 'Zlk
zf is the sum of the positive, negative and zero sequence impedances of the
faulty circuit, i.e.
zf = 'zf + 2zf + Ozf
'Z' = 'z + 'zf
OZ' = oz + Ozf
2 ~ =' 2z + 2zf
From the fault voltages the branch currents are obtained as follows:
O I f . = OY. .(O v{ - O V
)!
11 IJ (8.3.31)
'16 = yij(l V{ - v!) (8.3.32)
21; = zYij(2I/{ - vf) (8.3.33)
Where necessary the corrections for taps on the positive and negative sequence
networks are
I
II Calculate the fault currents
and print fault MVA
A Next fault
..
UIN
oa
*N
00
N9
me
..
d
00
-..
00
09
on
Cul
0-
00
N9
..
09
rO
00
00
00
00
00
..
00
0-
00
00
2s
..
-a
00
00
00
3"i
em
'2
..
00
%
W
O
a
J
.
C w
..................................
00000
00000
*0000
00000
0a000
00000
00000
00000
00000
00000
OOOOO
00000
0000
00-
OONO
0000
....
m-0
9010
0000
0000
N- 4-0
11)9U t U 0000
.....
0011)CC 0000
cnon* UUOO
OOOOO
O W 0 0
?Y.9
0000
ONbsO 00111)
00-C 4rlll11)
.....
DW--
0-,
roo00
00000
....
CC11)N
2!2zg
0000
Z
*
0
4 0 -4
OI+.IvoI? OC-
O O M I OS9.2-
44, I ? C O )
QOOOO QOO00 001)Oo
-d
*S-
02s-
-99-
*O.rD
O S I O O -DOD
ul*OFIN O O I L n I I - I(YOIII)
0-
x ...................a
n
C
Z
W
a
0
3
3
c
W
xx-
21-2
W
ww-=*
mats+-
0
C
Szaa IS
0333 00
z x z d s A Y ~ S D ax
S 3 a a I Z d3I15 3 3 - S
N'PI~ mnaa-
-@322 * - 1 S 1
Sooaa om--
a
-=co
cainxx arccx -c*a
t-m
*I 0
003 %
X
13.uId3
~ N Y ' O
e
-SS-
a - =
rc-
~a a ~ ~ ~
a m d a .&41'&T.
-2-
-a-
=ws&
w-gc~
cch
5
a- Ja
-
J-
'&r..C
$
ZS
m 9 r o o 1-0-
**Coo no-
r-noo - = w o o 0
aPB-00 no-
no000
c-
-9000
wmn
M-C-0
000ac
001m
.....
ocom
OOONO
oowo
IIIII
q OIDW
-nsm
n n n o e n r r r a r -QYI
8-9 NIIO~ s909
gcroo
nrnro n m o o o CIO~IYN
oe-0- n
- etnal-0
lor-
L ulm- o n m u r u o n e -un 4-9 uaerotu ON- oOr09
,
n
5
omo-
000,
-8000
0-0
O m 0 4 OOCOO
,000 ,NO0 o v;;c;,'f d;&G
==Y!!!!!!!!xE?~
&AS ;,'s;
g 42-5
00000 00000
t I
:?
' 2:4
II
439::
OOOOO
II
33z:
00000
I
,Z
0
/
n
0
C
C
3
*
2
x
n
I
nrr-
000.-
8 ss494 S?.SS
n 00000
N I 1
o
-NN
-
n n s v ' n -09--
-e
0000-
~s r o c aaoor
9.-D-
-I0
o 3 w o ooocuo 00000 00000
~
r N N o r moo90
1
O N 0 4
sssss sssss
00000 00000
I 1 Ill
00000
1
~
o
5
0
nrrmn
cV-)r,
4
c-.
I
-
-OD-
-,3r
p
11
oeono
*n.nrw
%?3S
%SSZ
-.*.*.
r rrod e w e
'7 ' 7
a
1
b
3 P
0 .c
4
f
s3
.
f
(I)
g
N
n
-
Z
rl
3
(1
rl
II)
m
3 I W r
Ulw
Z
0
>
x
f w
z L%tS 5.
a a
" O I -
023 3:=
C.2
9
k:&4gp2
3 -203
>-,t
-2-
s 4 4 332%
a-Ze
H-,z
-ax0
2
:
s
+
03=HA
=A= C*C.
& E 4 3 gg&
CC01oQ
.....
*CON
OIOOO
QINOI
O*N N
4 4 4
d
where p denotes the differential operator -.
dt
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
which 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 degree (electrical) from each pole and is referred to as the 'quadrature axis'.
In the dynamic state, each synchronous machine is rotating independently
from each other 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 which rotates at synchronous speed.
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.
Mechanical equations
The mechanical equations of a synchronous machine are very well estab-
li~hed(').(~)
and need be only briefly outlined. Three basic assumptions are made
in deriving the equations:
Machine rotor speed does not vary greatly from synchronous speed
(1.0 P.u.).
Machine rotational power losses due to windage and fric'tion are ignored.
Mechanical shaft power is smooth, that is 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:
Electrical equations
The derivation of equations to account for flux changes in a synchronous
machine has been given by C ~ n c o r d i a 'and
~ ) Kimba~-k.(~)
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 Chapter 10.
(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.
The per unit system adopted is normalized to eliminate factors of &a, n
and turns ratio, although the term 'proportional' should be used instead of 'equal'
when comparing quantities. Note that one 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 9.1 shows the flux and voltage phasor diagram for a cylindrical rotor
synchronous machine in which all saturation effects are ignored. The flux Ff
Direct axis
Quadrature
axls
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 degrees. The terminal voltage V is found from this
voltage El by considering the volt drops due to the leakage reactance XI and
armature resistance Ra. By similar triangles, the difference between Ef and El is in
phase with the I X l volt drop and is proportional to I . Therefore the voltage
difference may be treated as a volt drop across an armature reactance X a . The
sum of X a and XI 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 axes it is
convenient to resolve many phasor quantities into components in these axes.
The stator current may be treated in this manner. Although F , will be
proportional to I, 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 I , and F , is
unity but the proportionality between I, and F, is less than unity and is a function
of the saliency.
In Fig. 9.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
Direct axis
\A ,,
Quadrature
f axis
Fig. 9.2. Phasor diagram of a salient pole synchronous machine in the steady state
The phasor diagram of the machine operating in the transient state is shown in
Fig. 9.3.
\
Fig. 9.3. Phasor diagram of a synchronous machine in the transient state
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:
E: - V, = RaI, - X,"Id (9.2.15)
E,"- Vd = RaId + X i I , (9.2.16)
pEi = (Ei + (XI;- X i ) I d - Ei)/T;, (9.2.17)
pE," = (El;- ( X i - X,")I, - E:)/Tlj, (9.2.18)
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 Fig. 9.4. It
Fig. 9.4. Phasor diagram of a synchronous machine in the subtransient state
should be noted that equations (9.2.11) and (9.2.12) are now true only in the
steady-state mode of operation, although once subtransient effects have decayed,
the error will be small.
MACHINE MODELS
It is possible to extend the model beyond subtransient level but this is seldom
done in multi-machine programs. investigation^(^) 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 (Xi)
requiring no differential equations. Only the algebraic equations (9.2.11) and
(9.2.12) are used.
Model 2-D-axis transient effects requiring one differential equation (pEh).
Equations (9.2.1I), (9.2.12) and (9.2.13) are used.
Model 3-0- and q-axis transient effects requiring two differential equations
(pEb and pE;). Equations (9.2.1I), (9.2.12), (9.2.13) and (9.2.14) are used.
Model 4--0- and q-axis subtransient effects requiring three differential
equations (pEi, pE," and pEi). Equations (9.2.13), (9.2.15), (9.2.16), (9.2.17) and:
are used. This last equation is merely equation (9.2.14) with modified primes.
Whether it is a subtransient or transient equation is open to argument.
Model 5-D- and q-axis subtransient effects requiring four differential
equations (pE;, pE;, pEi and pE;). Equations (9.2.13),(9.2.14), (9.2.15), (9.2.16),
(9.2.17) and (9.2.18) are used.
The mechanical equations (9.2.5) and (9.2.6) are also required to be solved for
all these models.
Groups of synchronous'machinesor 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 (9.2.5) and (9.2.6) are not
required. This sixth model is thus defined as:
Model 0-Infinite machine-constant voltage (phase and magnitude) behind d-
axis transient reactance (Xi). Only equations (9.2.1 l ) and (9.2.12) are used.
Saturating
Regulator amplifter exciter
1 Ef
I
Other
sig&ls
- Kfp
(1+T/.p)
Feedback loop
voltage
Controlled
voltage
Other
signals Feedback loop
Fig. 9.5. Block diagrams for two commonly used AVR models. (6) (a) IEEE Type 1
AVR model; (b) IEEE Type 2 AVR model. (01982 IEEE)
Vt
Co$olled =
1 a,,' + -
Z -- 1
( U ~ T W') T e l d
N
busbar voltage
voltage -
signal
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.
The equations for the AVR model shown in Fig. 9.6 are as follows:
216
subject to
and
where Se =f ( E f )
V g = E f [unless IEEE Type 2 when V g = Val (9.3.10)
Vaux= A predefined signal
The IEEE(6)recommends that Se be specified at maximum field voltage (Se,,)
and at 0.75 of maximum field voltage (Seo.7,,ax).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 (9.3.9)may by transfor-
med to:
V e = ( k , .Ef - k2)Ef (9.3.11 )
where
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.
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 recom-
mended by the IEEE(7)are shown in Fig. 9.7. It will be noticed that if limits are
not exceeded, the two models are identical. The difference is due to the
Specifled Ps
power setting 7
4 v
(a)
Specified
power sett~ng 7"
M,~,,,,
speed power setting
277 f o b Ref. Speed
(b)
Fig. 9.7. Typical models of speed governors and valves. (7) (a) Thermal governor
and valve; (b) Hydro governor and valve. (01982 IEEE)
( 1 t T1.p) ( 1 t T3.p)
Power to
I turbine
assumption that, in a hydro 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 Fig. 9.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 Fig. 9.8 are:
The valvelgate position setting (Gv) is subject to opening and closing rate limits
(omaxand cmX respectively) and to physical travel limits so that:
pw
Power dehvered (1tT4p)
through valve
PI
I
Mechanical power
to generator
/P and L P
(b) turbine power
here Pm, is the initial steady-state mechanical power. Note that for this simple
model, the initial value of Pgv is PmIK 1, even though all the steam passes through
the valve.
Provided that the H P valve does not close fully, then, rather than inject the
power from the I P and L P turbines as shown, it is easier to treat it as a simple
turbine (PI = 0 ) but with the speed regulation modified by:
For the circuit shown in block diagram form in Fig. 9.10, the equation is:
K(1+ T2p)
vout =
+
( 1 T l p ) Vin
This can be transformed to:
+
vout =-(
K . T 2 (TlIT2) T I - p
T1 1 T1.p + ) 'in
and then to:
which can be represented by the block diagram in Fig. 9.1 1, and is a lag circuit in
parallel with a gain.
It is important to remember that the time constant TI, must be nonzero even if
the integration method can accommodate zero-time constants.
9.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 specific large loads which are considered
in the next chapter. 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 character is ti^(^) may be adopted such that the MVA loading at
a particular busbar is a function of voltage (V) and frequency (f):
where K p and Kq are constants which depend upon the nominal value of the
variables P and Q.
Fig. 9.12. Characteristic of different load mo-
dels. (a) Active and reactive power against
voltage; (b) Current against voltage
Load
Filament lamp
Fluorescent lamp
Heater
Induction motor half load
Induction motor full load
Reduction furnace
Aluminium plant
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 nondynamic
nature of the load model. From a purely practical point of view, then the load
characteristics are only 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.
Machine
A lmaglnary axis
direct axis
---- ----- v Machine
quadrature axis
VLW.-qI
This transformation is equally valid for currents as is the reverse transformation:
-sin 6 cos 6
When saliency exists, the values of X i and Xi used in equations (9.2.15) and
(9.2.16) and/or X i and X i used in equations (9.2.11) and (9.2.12) are different.
Therefore, the Norton shunt admittance will have a different value in each axis
and when transformed into the network frame of reference, will have time varying
components. However, a constant admittance can be used, provided that the
injected current is suitably modified to retain the accuracy of the Norton
equivalent.(lO)This approach can be justified by comparing the two circuits of
Fig. 9.14 in which Yt is a time varying admittance, whereas yo is fixed.
(bJ
Fig. 9.14. Method of representing synchronous
machines in the network. (a) Norton equivalent
circuit; (b) Modified equivalent circuit
At any time t, the Norton equivalent of the machine is illustrated in Fig. 9.14(a),
but the use of a fixed admittance results in the modified circuit of Fig. 9.1qb).
The machine current is:
and hence
where iadj accounts for the fact that the apparent current source is not accurate in
this case.
The injected current into the network which includes Yo is given by:
where
- -
Zunadj = YOEn
A suitable value for yo is found by using the mean of direct and quadrature
admittances, i.e.
where
Xdq= +(Xi + Xi)
The unadjusted value of current injected into the busbar is:
(9.6.6)
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 (9.6.3)
and also equations (9.2.15) and (9.2.16):
and transforming
The current injected into the network thus represents the deviation of the load
characteristic from an impedance characteristic.
FAULTS
Although faults can occur anywhere in the system, it is much easier com-
putationally to apply a fault to a busbar. In this case, only the shunt admittance at
the busbar need be changed, that is, 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 suitable 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
in infinite admittance in parallel with the existing shunt admittance. In practice, a
nonzero but sufficiently low-fault impedance is used so that the busbar voltage is
effectively 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 be modified.
Alternatively the factors can remain constant and diakoptical techniques(") 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 techniques. 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
interconnection 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 of sufficiently high impedance
that the power 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.
Where a fault occurs on a branch but very close to a busbar, nonunit 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 bareaker 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
voltages 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 1/100, the change
in branch impedance is very small compared to 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
9.7. INTEGRATION
Many integration methods have been applied to the power-system transient
stability problem and the principal methods are discussed in Appendix I. 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 studies. 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. It has only been with the introduction of
more detailed system component models with very small time constants that the
problems of stability 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 (1.4.1) to (1.4.3), for this method
the number of function substitutions is four ( v = 4) and:
The characteristic root of this fourth-order method, when applied to equation
(1.3.3), is:
and to ensure stability, the step length h must be sufficiently small that z, is less
than unity.
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(10),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
(1.3.2) with k equal to unity and is thus a single-step method. The solution at the
+
end of n 1 steps is given by:
It has second-order accuracy with the major term in the truncation error
being -&h3.
The characteristic root when applied to equation (1.3.3) is:
Z, = 1 -2bn+,
where
If Re(R)<O then 0 5 bn+,I 1.0 and lz, 1 I 1.0. The trapezoidal method is
therefore A-stable, a property which is shown in Appendix I 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 Z-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 1-
stable finite difference method possible.
The method, as expressed by equation (9.7.4), is implicit and requires an
iterative solution. However, the solution can be made direct by incorporating the
differential equations into equation (9.7.4). Rearranging forms algebraic equa-
tions as described in Appendix I.
For example, consider the trivial transfer function shown in Fig. 9.16. The
differential equation for this system is given by:
with the input variable being denoted by 'z' to indicate that it may be either
integrable or nonintegrable.
The algebraic form of equation (9.7.7) has a solution at the end of the (n + 1)th
step of:
where
and
Table 9.2
Table 9.3
Despite the three orders of accuracy difference between it and the Runge-Kutta
Gill, the backward Euler method compares well.
The results for the trapezoidal and backward Euler methods were obtained
using linear extrapolation of the nonintegrable variables at the beginning of each
step. This required the storing of machine terminal voltages and currents together
with other nonintegrable variables obtained at the end of the previous step.
For small step lengths the characteristic roots of both methods tend towards, but
never exceed, unity (positive).
ZlcBE) and z,(,,,~) + +1 as h l + 0
The effect of too large a step length can be shown in a trivial but extreme
example. The system in Fig. 9.16 and equation (9.7.7) 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, that is a
constant output of unity. In fact, the output oscillates with y, = 2, y, = 0, y, = 2,
etc.
Table 9.4 shows the effectof different step lengths on this simple system with a
nonzero-time constant T. This table shows that oscillations occur when h l is
smaller than - 2, that is, when the characteristic root z, is negative. The
oscillations decay with a rate dependent on hl, that is, the rate is dependent on the
Table 9.4. The effect of different step lengths on the solution of a simple system (Fig. 9.16) by the
trapezoidal method
where
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'' 2, are available to remove the oscillations but they require a lot of
storage and it is simpler to reduce the step length.
In only a few cases where the differential equation is complex need the value of yn
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 nonintegrable 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 can not
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, has been found to be satisfactory in
almost all cases.(13)
The two mechanical differential equations are given by equations (9.2.5) and
(9.2.6) and the algebraic form is:
where
C , = (1 - 2 . M;Da)o + M,(Pm - Pe + 4nlfo. Da)
and also
where
but machine speed o is a useful piece of information and would still require
evaluation in most problems.
It is also more convenient to retain the electrical power (Pe)as a variable rather
than attempt to reduce it to its constituent parts:
where
Cd = (1 - 2Md)E;- M d ( X q- X;)I,
M d = h/(2Ti0 h) +
also
Ei = C,, + M,,[Ei + ( X ; - XC)ld)
where
C,, = (1 - 2Mqq)E: M,,(Ei+ +( X i-Xi)ld)
M,, = h/(2Tio + h)
and also
E; = Cdd + Md(E;- ( X i - X p , )
where
cdd= (1 - 2Mdd)Ef;+ Mdd(E;- ( X i - X;)I,)
M,, = h/(2Tio + 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 Fig. 9.17, for which can be written:
subject to
The feedback loop can be rearranged to avoid the derivative of Vout being
explicitly required as described in Section 9.3 and this is shown in Fig. 9.18.
Equation (9.7.32) is now replaced by:
Equations (9.7.31) and (9.7.34) can be transformed into the algebraic form:
and
where Ef* is the value of Ef at the previous iteration and k: and kz are
determined from Ef *. The equation describing the saturating exciter is thus:
and applying the trapezoidal rule the algebraic form of solution is:
+
Ef = LCe,- Me,-(Vam)]/[l + Me,-(k: Ef * - kz - Ke,-)] (9.7.40)
where
+ Ke,-)Me,-)Ef + Me,-Vam
Ce,- = (1 - 2(Ke
Me,-= hl(2Te + h(Ke + Ke,-))
Ke,-=k,Ef-k2
Cef, Me,- and Ke,-are evaluated once at the beginning of the step and, hence, only
contain information obtained at the end of the previous step.
Overall structure
An overview of the structure of a transient stability program is given in Fig. 9.19.
Only the main parts of the program have been included, and as can be seen, the
Section 1
I I
1 Read in steady-statesystem data
Section 4
Section 5
I
r
1
1
Etther Solve for machinesand network iteratively
or recaIculatenoninte~rabIevariables
(initially and when iwftching has occurred 1
I
YES
I
1
KBI FA1 = 0
KBIFA3.0
I @%+@ data?
( a 1 Section 1
Switching
data input
and Is
initial conditions
conditions
YES initial
conditions
Perform
numerical part
of, bifactorization
( b ) Section 2
switching ops.
if branch change
set KBIFA3 = 1
if new branch
part of
( c ) Section 3
1 For further
r Print-out
busbar and
branch results
if required
Print-out
sync. machine
results
i
AVR results
if reouired
Print-out
speed gov.
results if
required
7
11
Set flag for step
( d ) Section 4
TIME <MAXTI >
1 Perform oower
balonc; check
to confirm
initial data
1
I
d
,
i End of case'
* KASE =I?
( steady-state
conditions
1
Reset initiol
Isteady-state )
conditions
(e Section 5
Fig. 9.20. Structureof transient stability program (a) Section 1 (b) Section 2
(c) Section 3 (d) Section 4 (e) Section 5
in case H is reduced
C and M for
nonintegrable ---------+c!-~~zi%tGt1
dusually required
L------.J
I
II alaebraic form of
tr~pezoidolmethod
-Same
-Lspeed
for each
gov.
1
I
ITR = 0
IHALF = 0
(a Section 1
ITR = ITR +1
NO
machine
non -integrable ---------~5;Gd;ov_colGot;
--------I usually required I
L ------- J
Evaluate integrable
variables using Same for each A V ~
algebraic form of =I Same for each 1
(
trapezoidal method
also determine
ERROR I+ speed gov.
1
Y
between iterations
YES
( b 1 Section 2
Re-evaluate
conditions
at beginmg of
*
T
0
Not converging'
TIME = MAXTIME
Exit
(c Section 3
Fig. 9.21. Structure of machine and network iterative solution (a) Section l
(b) Section 2 (c) Section 3
If convergence has not been achieved after a specified number of iterations, the
case study is terminated. This is done by setting the integration time equal to the
maximum integration time. The latest results are thus printed out and a new case
study is attempted.
9.10. REFERENCES
1. 0. I. Elgerd, 1971. Electrical Energy Systems 7heory: An Introduction. McGraw-Hill,
New York.
2. B. M. Weedy, 1979. Electric Power Systems. Wiley and Sons, London.
3. C. Concordia, 1951. Synchronous Machines. 7heory and Performance. Wiley and Sons,
New York.
4. E. W. Kimbark, 1956. Power System Stability :Synchronous Machines. (Vol. 3). Wiley
and Sons, New York.
5. P. L. Dandeno, et al. 1973. 'Effects of synchronous machine modeling in large-scale
system studies', I E E E Transactions on Power Apparatus and Systems, PAS92 (2), 574-
582.
6. IEEE Committee Report, 1968. 'Computer representation of exciter systems', IEEE
Transactions on Power Apparatus and Systems, PAS87 (6), 1460-1464.
7. IEEE Committee Report, 1973. 'Dynamic models for steam and hydro turbines in
power-system studies', IEEE Transactions on Power Apparatus and Systems, PAS-92
(6), 1904-1 915.
8. P. L. Dandeno and P. Kundur, 1973. 'A noniterative transient stability program
including the effects of variable load-voltage characteristics', l E E E Transactions on
Power Apparatus and Systems, PAS92 (9, 1478-1484.
9. G. L. Berg, 1973. 'Power system load representation', Proc. I E E 120 (3), 344-348.
10. H. W. Dommell and N. Sato, 1972. 'Fast transient stability solutions', IEEE
Transactions on Power Apparatus and Systems, PAS-91 (4), 1643-1650.
11. A. Brameller, et al. 1969. Practical Diakoptics for Electrical Networks. Chapman and
Hall, London.
12. L. Lapidus and J. H. Seinfeld, 1971. Numerical Solution of Ordinary Differential
Equations. Academic Press, New York.
13. C. P. Arnold, 1976. 'Solutions of the multi-machine power-system stability problem.
Ph.D. thesis, Victoria University of Manchester, U.K.
Power System Stability-
Advanced Component
Modelling
10.1. INTRODUCTION
This chapter develops further some of the component models described in
Chapter 9 and introduces new models needed to investigate the effects of other
a.c. system plant components. Turbine generator models are extended by
considering the effects of saturation in the synchronous machine and the response
of compound thermal turbines. Detailed consideration is also given to the
modelling of induction motors and static power converters. The chapter also
deals with protective gear modelling and unbalanced faults.
The induction motor model allows for a good representation over the whole
speed range so that motor starting can be investigated. The model can be created
in three ways depending upon the induction motor data available.
Basic formulation of three-phase bridge rectification and inversion has already
been considered in Chapter 3 and here it is extended so that the dynamic model
can include abnormal operating conditions encountered during stability studies.
It must be clarified however, that the controllability of h.v.d.c. links during large
disturbances in either the a.c. or d.c. system cannot be determined by the type of
transient stability program described in this chapter. These and other
problems associated with transient stability analysis involving h.v.d.c. links are
considered in Chapter 13.
The grouping of subroutines relevant to a particular component of the power
system or aspect of the study, as developed in Chapter 9, should be retained for
the models produced in this chapter. This ensures that additional models can be
incorporated easily and models removed when not necessary.
Figure 10.2 shows a typical voltage and MMF diagram for a round rotor
synchronous machine. Potier voltage Ep, the voltage behind Potier reactance,
characteristic
Ff
A third method'') distinguishes between the saturation in the rotor and stator,
and saturation factors based on Ep and Ep, are obtained. This method is dficult
to implement because it is necessary to ensure that saturation is not applied twice
to any part of the machine. That is, the saturation in the field poles must be
isolated from that of the armature, giving two saturation curves. The two
saturation factors for this case may be defined as:
iron MMF in the stator
S,=1+
total air gap MMF
iron MMF in the rotor
S,=1+ (10.2.11)
direct axis air gap MMF
where S,, acts equally on both d and q axes and S, acts on the direct axis only.
Figure 10.3, 10.4 and 10.5 demonstrate the differences between the three
methods of saturation representation.
Fig. 10.3. Salient pole synchronous machine with direct axis saturation only
Fig. 10.4. Salient pole synchronous machine with direct and quadrature axis saturation
Fig. 10.5. Salient pole synchronous machine with separate stator and rotor saturation
where n is normally either 7 or 9. Only one point is needed to specify the curve and
if If is always specified at a predetermined voltage, the data entries required per
curve are reduced to one, from which Cn may be readily determined.
Potier reactance
The Potier reactance of a machine is rarely quoted, although the open-circuit
saturation curve is normally available. In order to model the saturation effects, it
is thus necessary to estimate this reactance.
From knowledge of the leakage reactance XI, Be~kwith'~) calculated that:
Xp = XI + 0.63(X; - XI) (10.2.14)
for a salient-pole machine, as most of the saturation occurs in the poles, and
Xp = 0.7X; (10.2.17)
for a round rotor machine, as most of the saturation occurs in the rotor teeth and
the Potier and leakage reactances have similar values.
The effect of saturation on the synchronous machine model
Saturation effectively modifies the ordinary differential equations describing the
behaviour of the voltages used to model the synchronous machine. Equations
(9.2.13),(9.2.14),(9.2.17)and (9.2.18) become respectively:
+
PEI,= ( E f ( X d- X;)Id - SdEI,)/Ti0 (10.2.18)
PEA = ( - ( X , - XI,)I, - SqE;)/TI,o (10.2.19)
~ (SdEi
P E= + ( X , - X,")Id- S,E;)/T,", (10.2.20)
~ (S,E; - (XI, - X;)I, - SqE&')/Ti0
P E= ( 1 0.2.21 )
where S, and S, are the direct and quadrature axes saturation factors.
Where subtransients are considered then equations (9.2.15) and (9.2.16) are
replaced by :
Ei - Vt, = Ra-I, - Xis.Id (10.2.22)
Ei- Vt,=Ra.I,+XC.I, (10.2.23)
That is :
Ra Xdqs 1
-
(Ra2+ X,"X;)
- xaS fi
The current I!$, is similar to iadj
developed in Section 9.6:
The current injected into the network is given by equation (9.6.4) and as the
terms in the brackets contain no saliency then in the real and imaginary axis of the
network:
where i!& contains saturated but nonsalient reactance terms and I$:, contains
salient and saturated reactance terms.
Note that the third part of equation (10.2.27) is yo Vand not yo E". This part of the
injected current is merely the current flowing through yo and could be eliminated
if Yowas not included in the network. The conditioning of the network would be
reduced however, and in certain systems this could lead to numerical problems.
where
and
af i -[(n- ~ ) . C , , ( E ~ , ) " - ~ - X ~-
-- ( XXp)(E;
; - Vt,)]
-1
as, [(sd - l l X p + X&'l
a
- -f 2 - [ ( n - l ) . C n q ( E p d ) n - 2 . X p (-
~ ;Xp)(E," - V t , ) ]
-1 (10.2.34)
as, C(S, - ~ ) X + P xi1
This decouples the Newton method and each saturation factor may be solved
independentl~'~)
sy + 1) , -f 1(p)I(af l I a s d ) c p )
S(P)
d (10.2.35)
sr ) = S:) - f ~ l ( a2/asqyr)
+ f
Despite the advantages of a Newton form of solution, it can be found to be
divergent if too great an integration step length is used. Analysis of the functions
f 1 andf 2 show that they have discontinuities, when XC = (S, - 1)Xp and Xi =
(S, - 1)Xp respectively, although otherwise are almost linear. It is therefore
necessary to monitor this iterative procedure and modify the step length if
necessary.
The evaluation of S, and S, should be performed twice during each iteration.
Considering Fig. 9.21, this is during the calculation of the injected currents into
the network and the calculation of the nonintegrable variables. Provided the
discontinuity is not encountered, convergence is achieved in one or two iterations
at each re-evaluation especially if the saturation factors are extrapolated at
the beginning of each step.
Fig. 10.6. Generalized detailed turbine model including H.P. and interceptor valves
Some normal compound turbine configurations are shown in Fig. 10.7 and
Table 10.1 gives typical values for these configurations using the generalized
model. A hydroturbine can also be represented and the values given in Table 10.1
-- Shaft
To condenser
Control
,
Valve valves
position :eS;" Shaft
TO condenser
(b)
valve ,%
!:\:
position steam
chest
1Reheater
valve
position
w$kz!
steam HFf IP shaft
chest
1~rossoverj
I
I%--
TO condenser
L P shaft
Control
Valve valves,
position * steam
chest
I Reheater1 I Reheated ICrossover
Tandemcompound
singlereheat
Tandemcompound
doublereheat
Cross-compound
singlereheat
Cross-compound
single-reheat
Cross-compound
double-reheat
Hydro 9.9a 0 4Tw - - -2 0 3 0 0 0 0 0
are justified by the method of representing a lead-lag circuit described in
Chapter 9, with the time constant T , set at Tw in the case of the simplest model.
The full set of equations for the detailed turbine model is:
pG4 = (P,, - G4)/T4 (10.3.1)
pG5 = (G4 - G s ) / T ~ (10.3.2)
Piv = G5-Pvi (10.3.3)
pG6 = (Piv - G6)/T6 (10.3.4)
and that, in the initial steady state, the interceptor valve, if present, will be fully
open (Pvi = 1 ) in which case:
Mechanical equations
It is necessary to express the equation of motion of an induction machine in terms
of torque and not power. Also symmetry of the rotor makes its angular position
unimportant, and slip (S) usually replaces angular velocity (a)as the variable,
where:
Assuming negligible windage and friction losses and smooth mechanical shaft
power, the equation of motion is:
pS = (Tm - Te)/(2Hm) (10.4.2)
where Hm is the inertia constant measured in kWs/kVA established at
synchronous speed. The mechanical torque (Tm) and electrical torque (Te) are
assumed to be positive when the machine is motoring.
The mechanical torque Tm will normally vary with speed, the relationship
depending on the type of load.
A commonly used characteristic is:
where k = 1 for fan-type loads and k = 2 for centrifugal pumps. A more elaborate
torque/speed characteristic can be used for a composite load, i.e.:
which can include the effect of striction when start-up is being considered.
In terms of slip the torque is thus:
where
B a{b + 2c)
C a c
The values of A, B and C are determined from the initial (steady-state)loading of
the motor and hence its initial value of slip.
The electrical torque Te is related to the air gap electrical power by the
electrical frequency which is assumed constant and hence:
Electrical equations
A simplified equivalent circuit for a single-cage induction motor is shown in
Fig. 10.8, with R1 and Xl referring to the stator and R2 and X2 referring to the
Tb =
+
( X 2 Xm)
2nf OR2
and the open-circuit reactance X o is:
The reactances are unaffected by rotor position and the model is described in
the real and imaginary components used for the network, that is, in the
synchronously rotating frame of reference. Thus, for a full description of the
model, the following equations are used:
When a torque slip characteristic of the motor is available, then a simple solution
is to modify the torque-slip characteristic of the single-cage motor model.
Double-cage or deep-bar rotors have a resistance and reactance which varies
with slip. A cage factor Kg can be included which allows for the variations of
rotor resistance.
An alternative to the cage factor is the use of a better rotor model, though this is
often restricted by the unavailability of suitable data.
Induction motor loads having double-cage or deep-bar rotors can be
represented in a similar manner to a single cage motor.(10)*(")It is assumed that
the end-ring resistance and that part of the leakage flux which links the two
secondary windings, but not the primary, are neglected. The steady-state
At any instant during a transient stability study, the rotor impedance may be
assumed to be the steady-state value given above.
Analysis similar to that used in developing equations (10.4.10) to (10.4.13)
gives :
pE; = - 2nf o - S ( E-
~E 3 + pEi + (Ei - E: - (X' - X")Il,,,)/T;; (10.4.20)
PEL = 2rcf o.S(Ei - Ei') + pEk + (EL - EL + ( X ' - X")I 1,)ITs (10.4.21)
T;; =
X3 + ( X 2 .X m ) / ( X 2+ Xm)
2nf OR3
If the rotor is of the deep-bar type, then the parameters of the equivalent
double-cage type may be determined using equation (10.4.16)and (10.4.17).The
rotor parameters at zero slip are:
and
X2-Xx
X3=
(X2 - Xx)
where
---
from assuming R2 $R 3 and X 2 SX 3. --
-
the number of variables reduces to two and a simple iterative procedure yields a
result in only a few iteration^."^) A reasonable starting value is X 2 +x,derived
H 'inj,
'inj,
where the minus sign confirms the induction machine as being assumed to be
motoring.
Rectifier loads
Large rectifier loads generally consist of a number of series and/or parallel
connected bridges, each bridge being phase shifted relative to the others. With
these configurations, high-pulse numbers can be achieved resulting in minimal
distortion of the supply voltage without filtering. Rectifier loads can therefore be
modelled as a single equivalent bridge with a sinusoidal supply voltage at the
terminals but without representation of passive filters. This model is shown in
Fig. 10.10.
Rectifier loads can utilize a number of control methods. They can use diode
and thyristor elements in full or half-bridge configurations. In some cases, diode
bridges are used with tap changer and saturable reactor control. The effect of the
saturable reactors on diode conduction is identical to delay angle control of a
thyristor over a limited range of delay angles. All these different control methods
Fig. 10.10. Static rectifier load equivalent circuit
can be modelled using a controlled rectifier with suitable limits imposed on the
delay angle (a).'13'
STATIC LOADS
Operating under constant current control, the d.c. equations are
where A is the constant current controller gain and Ids is the nominal d.c. current
setting as shown in Fig. 10.11.
Natural voltage
characteristic
Constant current
control characteristic
Ids Id
Fig. 10.11. Simple rectifier control characteristic
DYNAMlC LOADS
The basic rectifer load model assumes that current on the d.c. side of the bridge
can change instantaneously. For some types of rectifer loads, this may be a valid
assumption, but the d.c. load may well have an overall time constant which is
significant with respect to the fault clearing time. In order to realistically examine
the effects which rectifiers have on the transient-stability of the system, this time
constant must be taken into account. This requires a more complex model to
account for extended overlap angles, when low commutating voltages are
associated with large d.c. currents.
When the delay angle (a) reaches a limiting value, the dynamic response of the
d.c. current (I,) is given by:
where L, represents the equivalent inductance in the load circuit. Substituting for
V, using equation (3.2.5) gives:
The slow response of the d.c. current when a large disturbance has been applied to
the a.c. system can cause the rectifier to operate in an abnormal mode.
After a fault application near the rectifier, the near normal value of d.c. current
(I,) needs to be commutated by a reduced a.c. voltage. This causes the
commutation angle (p)to increase and it is possible for it to exceed 60". This mode
of operation is beyond the validity of the equations and to model the dynamic
load effects accurately it is necessary to extend the model.
The full range of rectifier operation can be classified into four modes.(14)
Mode 1-Normal operation. Only two valves in the bridge are involved in
simultaneous commutation at any one time. This mode extends up to a
commutation angle of 60".
Mode 2-Enforced delay. Although a commutation angle greater than 60" is
desired, the forward voltage across the incoming thyristor is negative until
either the previous commutation is complete or until the firing angle exceeds
30". In this mode, p remains at 60" and a ranges up to 30".
Mode 3-Abnormal operation. In this mode, periods of three-phase short
circuit and d.c. short circuit exist when two commutations overlap. During this
period there is a controlled safe short circuit which is cleared when one of the
commutations is complete. During the short-circuit periods, four valves are
conducting. Commutation cannot commence until 30" after the voltage
crossover.
for the different modes of operation. Equations (3.2.5) and (3.2.7)do not apply
for a rectifier operating in Mode 3 and they must be replaced by:
and
a
I , = ---v,e,,(cos a' - cos y')
fixc
Fourier analysis of the waveform leads to the relationship between a.c. and d.c.
current given by equation (3.2.9) where the factor k is now:
3(/ - 2a' - / - 2y' -j2p)
k=
~ ( C Oa'S- COS y')
where
and
A graph showing the value of k for various delay angles (a) and commutation
angles (p) is shown in Fig. 10.13.
Commutation angle (degrees
20 40 60 80 100 120
I I 1 1
1.2-
11-
k 0
-1
.
0.9 -
0.8 -
The mode in which the rectifier is operating can be determined simply by use of a
current factor K I . The current factor is defined as:
Substitution in this, using the relevant equations, yields limits for the modes.
Mode 1:
K I I cos (60" - a)
and
K I _< 2 cos (a) for rectifier operation
Mode 2:
Mode 3:
2.
when a 1 30"
2
K I < -cos(u - 30") when u > 30" (10.5.13)
Jr
Mode 4:
when a I 30"
L
K I = -cos(u - 30") when u > 30" (10.5.14)
Js
This can be demonstrated in the curve of converter operation shown in Fig. 10.14.
Delay angle, a
D.C. link
Provided that it can be safely assumed that a d.c. link is operating in the Quasi
Steady State (QSS) Mode 1, the equations developed for converters in Chapter 3
can be used. That is, the converters are considered to be controllable and fast
acting so that the normal steady-state type of model can be used at each step in
the transient stability study.
The initial steady-state operating conditions of the d.c. link will have been
determined by a load flow and in this, the control type, setting and margin will
have been established.
CONSTANT CURRENT CONTROL
During the solution process at each iteration the control Mode must be
established. This can be done by assuming mode 1 (i.e. with the rectifier on C.C.
control) and by combining equations (3.3.2),(3.5.1)and (3.5.2)a d.c. current can be
determined as:
7
Assuming this current to be valid, then d.c. voltages at each end of the link can be
calculated using equations (3.2.5) and (3.3.2). The d.c. link is operating in mode 2
(i.e. with the inverter on C.C. control) if:
d' ,-1- 'dimodc,<O
The d.c. current for mode 2 operation is given by:
where Pd, is the setting at the electrical mid-point of the d.c. system, that is:
The correct value for Id-, can then be found from Table 10.3. Control mode 2 is
determined using equation (10.5.16) and in this case the following quadratic
Table 10.3. Current setting for constant power control from quadratic
equation
1,' 'A, 1,
Within Outside Id I
Outside Within 'A,
Within Within Greater of I,, and I,,
Greater Greater Idms
Greater Less 1,-
Less Greater l~...
Less Less 0
where
If the link is operating under constant power control but with a current margin
then for control mode 2
D. C. POWER MODULATION
It has been shown in the previous Section that under the constant power control
mode, the d.c. link is not responsive to a.c. system terminal conditions, i.e. the d.c.
power transfer can be controlled disregarding the actual a.c. voltage angles. Since,
generally, the stability limit of an a.c. line is lower than its thermal limit, the
former can be increased in systems involving d.c. links,'by proper utilization of the
fast converter controllability.
The d.c. power can be modulated in response to a.c. system variables to
increase system damping. Optimum performance can be achieved by controlling
the d.c. system so as to maximize the responses of the a.c. system and d.c. line
simultaneously following the variation of terminal conditions.
The dynamic performance under d.c. power modulation is best modelled in
three separate levels('5).These levels, illustrated in Fig. 10.15, are the a.c. system
controller (i), the d.c. system controller ii and the a.c.-d.c. network (iii).
!)
Fig. 10.15. A.c.-d.c. dynamic control structure. (i) A.c. system controller ; (ii) D.c.
system controller; (iii) A.c.-d.c. network
A c system
I controller
I
The a.c. system controller uses a.c. and/or d.c. system information to
derive the current and voltage modulation signals. A block diagram of the
controller and a.c.-d.c. signal conditioner is shown in Fig. 10.16.
The d.c. system controller receives the modulation signals AZ and A E and
the steady-state specifications for power Po current I, and voltage E,.
Figure 10.17 (a) illustrates the power controller model, which develops the
scheduled current setting; it is also shown that the current order
undergoes a gradual increase during restart, after a temporary blocking of
the d.c. link.
The rectifier current controller, Figure 10.17(b), includes signal limits
and rate limits, transducer time constant, band-pass filtering and a voltage
dependent current order limit (VDCOL).
The inverter current controller, Fig. 10.17(c), includes similar com-
ponents plus a communications delay and the system margin current (Im).
Current order
restart dynamics
IO2
Current rote Current
limit limit Regulator
delay
A, A Voltage rate
Fig. 10.17. D.c. system controller. (a) Power controller; (b) Rectifier current
controller; (c) Inverter current controller; (d) D.c. voltage controller
and
y- = -
S*
I Kcrm I
The static-load rectifier model does not depart greatly from an impedance
characteristic and is well behaved for low-terminal voltages, the injected current
tending to zero as the voltage approacheszero. Figure 10.18compares the current
l mpedance
Rectifier (C.C.C.)
due to a rectifier with that due to a constant impedance load. As the injected
current is never large, the iterative solution for all a.c. conditions is stable.
When the rectifier model is modified to account for the dynamic behaviour of
the d.c. load its characteristic departs widely from that of an impedance.
Immediately after a fault application, the voltage drops to a low value but the
injected current magnitude does not change significantly. Similarly on fault
clearing, the voltage recovers instantaneously to some higher value while the
current remains low.
When the load characteristic differs greatly from that of an impedance, the
sequential solution technique can exhibit convergence problem^,'^) especially
when the voltage is low. With small terminal voltages, the a.c. current magnitude
of the rectifier load is related to the d.c. current but the current phase is greatly
affected by the terminal voltage. Small voltage changes in the complex plane can
result in large variations of the voltage and current phase angles.
To avoid the convergence problems of the sequential solution, an alternative
algorithm has been developed.(13)This combines the rectifier and network
solutions into a unified process. It, however, does not affect the sequential
solution of the other components of the power system with the network.
The basis of this approach is to reduce the a.c. network, excluding the rectifier,
to an equivalent Thevenin source voltage and impedance as viewed from the
primary side of the rectifier transformer terminals. This equivalent of the system,
along with the rectifier, can be described by a set of nonlinear simultaneous
equations which can be solved by a standard Newton-Raphson algorithm. The
solution of the reduced-system yields the fundamental a.c. current at the rectifier
terminals.
To obtain the network equivalent impedance, it is only necessary to inject 1p.u.
current into the network at the rectifier terminals while all other nodal injected
currents are zero. With an injected current vector of this form, a solution of the
network equation (9.5.1) gives the driving point and transfer impedances in the
resulting voltage vector.
where
The equivalent circuit shown in Fig. 10.19 can now be applied to find the
rectifier current (I,) by using the Newton-Raphson technique.
The effect of the rectifier on the rest of the system can be determined by
superposition :
where
[m = [Val + [Z] I, (10.5.29)
If the network remains constant, vector [Z] is also constant and thus only
needs re-evaluation on the occurrence of a discontinuity.
Thus the advantages of the unified and sequential methods are combined. That
is, good convergence for a difficult element in the system is achieved while the
programming for the rest of the system remains simple and storage requirements
are kept low.
The equivalent system of Fig. 10.19contains seven variables (Vterm,1,,0, $, a, Vd
and I,). With these variables four independent equations can be formed. They are
equation (3.2.5) and:
When the delay angle reaches a specified lower limit (ahn), the control
specification, given by equation (10.5.33), changes to:
and equation (3.2.5) is no longer valid. The d.c. current (I,) is now governed by the
differential equation (10.5.5). If the trapezoidal method is being used, this
equation can be transformed into an algebraic form similar to that described in
Chapter 9. Equation (3.2.5) is replaced by:
The variables ka and kb contain information from the beginning of the
integration step only and are thus constant during the iterative procedure.
J2
kb = ( 1 - 2 kc. ka)Id(t)+ -a.ka. +
cos a(t) 2. ka. V,,,,
Vterm(t) (10.5.38)
71Tdc Rd .
Tdc Rd
where
and t represents the time at the beginning of the integration step and h is the step
length.
Commutation angle p is not explicitly included in the formulation, and since
these equations are for normal operation, the value of k in equation (3.2.11) is
close to unity and may be considered constant at each step without loss of
accuracy. On convergence, p may be calculated and a new k evaluated suitable for
the next step.
In Mode 3 operation, the value of k becomes more significant and for this
reason the number of variables is increased to six to include the commutation
angle p. The equations [ F ( X ) ]= 0 for the Newton-Raphson method in this case
are :
Js
-a.
71
V,er,.cos(O fi
- $ ) . f (p)- -a.
71
Vie,,cos a' 3X = 0
+Lid
71
(10.5.41)
element
I
without rectifier
I/ Extract lheverinvoltaqes
Thevenin impedance
vector for rectifier
I -
from [FI -
I - ' 1
- conditions
Calculate new initial
I Calculate K1 if mode 1 and mode
I ooeration
i 1 .
Y
I
I Calculate mismatches
I I
I YES
I
I NO
I Calculate Jacobianfor \I
II
I
II lter > :1 I(
I
I
I
I
I
b
I
Movement
Fig. 10.23. Dead band analogy and effect. (a) Physical ana-
logy of dead band; (b) The effect of a dead band on output
physical linkage problem as shown in Fig. 10.23(a).In this example, the input (x)
and output (y) move vertically. The diagram shows the initial steady-state
condition in which x and y are equal. The input (x) may move in either direction
by an amount Db/2 before y moves. Beyond this amount of travel, y follows x,
lagging by Db/2 as depicted in Fig. 10.23(b). The effect of a dead band can be
ignored by setting D, to zero.
Stepped output permits the modelling of SVS when discrete capacitor (or
inductor) blocks are switched in or out of the circuit. It is usual to assume that all
blocks are of equal size. During the study, the SVS operates on the step nearest to
the control setting. Iterative chattering can occur if the control-system output
(B,) is on the boundary between two steps. The most simple remedy is to prevent
a step change until B, has moved at least 0.55 B, from the mean setting of the
step.
The initial MVAR loading of the SVS should be included in the busbar loading
schedule data input. However, it is possible for an SVS to contain both
controllable and uncontrollable sections (e.g. variable reactor in parallel with
fixed capacitors or vice versa). It is the total MVAR loading of the SVS which is,
therefore, included in the busbar loading. Only the controllable part should be
specified in the SVS model input and this is removed from the busbar loading
leaving an uncontrollable MVAR load which is converted into a fixed
susceptance associated with the network.
Note that busbar load is assumed positive when flowing out of the network.
The sign is therefore opposite to that for the SVS loading.
2 0 MVAR
3 t
(a)
7 0 MVAR
( limits 3 0 and 1 0 0 MVAR
30MvARf 1480MVAR
( Limits 0 and 7 0 MVAR lt
In order to clarify this, consider an overall SVS operating in the steady state as
shown in Fig. 10.2qa). The busbar loading in this case must be specified as
- 50 MVAR and it may be varied between - 10 MVAR and - 80 MVAR
provided the voltage remains constant.
The SVS may be specified in a variety of ways, some more obvious than others,
the response of the system being identical. Three possible specifications are given
in Table 10.4. In the first example, the network static load will be 20 MVAR +
while in the second case the static load will be zero. The third example may be
represented by an overall SVS as shown in Fig. 10.2qb).
Table 10.4. Examples of MVAR loading specification
for SVS shown in Fig. 10.24(a)
where
Although not necessary for the solution process, the MVA output from the SVS
into the system is given by:
s = pi;vs
and hence
10.7. RELAYS
Relay characteristics may be applied to a transient stability program and the
effect of relay operation automatically included in system studies. This permits
checking of relay settings and gives more realistic information as to system
behaviour after a disturbance, assuming 100 percent reliability of protective
equipment. Reconstruction of the events after fault occurrences may also be
carried out.
Unit protection only responds to faults within a well-defined section of a power
system and as the faults are prespecified, the operation of unit protection schemes
can equally be specified in the switching data input. Thus, only nonunit
protection needs to be modelled and of these overcurrent, undervoltage and
distance schemes are the most common.
The integration then proceeds until the time-step nearest to tcb when circuit-
breaker opening is simulated by reducing the relevant branch admittance to a
very small value. Alternatively, the integration step length can be adjusted to
open the circuit at time t,,.
During the period between relay operation and tcbthe simulation of relay drop
off may be desired. In this case, if current falls below a prespecified percentage of
relay setting current then tcb is reset to a large value.
when
Setting multiple
(a)
Six segment
linear approximation
L I I I I
3 10 30 100
Setting multiple
(b)
Fig. 10.25. Inverse definite mini num time lag overcurrent relay (from BS142
(1966)). (a) Line xale; (b) Logarithmic scale
Initially, travel is set to zero and relay operation is assumed when D equals or
exceeds 1 p.u. If necessary, the relay operating time may be determined by linear
interpolation backwards over the last time interval.
( D , - 1.0)
t =t- .h
OP (4-4-h)
when
D, 2 1.0
and
D l - , < 1.0
and from this the circuit-breaker operating time is given by:
tCb = t o , + tdel
Many static relays have been designed which conform to mechanical
characteristics but they have also permitted different and more suitable
characteristics to be developed. These may be modelled in a similar manner.
Undervoltage relays
Apart from the fact that the relay-operating current is proportional to primary
voltage and not primary current, these relays should be modelled in the same
manner as instantaneous or fixed time-delay overcurrent relays.
Distance relays
As in practice, both busbar voltage and branch current signals are required, from
which an apparent impedance Z,/8, of the system at the relaying point can be
calculated. This is then compared with the relay characteristic to determine
operation.
A typical three-zone distance protection relay is shown in Fig. 10.26. Assuming
circular characteristics, then the settings of the relay may be identified by forward
reach Zrf/ measured in impedance (complex)coupled with backward reach R,
+
expressed as a per unit of forward reach. From this information, the centre ( p jq)
and radius (a) of each of the three circles in the impedance plane can be
established :
a = LZ
2 I/('+ Rb) (10.7.7)
p = iZr/(l - R,) cos Brf
q = i Z r/(1 - R,) sin Orf
In the example in Fig. 10.26, R, for zones 1 and 2 is zero.
where Rr is the rotor resistance. While R, and Ra will differ, the overall difference
between the negative-sequence impedance (2,)and positive-sequence impedance
representing the machine is so small as to be neglected in most cases. Further,
rotating machinery do not generate negative-sequence e.m.f.s and hence there is
no negative-sequence Norton injected current.
Thus, ignoring d.c. equipment, the overall negative-sequence network is
identical to the positive-sequence network with all injected currents set to zero.
Negative-sequence currents have a braking effect on the dynamic behaviour of
rotating machinery. For a synchronous machine, where torque and power may
be assumed to be equivalent, the mechanical breaking power (Pb) is(18):
Pb = I;(R, - Ra) (10.8.3)
which may be added directly into the mechanical equation of motion given by
equation (9.2.4).
A similar expression can be found for negative-sequence-breaking torque in an
induction machine where the speed of the rotor and the negative-sequence
currents are taken into account.
Type 4
P 0
o--=.o---O
where [if2]is a zero vector except at the point of fault. The vector [t2]contains
the negative-sequence voltages at all busbars from which the machine negative-
sequence currents are readily obtained.
If phase information is required, then the zero-sequence voltages at all busbars
also need to be determined, depending on the type of fault. This is done in an
identical manner to that used for the negative-sequence system.
10.10. REFERENCES
1. S. B. Crary, 1945. Power System Stability-Steady-State Stability, Vol. 1. Wiley &
Sons Ltd., New York.
2. D. W. Olive, 1966. 'New techniques for the calculation of dynamic response', IEEE
Transactions on Power Apparatus and Systems, PAS-85(7), 767-777.
3. T. J. Hammons, and D. J. Winning, 1971. 'Comparisons of synchronous-machine
models in the study of the transient behaviour of electrical power systems',
Proceedings I EE, 118 (lo), 1442-1458.
4. S. Beckwith, 1937. 'Approximating Potier Reactance', AIEE Transactions, 813.
5. A. H. Knable, 1956. Electrical Power Systems Engineering-Problems and Solutions.
McGraw-Hill, New York.
6. H. W. Dommell, and N. Sato, 1972. 'Fast transient stability solutions', IEEE
Transactions on Power Apparatus and Systems, PAS-91 (8), 1643-1650.
7. C. P. Arnold, 1976. 'Solutions of the multi-machine power-system stability problem.'
Ph.D. thesis, Victoria University of Manchester, U.K.
8. IEEE Committee Report, 1973. 'Dynamic models for steam and hydro turbines in
power-system studies', IEEE Transactions on Power Apparatus and Systems, PAS-92 .
(6), 1904-191 5.
9. D. S. Brereton, D. G. Lewis, and C. C. Young, 1957. 'Representation of induction
motor loads during power-system stability studies', AIEE Transactions on Power
Apparatus and Systems, PAS-76, 45 1-461.
10. H. E. Jordan, 1979. 'Synthesis of double-cage induction motor design', AIEE
Transactions on Power Apparatus and Systems, PAS-78, 691-695.
11. C. P. Arnold, and E. J. P. Pacheco, 1979. 'Modelling induction motor start-up in a
multi-machine transient stability program.' IEEE PES Summer Meeting. Vancouver,
B.C., Canada.
12. E. J. P. Pacheco, 1975. 'Induction motor starting in an electrical power-system
transient-stability programme.' M.Sc. dissertation, Victoria University of
Manchester, U.K.
13. C. P. Arnold, K. S. Turner, and J. Arrillaga, 1980. 'Modelling rectifier loads for a
multi-machine transient-stability programme', IEEE Transactions on Power
Apparatus and Systems, PAS-99 (I), 78-85.
14. D. B. Giesner, and H. Arrillaga, 1970. 'Operating modes of the three-phase bridge
converter', Int. J. Elect. Engng. Educ., 8, 373-388.
15. IEEE Working Group on Dynamic Performance and Modeling of DC Systems, 1980.
'Hierarchical Structure'.
16. K. S. Turner, 1980. 'Transient stability analysis of integrated a.c. and d.c. power
systems', Ph.D. Thesis. University of Canterbury, New Zealand.
17. CIGRE Working Group 31-01, 1977. 'Modelling of static shunt VAR systems for
system analysis', Electra, (51), 45-74.
18. E. W. Kimbark, 1956. Power System Stability: Synchronous Machines. Wiley & Sons,
New York.
Dynamic Simulation of
A.C.-D.C. Systems
11.1. INTRODUCTION
During disturbances, the interaction between supply systems and large power
converters, especially those used in h.v.d.c. transmission, generally results in large
waveform distortion and valve maloperations. Such behaviour cannot be
predicted with the quasi-steady-state models described in Chapter 3, because the
topological converter changes resulting from the disturbance can not be specified
in advance.
The simulation of transients associated with h.v.d.c. transmission is presently
carried out in physical models and simulators which provide the required
detailed modelling of the converter plant.
However, the pressing need to represent the behaviour of h.v.d.c. converters in
conventional a.c. power-system computer studies is exciting the development of
digital methods. Such methods are normally limited in their practical application
by the restricted size and accuracy of the a.c. system representation and by the
cost of running the programs. As a result, programs designed to analyse the
transient behaviour of static converters have tended to use a simplified a.c.
network representation.
This chapter describes a general dynamic model, formulated in terms of state-
space theory. It uses nodal analysis with diakoptical segregation of the plant
components undergoing frequent switching, to avoid involving the whole
network in unnecessary topological changes. The model accepts varying degrees
of a.c. system and converter plant representation to meet the requirements of any
particular study, e.g. small disturbances, faults, etc. In each case the individual
units of greater relevance are emphasised, whereas the rest of the power system is
synthesized into a simpler equivalent circuit.
Definitions
Let the network contain n nodes interconnecting 1inductive branches, r resistive
branches and c capacitive branches subject to the following restrictions('):
- one end of each capacitive branch (or subnetwork) is the common reference
point of the system;
- at least one end of.each resistive branch (or subnetwork) is the common
reference or a node also having a capacitive connection;
- at least one end of each inductive branch (or subnetwork) is the common
reference or a node having a resistive or a capacitive connection.
These restrictions are taken care of in the representation of the system
components.
The network nodes can thus be subdivided into three parts according to the
types of branches connected to them namely;
- a nodes: which have at least one capacitive connection.
- nodes: which have a least one resistive connection but no capacitive
connections.
- y nodes: which have only inductive connections.
The topological matrices Kf,, K:,, K:, which are the branch-nodes incidence
matrices (transposed) of the 1, r and c branches, respectively, are defined by their
general elements as follows:
K;, = 1, if node i is the sending end of branch p
= - 1, if node i is the receiving end of branch p
= 0, otherwise.
These matrices can be partitioned according to node type as follows:
Kfn = CKfaKfsKf,l
Kin = CK:aK:BK: ,I (11.2.1)
Kin = CK:aKBK:
where K:,, K:s and K:, are null by definition.
302
Voltage and current relations
In the absence of current sources, the following nodal equation applies
and
Ja = - KalIl - Karlr
= K,Ic
Evaluation of dependent variables
Premultiplying equation (11.2.9) by Kg, yields
+
V, = - RB(K,,Il K,,R; 'KL VJ (11.2.15)
where
R, =(K,,R,-'K~,)-~ (11.2.16)
Premultiplying equation (11.2.7) by K,, and substituting equation (1 1.2.3)
yields
where
I dl
1" "3
0 >= *>b -
and d.c. system
Ys
d2
Fig. 1 1 . 1 . Injected currents
which is not affected by topology changes, but is, in general, not sparse.
Inspection of equations (11.2.32) and (11.2.35) indicates that they are two
interdependent equations for the two vectors Vyand V,. To solve for these two
vectors would require an iterative process, but the associated computational
burden would be prohibitive. However, in a practical interconnected a.c.-d.c.
power system, the only interface between converter and nonconverter inductive
subnetworks occurs at the a.c. terminals of the convertor transformers. The usual
practice is to have a.c. harmonic filters, including a high pass, connected on the
a.c. side, which would by definition create a /?node at this busbar. Alternatively,
the connection of an a.c. line or static capacitors will define it as an a node. Thus, if
the restriction K,, = 0 is applied, the computational problem is removed without
imposing any real system representation restrictions. In the rare case of the
converter busbar being a y node, the formulation without 6 nodes would be more
economic.
By applying this network restriction, i.e. K y k null, to the above equations, a
direct solution is obtained, i.e.
v, = - L,KdkL; '( - RkIk+ Kia Va+ K:. v8) ( 1 1.2.37)
and
Vy= - L y K y l L ;'(El- RIIl- pLl.I1+ K h Va+ KfS Vg) (11.2.38)
Matrix L, in equation ( 1 1.2.37)is block-diagonal with the blocks correspond-
ing to node clusters which are separated by non-&nodes. This implies that a
switching operation in a converter will cause a change in only that part of L,
where the switching takes place.
As already explained, the primary windings of the converter transformers are
normally connected to either a or P nodes, and therefore these transformers can
be considered short circuited at their primary terminals for the purpose of setting
up L, or modifying it. For this purpose, the transformer phases can be simply
represented by their leakage reactances referred to the secondaries and connected
to the corresponding secondary winding ends.
The modification of L, to represent the various modes of operation can be
effected using a simple procedure which adds or removes single elements from L,
thereby avoiding re-formation and inversion for each mode.
When an additional link of inductance ( L )is connected between nodes i and j in
the network, L, is modified to the new matrix LL whose general element is given,
in terms of the elements of L,, by
where
+
S = Lii L j j - L i j - L j i + L
One of the nodes i and j may be external to the &network, in which case the
corresponding elements in equation ( 1 1.2.39)will be zero. Conducting valves are
represented by adding a link of zero inductance and the removal of an existing
branch is represented by adding a link of negative inductance equal in magnitude
to that of the removed branch and connected in parallel with it. In the latter case,
if the branch to be removed provides the only path from one part of the network
to the external system or to the common reference point, that part of the network
will become isolated. This condition is indicated by a vanishing denominator (S),
and it must be avoided.
Modes of operation
With reference to the three-phase bridge configuration of Fig. 11.1, the
interconnection of the converter transformer secondary windings and the d.c.
smoothing reactor, is determined by the number and positions of conducting
valves. For the converter bridge to be in a conducting state, two or more valves
must be ON in such a way that they provide a closed path for the current as
viewed from the d.c. terminals. With six valves in the bridge, conducting states
can be categorized into the following modes:
A- (Noncommutating mode): two valves on different sides and arms of the
bridge are ON.
B- (Normal commutation): three valves, one on each arm, two commutat-
ing and at least one conducting on each side of the bridge.
C- (Noncommutating arm short circuit): both valves on one of the arms are
ON, with no other valves conducting.
D- (Commutating arm short circuit): as in C with one or both valves on
another arm being ON. With one valve on the second arm conducting
(Dl), there is a single commutation process. When both valves on the
second arm conduct (D2), commutations are simultaneously taking
place on both sides of the bridge.
E- (A. C. short circuit): this is a case when the three secondary terminals of
the converter transformer are short circuited by the conducting valves.
This occurs when four or more conducting valves involve all three arms
of the bridge. With four conducting valves (El), two commutations occur
simultaneously which may be on the same side of the bridge. When five
valves conduct (E2),commutations take place on both sides and with six
valves conducting (E3), multiple commutations take place on both sides
of the bridge.
The bypass valve, used in many existing schemes, has been left out of the above
classification since its effect is implicitly included in Modes C, D and E.
Converter topology
The modes of operation stated above are represented by defining a set of three
nodes for each converter bridge nominally referred to as the top (a,), middle (6,)
and bottom (6,) as shown in Fig. 11.2. 6, and 6, correspond to dl and d,,
respectively, in Fig. 11.1 during normal operation and the presence of conducting
valves determines the connections of the transformer secondary terminals (a, b
and c) to any of the three nodes 6,, 6, and 6,.
Nominal
inverters
-
(d (e
Table 1 1 . 1
- -- --
shown in Fig. 11.6(b). The permanent branch list shown in Fig. 1 1 . q ~ )is
transformed to the working list, defined in terms of &nodes shown in Fig. 11.6(d).
Finally, the K,, incidence matrix, shown in Fig. 11.6(e), is derived from the
working branch list. The process is limited to the bridge nodes, and the external
connections of the converter network, established permanently at input, are
indicated by zero entries.
voltages becomes a matter of routine. The relations stated will apply regardless of
whether the valve is conducting or not, for, when the valve is conducting, its end
will be set to the same b-nodes and its voltage will correctly be calculated at zero.
The current relations, on the other hand, are entirely determined by two
factors: the currents in the converter branches, i.e. the smoothing reactor and the
transformer windings, and by the prevailing combination of conducting valves.
Considering the set of injected currents shown in Fig. 11.7, conducting valves
provide the paths for these currents to flow in the system and the instantaneous
values of valve currents can be obtained from these injected currents. The injected
currents themselves are obtained directly from the branch currents. J,, and J,,
are equal and numerically equivalent to the d.c. current. The phase injected
currents are obtained from the current in the transformer secondary windings,
e.g. for a star-connected secondary, J, is the negative of the current in the
secondary winding connected to node a and for a delta-connected secondary, the
injected currents are given by the difference between the currents in the windings
connected to the node.
The relationship between the valve currents in the top and bottom sides of the
bridge (I,, and I,,) on arm p, with the injected current at p, Jp, is
Jp =I v l -I v b
when the bottom valve (vb) is OFF:
Initial conditions
Conducting valves in each bridge are determined from knowledge of phase 'a'
voltage phase angle 8,, measured with respect to the steady-state load-flow
reference, and taking into consideration the inherent phase shift in the bridge
+
transformer 8, (0 or 30") and the delay angle a.
RepreSenting the a.c. voltage as a cosine function, the voltage crossover for
valve vl in the bridge, at time t = 0, will have occurred 'A' degrees earlier.
Referring to Fig. 11.8, A is given by:
In this condition valves v, and v, are conducting and valve v2 is due to be fired at
time F,,, which is 60 degrees later than the firing of v,, hence,
greater than 60 degrees then v , will have been conducting for more than one firing
period indicating that u, is already conducting. In either case, the sequence is
changed until the correct pair of conducting valves has been defined. This
approach ignores the commutation angle which, if included, would probably
complicate the issue without offering a substantial advantage since the initial
system conditions are approximate in the first place.
A refinement is to adjust the load-flow reference angle so that commutation
periods are avoided at the converter. This can be effected by choosing the firing
time of the next valve, k + 1, to be
317
where in the above analysis k = 1. Then if the steady-state commutation angle
is p.
C = 60 - p for six-pulse operation
and
C = 60 - 8, - p for twelve-pulse operation
This will ensure maximum period exists before any alteration in the converter
topology is effected by a valve commutation. Therefore
B=k.60-C
and
A=k.60+a-C
and the required value for 8, will be
&P=(k- 1 ) . 6 0 + 8 , + a - C
which can be achieved by altering the load-flow reference angle by ( e P- 8,)
Having established the pair of conducting valves in the bridge, the currents in
the transformer windings must be established on both sides. For a star bridge,
shown in Fig. 11.8(b), the secondary winding currents of each phase are set equal
to the given d.c. current if its bottom valve is conducting and to the negative of
that if the top valve is conducting, otherwise it is set to zero. Primary (a.c. side)
winding currents are set equal to the corresponding secondary currents with the
signs reversed.
Current is forced to flow in a secondary winding of a delta bridge even though
both valves corresponding to its phase are not conducting. In the case shown in
Fig. 11.8(c),the secondary winding of phase (c) carries a third of the d.c. current
with valves v, and v , both not conducting. In general, one winding carries two-
thirds of the d.c. current in one direction while the other two carry one-third each,
flowing in the opposite direction. The table shown in Fig. 11.8(d)gives the current
distribution in the delta winding for all combinations of conducting pairs of
valves. The primary winding currents are set to ( - fi) times the corresponding
secondary currents.
where :
w = dO/dt is the angular velocity
and
[ L ]= [y
]Lrs Lrr
Lsr =
Torque matrix:
where
and
The model described contains two damper circuits aligned with the direct and
quadrature axes respectively. The machine inductance matrix is thus 6 x 6
square. The fourth harmonic of inductance variation, although not normally
known, is considered to be of relevance and has been included.
The AVR and governor representation when included, should be treated as a
separate set of equations. This strategy is recommended due to the wide variation
of the characteristic eigenvalues of such subsystems allowing a consequential
saving in computation time.
A detailed analysis of the electric network without generator control is
perfectly adequate for a number of cycles after a disturbance or for disturbances
lasting only for a few cycles, like a temporary valve mal-functioning of a stitic
converter. These are carried out over a few cycles of fundamental frequency only
and the generator controls (AVR and governor) can either be left out or
approximated by constrained algebraic relations. When specially required, both
the AVR and governor representations can be added to the present model as
separate subsystems taking the applied field voltage and the rotor angular speed
as the interface variables, or in the case of thyristor-controlled excitation it is
possible and desirable to represent the dynamic behaviour of the thyristor bridge.
and
where I , and I , are taken as positive when they flow into their respective
windings.
Numerically (in p.u.) the self inductances are equal, as are the mutuals. Their
relationship is approximately
where L,, is the total leakage inductance of the transformer. The value of the self-
term in modern transformers is very high (e.g. 1 percent magnetizing current) and
the leakage reactance is very small (e.g. 5 percent) so that in p.u.
L , , = L,, = 100 p.u.
and thus
L , , = L,, = 99.975 p.u.
The numerical value of the off-diagonal terms is therefore close to that of the
diagonal terms. Since the formulation requires inversion of the transformer
inductance, numerical instabilities could be created by the use of inaccurate data.
If there are further windings on the same core the matrix equation can be
extended to include them.
Three-phase transformers can be treated in the same way with the possible
simplification resulting from ignoring interphase mutuals and reducing the
equations to three independent sets for the three phases. This is accurate when
banks of single-phase units are used and acceptable with five-limbed units.
Effects like phase-shifts, zero-sequence circulating currents, neutral earthing
etc. are automatically catered for by the terminal connections. Off-nominal taps
can be represented by suitable modification of the self and mutual elements of the
tap-side winding, as described in Chapter 2, noting that the current in the tap-
winding is identical to that in the main winding and that the terminal voltage is
the algebraic sum of the main and tap-winding voltages.
If the converter transformers are subject to off-nominal tap change, it is
necessary to recalculate the elements of the inductance matrix. When the
nonlinear relationship between the derived equivalent circuit and the magnetiz-
ing characteristics are known, the elements of the coefficient matrix can be
continuously updated to correspond with the instantaneous values of the
currents, voltages or flux linkages. In such cases, the formulation using 6, and L,
suggested in Section 11.2 increases the programming complexity and it is easier to
use only the L, formulation.
Transformer saturation
In the coupled-circuit model used to represent the transformers, the magnetizing
current may be directly formed by summing the instantaneous primary and
secondary currents, with the appropriate sign. Multiplication by the magnetizing
reactance will result in the instantaneous magnetizing flux linkage, i.e.
Circuits (a) and (b) provide representation of local loads at the converter
terminals in order to investigate their damping effect and harmonic penetration.
Static capacitors (c) may be present for compensation and also to represent a
d.c. cable.
Circuits (4,(e) and (f) are typical of h.v.d.c. harmonic filters and line dampers.
Owing to their large-time constant, the converter harmonic filters have a major
influence on the converter waveforms following a disturbance and they have to be
accurately represented in the small-step dynamic simulation.
Static series elements consist of linear, uncoupled elements and are represented
by series inductances and capacitances. The latter are restricted, in that they must
be part of a capacitive network which provides a path to the common reference
point of the system; in practice, these capacitors are connected in series with
transmission lines and the restriction is met.
Converter
Fig. 11.13. (a) nth order network equivalent-to represent the har-
monic impedances of n resonant frequencies; (b) (Modified) Thevenin
equivalents of network
Usually a locus would show a tendency for the system to be inductive for the
important lower harmonics (up to fifth). The a.c. filters will be capacitive below
the fifth harmonic and so a parallel resonant condition may exist around the third
or fourth harmonic, depending on the strength of the a.c. system. The stronger the
system the higher the resonant frequency will be in general. An equivalent circuit
which maintains a constant impedance angle over the lower frequency range and
provides realistic damping of harmonic voltages has been proposed by B~wles.'~
Giesner@)proposed an adapted version of this (Fig. 11.13(b))in which the source
impedance (which may be calculated from the short-circuit level at the converter
transformer busbars) is split so that there are two, not necessarily equal, reactive
elements with the resistive elements designed to give the required impedance
angle. The representation of a simple e.m.f. behind a series (source) impedance is
also provided for, although as explained, this will only be accurate at a single
frequency. The impedance angle varies markedly from system to system. Its
correct representation is important as its effect on damping of resonant
conditions such as overvoltages is critical. Typical values of 75 and 85 degrees for
receiving and sending end a.c. systems respectively have been ~uggested.'~) For
frequencies at or above the resonance of the high pass filter, the ax. system
impedance will be swamped. Thus the only assumption made is in the mid-range
of frequencies where the likelihood of a resonant condition in realistic strength
systems is small.
The use of an equivalent representation for an a.c. system therefore assumes
there are no network resonances yet gives an adequate treatment to possible a.c.
filter-system resonances. The possibility of associated over-voltages and wa-
veform distortion and their effect on the converter controllers is then catered for
adequately. Explicit busbar information would not be available at all busbars,
but the extremely efiicient simulation achieved in comparison with full a.c. system
representation justifies the use of such a model. However, for this equivalent to be
a realistic representation of the a.c. network when the study duration is long
and/or the disturbance severe it must include all the associated generator effects,
i.e. must be time-variant under disturbed conditions.
LdRdrI r I r I r I r I r I
Fig. 11.14. Monople equivalent for New Zealand h.d.v.c. link in d.c. line fault
studies
disturbances. Of course, the analysis of a single-pole fault on a bipole system
would require explicit representation of each pole and any mutual effects.
Numerical integration
The choice of integration method is limited to the so-called single-step methods
which are self-starting. This is because multistep methods require information
store from previous steps and they are normally started by establishing a few
points using a single-step procedure. Taking into consideration the frequency of
the occurrence of discontinuities in converter bridges, it can be seen that a
multistep method would hardly start to operate before a discontinuity occurs and
the single-step procedure called upon again to restart the process.
Two methods for the solution of the differential equations namely, the fourth-
order Runge-Kutta method and the implicit inte'gration procedure based on a
trapezoidal approximation are explained in Appendix I.
Discontinuities
The valve-switching instants must be determined as accurately as possible, so
that the relevant network changes may be carried out.
Valve-firing discontinuities are decided beforehand by the control system and
the integration step-length can be changed so that the firing coincides with the
end of a step.
Valve extinctions, on the other hand, can only be predicted at the expense of
slowing down the computation. It is preferred to detect extinctions after they
have occurred, when valve currents become negative. At this point, the time of
extinction can be assessed by linear interpolation. The time thus obtained is then
used to interpolate all the other quantities in the network. This is a first-order
approximation assuming that the variables change linearly over the step. If a
more accurate assessment is required then a backward integration to the
estimated time of discontinuity may be carried out with a further interpolation.
Voltage crossovers do not result in discontinuities but their occurrence must be
detected to provide a reference for the firing control system. This is achieved by
linear interpolation in a manner similar to the calculation of valve extinction
instants.
relevant blocks in L f .
Recalculate L,-'
Is Ly time variant?
Form L y
write output to
state information on
t
Read progmm control data
I
).
I System already established? 1
Type o f component
End Line Transformer Synchronous Filter Converter Reduced
machine system
0 1 2 3 4 5 6
I I I
Follow relevant procedure to: read control and
circuit porameters,establish new nodes and
branches, assign node types and sequences, set
initial currents and voltages, set up primitive
coefficient matrices,set up branch lists, etc.
Fig. 11.16. Flow diagram of data input and initial value establishment
new network nodes created by the procedure, together with the type
designation and cross-referencing information;
(iii) the initial values of branch current and network node voltages;
(iv) for converter bridges, the basic (noncommutating) nodal solution matrix
L, and its indexing information, and the initial valve states in each bridge.
Alternatively, network data and initial conditions, stored from a previous
(initializing) run may be input in its final expanded form.
The network equations are established by defining all the coefficient values of
the variables $,, $, and Q,. This process involves the inversion of the matrices L,,
L,, R, and C,, the formation of permanent incidence matrices, including an
auxiliary connection matrix used in the calculation of valve currents, and the
formation of Rg, L,, L,, and K,,, i.e. all the variable matrices. All dependent
variables and the derivatives of state variables are then calculated so that the
solution of the network equations may commence.
Figure 11.16 shows a general flow diagram of the process.
External disturbances
Externally applied disturbances, as opposed to those inherent in the network as a
result of valve switchings, can be identified by a set of data defining the
disturbance and indicating the time of application. Besides the application of a.c.
and d.c. faults other typical requirements for the dynamic simulation are:
(i) a.c. fault clearance;
(ii) d.c. line charge/restart/reversal;
(iii) if a d.c. fault is to occur, detection parameters need to be established.
Simultaneous application of disturbances should be possible.
Output
The main output stream dumps the calculated results onto a mass storage
medium at specified intervals, say 10-20 degrees of fundamental frequency.
However, output immediately before and after switching operations is also
necessary. A separate program is recommended to retrieve the dumped results
and produce time-response graphs of the required variables. It is also advisableto
produce graphs of variables which are not explicitly defined in the program, such
t
I Calculate valve currents 1
J.
IExtinction detected 71
YES
Interpolate all state
variables back to extinction
~nstant.Calculate
dependent variables
t
For each bridge
check for cross - overs,
:
T
If time to next firing
+
output time
v
If observation of dynamic
analysis indicates mismatch with
quasi-steady state d.c. link
model intmnsient stability
study adjust quasi -steady state
perf0;mance occordin ly New
tran~ientstability will?hsn
prov~demore accurate Thevenin
equivalent sequences
FStage i
where
and
The dynamic analysis then calculates instantaneous phase values for each
Thevenin equivalent source voltage at t, from these quantities, by
e, = Ex cos (8,+ wt)
eb = Excos (0, + wt - 120')
7
Monitor and control converter operation:
any switching ?
Modify equations,
produce output
I ? i.e. L1 = f ( f ) e.g.
Time variant ~~..'.rctances
saturation or synchronous machine models I
Calculate new values for
relevant blocks in L!
recalculate ~ 1 - I
+
Advance time by h
where At3 = Oj+, - Oj is the change in Thevenin source angle over the relevant
transient stability step. This may also be used to yield a value of a relative to the
actual a.c. voltage waveforms.
Scope
The computer program must handle a combined a.c.-h.v.d.c. system of
reasonable size and flexible configuration. Hence it must be capable of referring
to a h.v.d.c. converter not only in the special context of a rectifier-line-inverter
arrangement, but also in terms of a multi-terminal h.v.d.c. system, for example.
The number of transmission lines, transformers, etc., are limited only by the
storage allocated to the combined total of their equivalent circuit elements.
Provision should be made for different transmission line configurations, for both
a.c. and d.c. operation, and the number of segments used to represent the lines
should be variable within the limits of the total storage available.
Sparsity storage
The coefficient and incidence matrices in equations (11.10.1) to (11.10.10) contain
a large proportion of zero elements (e.g. Fig. 11.20). Only the nonzero elements
need be stored in order to save computer storage and time in performing matrix
operations. Further reduction in storage is achieved by storing only a sample of
any constant sets which are repeated, i.e. assuming linearity and balance, all
three-phase matrix elements, e.g. multiple-segment transmission line series
inductance matrix, and each phase of a bank of single-phase transformers. This
will not apply to transmission line nodal capacitances nor to resistances, the
latter (R,, R,, R,, RB) are stored as simple diagonal vectors whose indexing
information is inherently given by the relative position of the elements.
Branch No 1 2 3 4 5 6 7 8 9 10 11 12
From 1 4 2 5 3 6 19210311
To 0 7 0 7 0 7 8 1 1 8 9 8 1 0
(e)
Fig. 11.20. Storage of coefficient and incidence matrices. (a) Equivalent
circuit; (b) Full inductance matrix; (c) Compact inductance matrix;
(d) Branch list; (e) Node-branch incidence matrix
Coefficient matrices
The representative samples of nonzero matrix elements can be stored in a vector
of real constants C, and the indexing information in three integer vectors LC,
KC, LK, as in the example in Fig. 11.19(c).C contains the nonzero elements of
the full matrix, scanned row-wise, ignoring repeated sets. The vector LC gives the
relative address in C of the start of information of each row that is defined in the
full matrix. KC contains the column numbers of all nonzero elements in the full
matrix while LK gives the relative address in KC of the start of information for
each row.
The branch and node coefficient matrices (L,, L,, L,, C,) are also ordered in
such a way that they are composed of block-diagonal submatrices of low order.
Each block is a square matrix representing a set of mutually coupled branches
with few, if any, null elements. The full matrix inverses are therefore identical to
the original matrix in structure. This makes it relatively easy to calculated inverse
and matrix-by-vector products in a systematic manner.
Incidence matrices
Incidence matrix information can be stored in two integer vectors: the coefficient
vector (K) and the associated index vector (L).
Taking a node-branch incidence matrix as an example, L gives the relative
address of the start of information in K for each node. The coefficient vector K
contains the numbers of branches that are connected to each node, preceded by a
minus sign if the node happens to be the receiving-end of the branch in question.
This is effectively storing the matrix by rows, replacing a + 1 by its column
number, a - 1 by the negative of the column number, and completely ignoring
zeros. Figure 11.20(e) shows the node-branch incidence matrix for the example
given.
Implementation
It is recommended to use subroutines to perform the tasks of forming, storing and
indexing the coefficient and incidence matrices. As well as these tasks, which are
basically required at input, the solution of the network equations is best achieved
by using special routines to perform the coefficient matrix inversion, the
evaluation of matrix-by-vector products, etc.
The product of a compactly-stored incidence matrix D and a vector V can be
obtained by algebraically adding the elements of V whose relative addresses and
signs are given by the elements of K pertaining to each row in D. The product of
the transpose of D and a vector can be obtained by a similar modified procedure,
and it is therefore only necessary to store one form of incidence matrices. The real
coefficient matrix-by-vector product can be achieved by summing up the
products of the nonzero elements in each row and the elements of the vector
whose relative addresses are given by the column numbers stored in KC.
Multiple products, such as those used in equations (11.10.5) and (11.10.6),
can be formed by scanning the expressions from right to left, progressively
carrying out single matrix-by-vector multiplications until the final product is
reached.
11.14. REFERENCES
1. J. Arrillaga, H. J. Al-Khashali and J. G. Campos Barros, 1977. 'General formulation
for dynamic studies in power systems including static converters', Proc. ZEE, 124
(1 1), 1047-1052.
2. J. D. Ainsworth, 1968. 'The phase-locked oscillator-A new control system for
controlled static converters', Trans. ZEEE, PASS7 (3), 859-865.
3. J. Arrillaga, J. G. Campos Barros and H. J. Al-Khashali, 1978. 'Dynamic modelling of
single generators connected to h.v.d.c. converters', Trans. ZEEE, PAS97 (4), 1018-
1029.
4. H. W. Dommel, 1971. 'Nonlinear and time varying elements in digital simulation of
electromagnetic transients', Trans. ZE E E, PAS-90, 2561-2567.
5. H. W. Dommel, 1975. 'Transformer models in the simulation of electromagnetic
transients.' Paper 3-114, 5th Power System Computation Conference (PSCC),
Cambridge, England.
6. N. G. Hingorani and M. P. Burberry, 1970. 'Simulation of a.c. system impedance in
h.v.d.c. system studies', Trans. ZE E E, PASS9 (5/6), 820-828.
7. J. P. Bowles, 1970. 'A.c. system and transformer representation for h.v.d.c. trans-
mission studies', Trans. Z EE E, PAS-89 (7), 1603- 1609.
8. D. B. Giesner and J. Arrillaga, 1971. 'Behaviour of h.v.d.c. links under balanced a.c.
fault conditions', Proc. ZE E, 118 (3/4), 59 1-599.
9. M. D. Heffernan, K. S. Turner, J. Arrillaga and C. P. Arnold, 1981. 'Computer
modelling of a.c.-d.c. disturbances-Part I. Interactive Coordination of generator
and converter transient models', Trans. ZEEE. Winter Power Meeting, Atlanta.
10. K. S. Turner, M. D. Heffernan, C. P. Arnold and J. Arrillaga, 1981. 'Computer
modelling of a.c.4.c. disturbances- Part 11. Derivation of power frequency variables
from converter transient response', Trans, ZEEE. Winter Power Meeting, Atlanta.
A.C.-D.C. System Disturbances
where
Z is the bus impedance matrix,
1 is a vector of nodal injected currents due to a.c. generators,
and
J is a vector of nodal currents injected by the d.c. converters.
The model of the d.c. link assumes that the converters retain controllability
during the disturbances, and thus maintain the current setting. When a three-
phase fault occurs in one of the a.c. systems, represented by a current If flowing
out of the faulted bus p, the quasi-steady-state model resolves the d.c. link
equations. With the d.c. link on constant-current control (at one of the two
converter stations), the magnitude of the current injections, J f , will remain fixed,
but the respective phase angles of the injections will depend on the terminal
voltages and the d.c. link operating conditions (a,, 4).
From the matrix equation (12.2.1), the following expressions are extracted for
the nodal voltages at the converters and fault buses (Fig. 12.1).
L.-.-.-.-.-.-.-.-.-.-.-i.
Fig. 12.1. Reduced equivalent network
vi= vm+ z m m ( lfi Ji)+ Zmn(lf+ J f ) + ZmpIf (12.2.2)
Vf = Vn + Znm(li+ Ji)+ Znn(I;C+ J:) + z,,I~ (12.2.3)
i' = - zfIf = vp + zpm(I/,+ J/,) + zpn(I{+ J{) + zppIf (12.2.4)
where
Zf is the fault impedance,
If is the fault current, and
Vm,Vn,Vp are the prefault voltages at the rectifier, inverter and fault buses
respectively, in the absence of d.c. transmission.
Equations (12.2.2) to (12.2.4)describe a matrix set containing three interrelated
Thevenin equivalents. Each equation contains two terms which couple that
equivalent to the other two equivalents, viz.
Z, j(If + J!) for k~ m, n, p
where
Ji = 0 and I; = If
This equation set can be easily incorporated into the transient converter
simulation of Chapter 11 by modifying equations (11.lo. l), (11.lOS) and (11.10.8)
to allow for the mutual coupling term, which has both real and imaginary
components of impedance, i.e.
(i) R, now contains off diagonal terms, and
(ii) L, contains additional off diagonal terms.
The initial conditions are obtained using the process described in Section
11.11. The quasi-steady-state fault study, which provides the Thevenin equiva-
lents, also provides the steady-state prefault data necessary for the dynamic
analysis to begin the establishment of correct dynamic initial conditions.
The applicability of the equivalent circuit at harmonic frequencies has been
NO
simulation i.e. Fault clearance
4 YES
fault branch:
1
For each node of
fault location
Fault already cleared apply fault
on this phase
_I
so that if h' is less than the next integration step (the nominal value as calculated
in the integration routine), the next time step will be
which, assuming linear extrapolation, will coincide with the zero crossing of fault
current. Fault clearance will then be effected as described above by alteration of
an elemental impedance value to re-establish the prefault circuit.
This process is repeated for each of the phases which are faulted. A flow
diagram depicting the actions taken for fault application and removal is
illustrated in Fig. 12.3.
Fig. 12.4. (a) Full system representation of a line fault and circuit breakers; (b) Inclusion
of circuit breaker resistances
(iii) After a specified time, the breakers are reclosed and the fault to ground is
removed by reforming the original connection matrices.
Each of the breaker operations can simulate the inclusion of resistance in the
makebreak sequences. A series resistive element may be included at the relevant
line end (r,, r, in Fig. 12.4(b)),the value of which may be practically zero without
affecting convergence patterns, since the line ends are generally connected to
other lines, and are therefore in general a nodes.
Opening of the circuit breaker thus involves
(i) increasing r, and r, at a rate defined by the breaker contact resistance us
time behaviour ;
(ii) at a certain resistance, setting elements a,r, and a,+,r, to zero in the
connection matrix K,,.
Reclosing the circuit breaker involves:
(i) removing Zf by setting element a, Zfto zero in K,,;
,
(ii) reinstating elements a,r, and a, + r, to their original values (1, - 1)in K,,;
and
(iii) decreasing r,, r, to comply with the preinsertion resistor values.
In cases where breaker contact resistance is not modelled, only step (ii) in the
opening action is simulated. Similarly in cases where no preinsertion resistance is
used, only steps (i) and (ii) are performed, with r, and r, equal to zero.
-125.0
0.0 1.0 2.0 3.0 4.0 5.0 6.0 10 80 9.0 10.0 11.0 120
Time (cycles)
(a)
25
Rectifier
0.0
Inverter
0.01 I -25.01 , , , , , 1
0.0 2.0 4.0 6.0 8.0 10.0 12.0 0.0 2.0 4.0 6.0 80 10.0 12.0
Time [cycles ) Time ( cycles )
(c)
-1.5
0.0 1.0 2.0 3.0 4.0 50 6.0 7.0 8.0 9.0 10.0 11.0 12.0
Time (cycles
(d)
Fig. 12.5. Typical results for a three-phase short circuit at the inverter end. (a) A.c. fault currents; (b)
Thevenin voltages; (c) Thevenin angles; (d) Inverter a.c. voltages
The behaviour of the d.c. link is indicated in Fig. 12.6, which shows the d.c.
powers at the rectifier and inverter ends.
A second example, involving a less severe disturbance, is the case of a three-
phase fault sufficiently distant from the rectifier end to maintain about 65 percent
of the nominal voltage at the rectifier terminal. Results for the a.c. current and
0.8
0.0
-0.7
-1.5
0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0 11.0 12.0
Time (cycles
Fig. 12.6. Converter powers for a three-phase short circuit at the inverter end
-250.01
0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0 11 0 12.0
Time (cycles
(a)
Time ( cycles)
Fig. 12.7. Typical a.c. results for a threephase short circuit at the rectifier end. (a) A.c. fault
currents; (b) Rectifier a.c. voltages
0.01
0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0 11.0 12.0
Time (cycles)
Fig. 12.8. Converter powers for a three-phase short circuit at the rectifier end
voltage waveforms are illustrated in Figs. 12.7. The variation of d.c. powers is
plotted in Fig. 12.8.
2 .o
Fig. 12.9. (a) Effect of instant of fault (and a,) on Id:; (b) Effect of
detection delay on Id:
retard action of the rectifier firing angle can be initiated. Again, instantaneous
The cumulative
detection is assumed, this time with a constant firing angle of amin.
effects of initial firing angle and elapsed time before protective action becomes
effective can be obtained from these curves.
The effect of noninstantaneous detection on the level of overcurrent is shown in
Fig. 12.9(b) for a zero resistance earth fault on the line side of the rectifier end
terminal. The nominal firing angle was 25 degrees with the fault applied 10
degrees after the last firing. Calculations based on the above idealized formulae
(12.6.1) and (12.6.2) indicate levels much higher than those obtained by
the dynamic simulation, for both instantaneous detection and for detection 30
degrees after fault application.
This indicates the need for the more accurate simulation in which the a.c.
system, filters, d.c. system parameters and converter controls are accurately
modelled, in order to assess the effect of voltage regulation, commutation, and
waveform distortion. Unless the simulation can accurately represent the actual
timing ofvalve firings in relation to fault occurrence,and the subsequent effects of
rapid control action associated with normal or protective action, inaccurate
results will be obtained.
The interactive use of the state variable dynamic analysis and a transient
stability program described in Chapter 11 provides accurate simulation of the
a.c.-d.c. system behaviour over the entire fault and recovery period.
The dynamic analysis requires suitable representation of practical d.c. line fault
detection and recovery schemes. Similarly, accurate representation of the
expected fault behaviour, including the arc characteristics, must be incorporated
into the model.
which will be directly related to the travelling waves initiated by the fault, and
contains information relating to the fault type, size and location. Since the
terminations of the d.c. line will affect the gradients of the travelling waves, this
must be considered in setting suitable detection levels.
By suitable choice of K, and K,, restraints may be put on q to ensure that a
monopolar fault is not seen as a fault on the healthy pole of a bipole system, due to
the voltage which may be induced thereon by mutual coupling action. Of course,
if this induced voltage is sufficient to initiate a flashover on the initially unfaulted
pole, an extremely undesirable condition, then q must be such as to initiate the
detection/protection scheme associated with that pole as well. On detection of a
line fault, the converter bridges associated with the faulted pole(s) will then
undergo the relevant protective action as described below.
Converter restart
On completion of the deionization period, the program initiates the restart
procedure for the link. In most cases, the line is restored to normal voltage and
prefault power. If re-energization at full voltage is not acceptable (e.g. due to wet
or dirty insulators), then a lower voltage may be used by bypassing one or more of
the converter bridges, or by the use of constrained control action, i.e. larger
minimum limits on a, and y,.
Transmission cannot restart itself, since operation at the completion of
deionization is with both converters in inversion, i.e. a, = 125" and y, =yo. A
starting order is therefore needed to remove the protective control action and to
release the normal control systems of the converters. The actual time required for
the d.c. link to regain nominal voltage and current will depend on the properties
of the d.c. line and the converter controls. Generally, fast restart with a speed of
power recovery conducive to stable operation, but with avoidance of overswings
in both a.c. voltage and current, is required. Oscillations of these parameters may
otherwise occur at the natural frequency of the line, and may be sufficient to
cause arc restrike or to put additional stress on other plant. For this reason,
restart by an instantaneous release of the normal controllers is not practical.
To avoid overswings in d.c. voltage (and current), it is recommended to
increase the direct voltage exponentially to the nominal value, with minimum
overshoot, rather than as a sudden step. This unit usually operates only for start-
up or re-energization following a line fault. The time constant of the exponential
increase must be long compared to the natural frequency of the line (e.g. 100 ms,
cf. 20 ms), to ensure stable, controlled recovery of the d.c. system from the fault
condition is obtained. A similar method is used in practice on the New Zealand
scheme with an equivalent time constant near 100 ms.
If the line fault persists, or the arc restrikes due to inadequate deionization, the
above process may be repeated a few times with increased deionization periods
and/or reduced restart voltages, and then, if it still persists, the faulted pole
removed from service.
12.8. RESULTS
Fault application
Once satisfactory initial conditions for the test system are achieved (requiring
about two cycles of dynamic simulation) a pole-to-pole short circuit adjacent to
the line side of the rectifier smoothing reactor is applied, immediately after a valve
firing in one of the rectifier bridges. Thus, for twelve pulse operation the fault
occurs 30 degrees before any retard action becomes effective. Initiation of the
fault is implemented by alteration of the relevant connection matrix element (K,,
or K,,) for an already established branch element, or by insertion of a new branch
element (R,) and reformation of K,,.
Fault detection
In the state-variable formulation, the detection principle is implemented by
monitoring the rates of change of the state variables I//, and Q, on each pole of the
d.c. line at each converter terminal, where 1 is the smoothing reactor branch and a
is the node relating to the converter's line side termination. Equation (12.6.1) may
then be rewritten as
An indication of the relative detection times for the two converter units is
shown in Figs. 12.10(a) (converter direct currents) and (b) (line side converter
voltages). An indication of the travelling wave initiated by the fault can be
obtained from this graph, as can the natural frequency of the d.c. line. The rectifier
detection unit is seen to respond almost instantaneously to the fault, as the
voltage gradient is very high. The protective action initiated, with the next firing
pulse retarded, is seen to be effective after 30 degrees, as predicted. This is
indicated by a change in the rate of rise of rectifier current. The time delay
associated with the travelling waves on the line is indicated by the inverter
detection delay of around 5 ms between fault initiation and identification.
Inverter current is seen to extinguish after approximately 8 ms, and the rectifier
current 6 ms later. The oscillation of inverter line voltage is seen to be of about 50
Time (cycles)
(a)
-0.5L 1 1 I I I I I
2.0 2.1 2.2 2.3 2.4 2.5 2.6 2.7 i
Time (cycles
44 C
$.$ g .$E%
-s'esg'=
L0
s'=
zag$
- 0:
""'0
r
w c, (b)
Hz frequency. The envelope of the fault voltage is seen to have a peak of around
0.15 p.u. or about 30 kV, which will promote a healthy arc.
2 5 10 20
It is apparent that each halving of the earth resistance increases arc extinction
by approximately a half cycle. The arc extinctions occur at the zero crossings of
arc current, provided the peak arc voltage is below the arc sustain level
(approximately 2 kvlmetre). The arc path resistance then increases effectively to
infinity, since no current flows, and provided restriking does not occur, the fault
branch may be removed or disconnected. This action is invoked at the
completion of the specified deionization time, when it is assumed the arc path is
completely deionized.
00.0
0.0 1.0 2.O 3.0 4.0 5.0 €
Time (cycles)
where
a, is the control angle which will give nominal line voltage, and
At is the elapsed time since transfer of the converter control to the restart
unit (i.e. from when a, = 90").
The constant k will determine the response rate.
Figure 12.12 illustrates the variation of a, for various values of k, as a function
of elapsed time (At).In this case a, = 10". For a d.c. line which has a time constant
of approximately 20-25ms, it is desirable to have no more than 63 percent of
nominal voltage when At = 25 ms, in order to obtain critical damping of the no-
load voltage rise. Curve (ii) exhibits this characteristic, and actually achieves
nominal voltage when At is approximately 100 ms.
The recharge process is thus implemented by adjusting the firing instant F j + l
by
where
Detection
I Fault present
1
V NO
4
{ ~ i n erecharg! in pragr
YES
1
.t-
I Set a, = 125' 1
or Y l = YO
render
controllers
inoperative
1
NO
{
1 'I
Extinction
Deionization
Fault removal
i
1 Remove fault
JYES
1
1
NO
r
-
YES
a, > 90'
t
NO Set a,= 90'
D.C. line 1 1
Line rechar ed -
restore CC! regulators
to both rectifier ond
mverter
160 18.0 2
Time (cycles)
(b)
0
Time (cycles)
(c)
Time (cycles)
(d)
-1.51
0.0 2.0 4.0 6.0 8.0 10.0 120 14.0 16.0 180 2
Time (cycles)
(el
25.0
12.5
0.0
-12.5
-25.0
Time (cycles)
(f)
12.5
0.0
-12.5
- 25.0
0.0 2.0 4.0 6.0 8.0 10.0 12.0 14.0 16.0 18.0 20.0
Time (cycles)
(g)
Fig. 12.14. (a) Converter d.c.currents ;(b) Converter d.c. voltages ;(c) Converterd.c. powers ;
(d) Rectifier a.c. voltages; (e) Inverter a.c. voltages; (0 Rectifier a.c. currents;(g) Inverter a.c.
currents
The fault was applied and detected using the process described in Section 12.8.
As a, is retarded to 125", Fig. 12.14(f) shows that two phases at the rectifier end
continue to conduct current for an extended period. The initial rectifier current
peak, which occurs 4 ms after the fault is applied, is only 80 percent above the
setting and decays to zero, 15 ms after fault initiation. The inverter has already
shut down by this time. Figure 12.14(b) illustrates the existence of voltage
oscillations, with the inverter voltage reaching 50 percent of the rated voltage
with reverse polarity. The oscillations die out within approximately 80 ms. Due to
the continued presence of the harmonic filters, a residual current level remains
throughout the fault (Figs. 12.14(f)and (g)), with considerable initial harmonic
content causing corresponding a.c. voltage waveform distortion (Figs. 12.14(d)
and (e)).This is more evident at the rectifier end due to the large initial increase in
direct current and also due to the large change in control angle, which does not
occur at the inverter end under E.A. control. The voltage distortion is seen to
reduce quickly as the residual (filter) currents become settled. In each case this
takes about five cycles, giving an indication of the relatively slow time response of
the harmonic filters.
During the same period the line voltages are seen to decay, and as indicated in
Fig. 12.14(b), extinction of the arc occurs some 80 ms after fault initiation.
During the deionization period, which extends for approximately 100 ms, the
a.c. voltage waveforms become sinusoidal but some amplitude variation is
noticeable in Figs. 12.14(d) and (e) as a result of the time-variant Thevenin
sources.
Some 200 ms after fault initiation, the arc is expected to be completely
deionized, the fault branch is removed and fault recovery commences. The
rectifier direct current and both line voltages indicate the effect of the exponential
recharge function. The current involved in line recharge is seen to be relatively
small. Because the inverter is o n E.A. control, the inverter does not conduct
current until the line is charged to about 90 percent of nominal. The recharge
time, measured from deionization until the line is recharged and converter
current has begun to increase, is seen to be about 80 ms.
Following the recharge period, which is completed 280 ms after fault initiation,
the converter currents are seen to increase quickly and steadily to the set values.
This increase takes place with the rectifier on C.C. control, and the inverter on
E.A. control. The current is seen to reach its setting after 60 ms with only a minor
initial overshoot, indicating that the C.C. control constants used are adequate. At
this time (cycle 19), the d.c. power levels are not a t their prefault values since the
inverter a.c. voltage is still depressed, although it is increasing steadily with the
reinjection of link power.
12.9. REFERENCES
1. J. Arrillaga, M. D. Heffernan, C. B. Lake and C. P. Arnold, 1980. 'Fault studies in a.c.
systems interconnected by h.v.d.c. links', Proc. IEE, 127, Pt. C (l), 15-19.
2. H. A. Peterson, A. G. Phadke and D. K. Reitan, 1969. 'Transients in e.h.v.d.c. power
systems: Part I-rectifier fault currents', Trans. IEEE, PAS-88 (7), 98 1-989.
3. M. D. Heffernan, J. Arrillaga, K. S. Turner and C. P. Arnold, 1980. 'Recovery from
temporary h.v.d.c. link faults', Pans. IEEE, Paper 80SM 67.5-679. PES Summer
Meeting, Minneapolis.
4. N. G. Hingorani, 1970. 'Transient overvoltages on a bipolar h.v.d.c. overhead line
caused by d.c. line faults', Trans. IEEE, PAS-89 (4), 592-610.
5. E. W. Kimbark, 1970. 'Transient overvoltages caused by a monopolar ground fault on
bipolar d.c. line: theory and simulation', Trans. IEEE, PAS-89 (4), 584-592.
6. A. Erinmez, 1971. 'Direct digital control of h.v.d.c. links.' Ph.D. thesis, University of
Manchester.
7. A. Kohler, 1967. 'Earth fault clearing on an h.v.d.c. transmission line with special
consideration of the properties of the d.c. arc in free air', Trans. IEEE, PAS-86 (3), 298-
304.
8. E. Uhlmann, 1960. 'Clearing of earth faults on h.v.d.c. overhead lines', Direct Current, 5
(2), 45-57, 65-66.
Transient Stability Analysis of
A.C.-D.C. Systems
13.1. INTRODUCTION
Having developed a detailed transient converter simulation program (TCS)
suitable for analysing the behaviour ofan h.v.d.c. link just after a disturbance, it is
only necessary to merge this model with a transient stability (TS) program to
complete the suite of programs. It is then possible to examine the response of all
system components over a period of several seconds.
Many of the problems associated with the merging of the two programs have
been discussed in Chapter 11. This was done to ensure that acceptable a.c. system
information was made available to the TCS programme. This chapter firstly deals
with the processing of the information, as obtained from the TCS, into a form
suitable for inclusion in the TS programme.
The second part of this chapter is concerned with three transient stability
analysis studies of the New Zealand a.c.-d.c. power systems. These studies are
used to illustrate principles as well as to justify the use of the integrated programs.
Time (cycles)
(a)
-20 0
:1 5 13.5 145 15.5 165 115 185 19.5 205 215
I
22.5
Time (cycles)
(b)
Fig. 13.1. Comparison between voltage and current waveforms at the rectifier terminal
during current setting changes (from a TCS study). (a) Rectifier voltages; (b) Rectifier
currents
The spectral analysis of the voltage waveform is, therefore, less prone to error
and the voltage magnitude is a better choice of variable.
A specification of voltage and real power is required at each integration step of
the TS study when the quasi-steady-state d.c. link model is replaced by the results
of TCS.
To obtain fundamental components of voltage and power, the TCS waveforms
must be subjected to Fourier processing over a discrete number of fundamental
cycles of the waveforms. A component of voltage and power can be obtained at
any sample point using a complete cycle (an aperture size of period T) spanning
T/2 both sides of the sample point as illustrated in Fig. 13.2(a). The sampling
points must be synchronized with the TS study integration steps and, in order to
ensure that sufficient samples are usedJ4)the TS program step length should not
exceed half of one fundamental cycle. It is convenient, therefore, to fix the step
length and sample rate at half a cycle of the fundamental frequency as illustrated
in Fig. 13.2(b).
Cycle of ~ n t e r e s t ~
t it--
: , TCS waveform
\
(a)
Seauence of
I
I
I
I
Aperture for
somple points
(b)
-1 21
00 06 12 18 24 30 36 42 48 54
Time (cycles)
(b)
Fig. 13.3. Typical waveforms obtained from a TCS study. (a) Inverter voltage
waveform for a three-phase fault at the inverter Terminal; (b) Rectifier voltage
waveforms for a single-phase fault in the inverter a.c. system
(ii) The nonlinear behaviour of the converter excites disturbances which are
not periodic in T.
(iii) The network equivalent source is given in a constant frequency frame of
reference but the TS study deviates from constant frequency after a
disturbance, effectively altering the fundamental period.
Spectral analysis
The above effects can be classified in terms of modulation, frequency mismatch
and nonperiodic noise. Their presence means that a simple fundamental
component of voltage (or current) cannot be accurately obtained merely by using
a standard discrete Fast Fourier Transform (FFT),(5)on the samples which make
up the TCS waveforms of Fig. 13.3.
-1.5 u
2.0 2.2 2 2.6 2.8 3.0
Time (cycles) Time (cycles)
(a) (b)
Fig. 13.4. Single cycles of one-phase obtained from TCS waveforms of Fig. 13.3 (a)
Inverter voltage; (b) Rectifier voltage
-8.0 I I I I I I I I
0.0 1.2 2.4 3.6 4.8 6.0 7.2 8.4 9.6 10.8 1:
Time (cycles)
-2.0 -
0.0 0.6 1.2 1.8 2.4 3.0 3.6 4.2 4.8 5.4 r
Time (cycles)
(b)
Fig. 13.5. Comparison of RMS and fundamental powers for the cases
shown in Fig. 13.3. a) Three-phase fault at the inverter end b) Single-
phase fault at the inverter end
harmonic currents and voltages. Considering the relatively low resistive com-
ponent of the converter transformer and the large X I R ratio of the a.c. system, the
in-phase component of the harmonic currents and voltages are very small.
Therefore the use of rms power is a good approximation and spectral analysis is
not needed to derive the power variable.
The rms value of a voltage, V(t), is defined over a period T as
Time (cycles)
(a
Time (cycles)
(b)
Fig. 13.6. Comparison of RMS and fundamental voltages for the TCS
waveforms of Fig. 13.3. (a) Inverter voltage; (b) Rectifier voltage
severely distorted case with a 5 percent difference observed during the period
immediately after disturbance initiation. This coincides with the period of
maximum distortion in Figure 13.3(b) at which time harmonic content is high.
This error is too large to be considered insignificant when transferring data from
the TCS to the TS programme.
It is therefore necessary to retain the spectral analysis to determine a
fundamental voltage magnitude from the transient converter simulation wa-
veforms but the rms approximation for power avoids the need for spectral
analysis of current.
The information from the TCS is embodied in real power (P,,,), and voltage,
(Ed,& variables obtained from processing of the TCS waveforms. These become
the specifications at the link terminal node when the QSS model is removed.
Using the unified algorithm developed in Section 10.5 the system equivalent
shown in Fig. 13.7 is solved, given the P,,, and Ed,, specifications.
This system can be used to solve directly for the a.c. currents which can then be
superimposed on the a.c. network in the same way as the a.c. current from the
QSS model.
The complex power Sf& is given by
Taking real parts of equation ( 13.4.2)and rearranging, 8 can be obtained directly
by:
~=~-P+COS-
The a.c. current can now be obtained using
This solution requires little computational effort when compared to the QSS
link model. With very low converter terminal voltages, the first term of equation
(13.4.3) becomes very large and any mismatch between Ed,, and P d y ncan cause
the cosine function to exceed unity. This difficulty is overcome by releasing the
power specification and allowing it to change until the equation can be solved.
The small adjustment of power setting has little effect on the transient stability
analysis, since the converter power involved is small under very low voltage
operating conditions.
Alternatively, the fault can be separated from the admittance matrix to obtain a
solution for the fault current using the Thevenin equivalent at the faulted node.
The fault current can then be superimposed onto the network in the same way as
the injected current obtained at the d.c. link nodes.
It is necessary to solve both the fault and the d.c. link simultaneously because of
their mutual coupling'shown in equations (12.2.2) to (12.2.4).
The simultaneous solution can be achieved by reformulating the d.c. link
.1
due to d.c. link obtained
Estimateof Theveninvdtoges
For this to be a correct solution, the Thevenin theorem states that all other
voltage and current sources must be present, i.e., the network equivalent at the
fault bus must be obtained with the d.c. link injected currents included. The
converse also applies in obtaining the correct solution for the d.c. link. The
solutions of both link and fault are combined in the algorithm of Fig. 13.10.
The algorithm relies on a good estimate of the d.c. injected currents at the first
iteration of each step if the speed of convergence of the TS program is not to be
impeded. However, the link solution obtained at the previous step is suitable
except for the first iteration after fault application. In this case, since the fault has
a dominant effect on the network, it is solved alone at the first iteration. A new d.c.
link solution can then be obtained for entry into the algorithm (at A) of Fig. 13.10.
This starting process is illustrated in Fig. 13.11.
The provision of a time variant system from the TS program is much simpler
with a d.c. fault. For this case, it is not necessary to provide explicit fault bus
representation since the fault is contained within the d.c. system being simulated.
The equivalent provided for the TCS program in this case is exactly the same as
that obtained for the QSS d.c. link model.
t
Solve for fault current 1
using [v,] [ 11
-
Filte
In this representation, the transient stability program, using a QSS d.c. link
model, provides the time variant system equivalents for the TCS. In turn the
instantaneous values of I,, I , , E l and E , of Fig. 13.12 are obtained from the TCS
to provide r.m.s. power and fundamental frequency voltage at the converter
terminal bus, for use in the transient stability programme.
Three worst case studies consisting of a rectifier a.c. system (South Island) fault,
an inverter a.c. system (North Island) fault and a d.c. system fault are e~aluated!~)
( seconds )
h v ( 1 1QSS
-5 0 Foult off ( 2 ) TCS fixed system equ~valent
( 3 ) TCS tlme var~ontsystem equ~volent
Fig. 13.13. Voltage and reactive power at the rectifier for the
rectifier a.c. system fault. (a) Terminal ax. voltage;
(b) Reactive power
Three sets of results are shown, the first utilizes the QSS model alone without any
use of TCS. The second uses the TCS together with the TS simulation but there is
no coordination between the two parts of the program and thus fixed a.c. system
equivalents are used in the TCS. The third set of results is obtained when the two
simulations are coordinated and the TCS has a time variant a.c. system
equivalent.
For the second case, the power factor of the injected current is forced to change
from lagging to leading as a result of specifying voltage and real power at the
converter node. This results in the impossible condition of the converter
supplying reactive power to the system. The unsuitability of this form of
modelling is further emphasized during the postfault period where the TCS
recovery voltage considerably exceeds that of the dynamically responsive system.
The third set of results using interactive coordination gives a much better
match between the QSS model performance and the TCS results. The departure
between the two models immediately after fault inception and fault removal is
due mainly to the different filter models. In the QSS model, filters are represented
by a shunt admittance while the TCS reflects the full dynamic response of the
filters. For the TCS representation, the combination of the dynamic filter
response and control action results in an initial drop in reactive power demand
from the a.c. system. The mismatch is more pronounced following fault clearance.
A further cause for mismatch is the requirement to align the TCS results with
the TS integration steps, at the disturbance times, as shown in Fig. 13.8. This
distorts the behaviour of the TCS results immediately after fault removal.
Although the reactive power profile is useful in indicating the accuracy of
matching between the two programes, it is the real power transferred by the d.c.
link which affects transient stability analysis. Fig. 13.14 shows the real power flow
at the link terminals. The inverter power (Fig. 13.1qb)) is shown with positive
power flow out of the link terminal to facilitate a comparison between the two
ends of the link.
5.0.
t
4 5-
4.0 -
-
a
3.5-
- ( 1 ) QSS
3 3.0 ( 2 ) TCS fixed system equivalent
( 3 ) TCS time variant system
~'25- equivalent
eauivalenl
--
I I
I
-5.0
Inverter I
Time (cycles
Fig. 13.15. Converter d.c. currents for the rectifier a.c. system fault
In the postfault period in particular, the fixed equivalent simulation (curve 2),
departs considerably from the time variant case (curve 3) and may be omitted
from further discussion.
Immediately after fault inception, the TCS indicates that the rectifier d.c. power
flow is less than the level predicted by the QSS model. This is caused by the
sudden drop in the rectifier d.c. voltage, giving rise to a drop in d.c. current as
illustrated in Fig. 13.15. Control response at the inverter restores the d.c. current
and the link recovers from the initial disturbance quickly.
At the inverter end the immediate power response predicted by the QSS model
deviates considerably from the TCS solution, which reflects the delay effect
caused by the d.c. line.
Although instantaneous power differences are observed, the nett energy
transfer as calculated by the TCS and QSS models are similar, thus causing
relatively small first swing maximum angle differences as shown in Table 13.1.
For the high-inertia machines at the rectifier terminals (Benmore) the differenceis
3 percent, while for the low-inertia synchronous condensers at the inverter
(Haywards), the difference is 12 percent. At both terminals the QSS model gives
the larger maximum angles, and thus provides pessimistic but safe results.
The solutions of the TCS and QSS models are not exactly coincident at the
beginning and end of the TCS study as shown in Table 13.2. The results are
Table 13.1. Maximum swing angles of synchronous machines at the d.c. link terminals for the
rectifier a.c. system fault
Benmore Haywards
Case (rectifier) (inverter)
(1) QSS
(2) TCS fixed equivalent
(3) TCS variant equivalent
Table 13.2. Coincidence of solutions at start and end of TCS for the rectifier a.c. system fault
obtained at the last integration step of the TS program for which a TCS data
point is available. A QSS solution of the link is obtained at this step but the
injected a.c. currents are derived from the final TCS data point (refer to Fig. 13.8).
In this case, the mismatches are small and maximum mismatch occurs for the
inverter power and is approximately 3 percent of the QSS value.
-8.0I I I I I I I I I I
0.0 1.2 2.4 3.6 4.8 6.0 7.2 8.4 9.6 10.8 12.0
Time ( cycles
Fig. 13.16. A.c. and d.c. apparent power at inverter terminal for the
inverter a.c. system fault
0.0 1.2 2.4 3.6 4.8 6.0 7.2 8.4 9.6 10.8 12.0
Time (cycles 1
Fig. 13.17. A.c. and d.c. apparent power at rectifier terminal for the
inverter a.c. system fault
1 1 I I I I ,
0.1 0.2 0.3 0.4 0.5
Time (seconds)
(a
A
4.0 - Start TCS End TCS
3.0- ( 1 QSS
--z
U
2.0-
4I!
( 2 ) TCS
0- 1.0-
ii
0 I I I ,
0.3 0.4 0.5
-1.0' Time (seconds)
Fig. 13.18. Voltage and reactive power at the inverter for the inverter a.c.
system fault. (a) Terminal voltage; (b) Reactive power flow
the real power obtained in deriving the TCS equivalents, includes .the real power
flowing into the fault. The difference between the a.c. power flowing into Bus 2,
and the d.c. power flow at the inverter is illustrated in Fig. 13.16.
Immediately after the fault is applied, the inverter experiences repeated
commutation failures which persist during the fault period. The commutation
failures produce oscillations in the d.c. power at the inverter terminal, as shown in
Fig. 13.16, illustrating the effect which they have on the d.c. link performance. A
comparison between Fig. 13.16 and 13.17 shows that thed.c. power at the rectifier
terminal is smoothed by the d.c. line response. The oscillations in Fig. 13.16 are
therefore not due to the line.
Figure 13.18 illustrates the difference between the QSS and TCS a.c. voltage
and reactive power at the inverter terminals. This case has a similar filter response
to that for the rectifier a.c. fault study. In the post fault period the reactive power
demands of the link are lower in the TCS results because of the slower link
recovery for this model.
A good match is obtained for the voltages, during the fault period becaclse of
the dominance of the fault in the TCS results. In the postfault period the match is
not so accurate because of the slower link response determined by transient
converter simulation.
The active power flows at the two terminals of the link are shown in Fig. 13.19.
The QSS model shows significant departures from the results obtained using
TCS. The differences between the results are caused by three factors:
(a) Due to the d.c. line response, power continues to be fed into the rectifier end
of the line for one cycle after the fault occurs at the inverter bus.
I .
i!
-3.0
- li ( b )
il
Fig. 13.19. Real power flows at the d.c. link terminals for
the inverter a.c. system fault. (a) Inverter; (b) Rectifier
At the inverter terminal, immediately before fault clearance takes place (i.e.
the first detected current zero after breaker operation), a commutation
failure occurs. This converter disturbance delays the clearance time of the
fault. It is impossible for a QSS model to detect this sequence of events.
In addition to the delay due to the commutation failure, inverter recovery is
also affected by the response of the d.c. line and the controller at the rectifier
end. As shown in Fig. 13.19(b)the TCS rectifier response is slower than the
QSS response.
The accurate response obtained from the transient converter simulation affects
the results of the TS analysis to a greater degree than for the rectifier a.c. system
fault study, as shown by the swing curves in Fig. 13.20. In this case, the maximum
swing angle obtained using the QSS model is less than that obtained using the
TCS model (i.e. the results are optimistic) and this is especially noticeable for the
low-inertia machines at the inverter terminal.
-901 I I I I I
0.1 0.2 0.3 0.4 0 . 5 ~
Time (seconds)
Table 13.3. Coincidence of solutions at start and end of TCS for the inverter a.c. system fault
4.0
0.0
-4.0
-
-
-
-
D.C.-
L
-80 I I I I I I I I I
0.0 20 4.0 6.0 8.0 10.0 12.0 14.0 16.0 180 20.0
Time (cycles)
Fig. 13.21 . A.c. and d.c. apparent power at rectifier terminal for d.c. line fault
1
Start TCS
5.0
I End TCS
-
E
-3
, 3.0
4.0
I
O?' 0.2 0.3 0.4 0.5 0.6
Time (seconds)
Fig. 13.22. Power flow to rectifier terminal for the d.c. line fault
However, as the TCS progresses, it becomes clear that transmission will not be
re-established at this time and matching between the two programs will be poor.
A number of subsequent equivalents may be obtained. The one adopted in this
case has the restart delayed and is ramped to nominal d.c. current linearly over a
two-cycle period. This equivalent provides excellent matching between the two
models as illustrated by curve 3 in Fig. 13.22.
As shown in Fig. 13.22, there is a 7 percent overshoot in the TCS link power at
full recovery. This is due to the finite rate of response of the controller. The
simulation is continued after recovery for a further four cycles to ensure that the
oscillation is damped and converges to the same conditions as the QSS model.
Throughout the d.c. fault there are no consequential converter disturbances
and because of this, and the accurate matching achieved with the QSS model, the
differences between the QSS and the TCS models are small.
Start TCS End TCS
.-. 'V
(1 QSS Matched
0.9 ( 2 ) TCS
I I I I I I
0.8
0.1 0.2 0.3
Time (seconds)
0.4 0.5 0.6 *
(0)
A
Start TCS
-
2.0 -
1.0-
-. -. '1.\
1
!
( 1 QSS Matched
(2)TCS
n
a
67 I
0 I
I 0.2
I
1.0 - I
.A-
I /.'
1. (b)
Fig. 13.23. Voltage and reactive power at the rectifier for the d.c. line fault. (a) A.c.
terminal voltage; (b) Reactive power flow
Benmore Haywards
(Rectifier) (Inverter)
However, as the matching between the QSS and TCS models is improved, the
differences in peak swing angle are reduced. For curve 2 of Fig. 13.22, the
difference is 2 percent and for curve 3 it is less than 1 percent. The maximum swing
angles for these cases are included in Table 13.4.These differences are particularly
small and are due in part to the improved matching between models and also to
the fact that the peak swing of both machines occurs very close to the d.c. fault
clearance time (0.43 sec. and 0.41 sec. respectively). Hence, any differences during
the link restart period, i.e. ramped or instantaneous in the case of the QSS model,
have little effect on the peak swing angle.
Conclusion
The three case studies presented in this section demonstrate the feasibility and
effectiveness of combining both the TCS and TS programmes. Together they
provide more accurate modelling of a.c.-d.c. systems disturbances than either
programme can provide on its own. Detailed information can be obtained from
the TCS regarding fault development, d.c. link response, converter performance,
over voltages, over currents, e t ~ . ( ~ l))" l
The different results obtained from the TS analysis of the three cases
investigated, indicates that a QSS model cannot be classified as giving either
pessimistic or optimistic results for peak swing angles. This emphasizes the need
to perform worst case transient converter simulations during investigations of a
system so that the QSS link model performance can be modified to more
accurately represent d.c. link performance for further studies. The system used for
these case studies is relatively strong and a weaker system, where stability is more
marginal, may produce much greater differences between the two models.
13.8. REFERENCES
1. W. Shepard and P. Zakikhani, 1972. 'Suggested definition of reactive power for
nonsinusoidal systems', Proc. IEE, 119 (9), 1 361 -1 362.
2. D. Sharon, 1973. 'Reactive power definitions and power factor improvement in
nonlinear systems', Proc. IEE, 120 (6), 704-706.
3. E. Micu, 1973. 'Suggested definition of reactive power for nonsinusoidal systems',
Proc. IEE, 120 (7), 796-797.
4. B. P. Lathi, 1968. Communication Systems. Wiley and Sons, New York.
5. W. T. Cochran, et al., 1977. 'What is the fast fourier transform? Proc. IEEE, 55 (lo),
1664- 1674.
6. F. J. Harris, 1978. 'On the use of windows for harmonic analysis with the discrete
Fourier transform', Proc. IEEE, 66 (1).
7. K. S. Turner, M. D. Heffernan, C. P. Arnold and J. Arrillaga, 1981. 'Computation of
a.c-d.c. system disturbances-Part 11. Derivation of power frequency variables from
converter transient response.' IEEE Winter Power Meeting, Atlanta, Ga.
8. K. S. Turner, 1980. 'Transient stability analysis of integrated a.c. and d.c. power
systems.' Ph.D. thesis, University of Canterbury, New Zealand.
9. J. Arrillaga, M. D. Heffernan, C. B. Lake and C. P. Arnold, 1980. 'Fault studies in a.c.
systems interconnected by h.v.d.c. links', Proc. IEE, 127, Pt. C (I), 15-19.
10. M. D. Heffernan, K. S. Turner, J. Arrillaga and C. P. Arnold, 1981. 'Computation of
a.c-d.c. system disturbances. Part I-Interactive coordination of generator and
converter transient models.' IEEE Winter Power Meeting, Atlanta, Ga.
11. M. D. Heffernan, J. Arrillaga, K. S. Turner and C. P. Arnold, 1980. 'Recovery from
temporary h.v.d.c. line faults.' IEEE Summer Power Meeting, Minneapolis, Ma.
Appendix I
Numerical. Integration Methods
1.1. INTRODUCTION
Basic to the computer modelling of power-system transients is the numerical
integration of the set of differential equations involved. Many books have been
written on the numerical solution of ordinary differential equations, but this
appendix is restricted to the techniques in common use for the dynamic
simulation of power system behaviour.
It is therefore appropriate to start by identifying and defining the properties
required from the numerical integration method in the context of power system
analysis.
Stability
Two types of instability occur in the solution of ordinary differential equations,
i.e. inherent and induced instability.
Inherent instability occurs when, during a numerical step by step solution,
errors generated by any means (truncation or round-off) are magnified until the
true solution is swamped. Fortunately transient-stability studies are formulated
in such a manner that inherent instability is not a problem.
Induced instability is related to the method used in the numerical solution of
the ordinary differential equation. The numerical method gives a sequence of
approximations to the true solution and the stability of the method is basically a
measure of the difference between the approximate and true solutions as the
number of steps becomes large.
Consider the ordinary differential equation:
with the initial conditions y(0) =yo which has the solution
Note that R is the eigenvalue'') of the single variable system given by the ordinary
differential equation (1.2.3). This may be solved by a finite difference equation of
the general multistep form:
and constraining the difference scheme to be stable when R =0, then the
remaining part of (1.2.5)is linear and the solutions are given by the roots zi (for
i = 1,2,. ..,k) of m(z) = 0. If the roots are all different, then
A system of ordinary differential equations in which the ratio of the largest to the
smallest eigenvalue is very much greater than one, is usually referred to as being
stiff. Only during the initial period of the solution of the problem are the largest
negative eigenvalues significant, yet they must be accounted for during the whole
solution.
For methods which are conditionally stable, a very small step length must be
chosen to maintain stability. This makes the method very expensive in computing
time.
The advantages of &stability thus become apparent for this type of problem as
the step length need not be adjusted to accommodate the smallest eigenvalues.
In an electrical power-system the differential equations which describe its
behaviour in a transient state, have greatly varying eigenvalues. The largest
negative eigenvalues are related to the network and the machine stators but these
are ignored by establishing algebraic equations to replace the differential
equations. The associated variables are then permitted to vary instantaneously.
However, the time constants of the remaining ordinary differential equations
are still sufficiently varied to give a large range of eigenvalues. It is therefore
important that if the fastest remaining transients are to be considered and not
ignored, as so often done in the past, a method must be adopted which keeps the
computation to a minimum.
Basically the methods consist of a pair of equations, one being explicit (BO= 0) to
give a prediction of the solution at tn+, and the other being implicit (Po # 0)
which corrects the predicted value. There are a great variety of methods available,
the choice being made by the requirements of the solution. It is usual for
simplicity to maintain a constant step length with these methods if k > 2.
Each application of a corrector method improves the accuracy of the method
by one order, up to a maximum given by the order of accuracy of the corrector.
Therefore, ifthe corrector is not to be iterated, it is usual to use a predictor with an
order of accuracy one less than that of the corrector. The predictor is thus not
essential, as the value at the previous step may be used as a first crude estimate,
but the number of iterations of the corrector may be large.
While, for accuracy, there is a fixed number of relevant iterations, it is desirable
for stability purposes to iterate to some predetermined level of convergence. The
characteristic root (2,) of a predictor or corrector when applied to the single
variable problem:
py = Ay with y(0) = yo (1.3.3)
may be found from:
When applying a corrector to the problem defined by equation (1.3.3) and
rearranging equation (1.3.2) to give:
the solution to the problem becomes direct. The predictor is now not necessary as
the solution only requires information of y at the previous steps, i.e. at yn for
++
i = l,2, ...k.
Where the problem contains two variables, one nonintegrable, such that:
py = Ay + px with y(0) = yo (1.3.6)
then:
where
and
, ,
Although cn+ and mn+ are constant at a particular step, the solution is iterative
using equations (1.3.7) and (1.3.8). Strictly in this simple case, xn+, in equation
(1.3.8) could be removed using equation (1.3.7) but in the general multivariable
case this is not so.
The convergence of this method is now a function of the nonlinearity of the
system. Provided that the step length is sufficiently small, a simple Jacobi form of
iteration gives convergence in only a few iterations. It is equally possible to form a
Jacobian matrix and obtain a solution by a Newton iterative process, although
the storage necessary is much larger and as before, the step length must be
sufficiently small to ensure convergence.
For a multivariable system, equation (1.3.1) is coupled with:
and the solution of the integrable variables is given by the matrix equation
K + I = Cn+1+ M n + l . C Y n + d n + l l t (1.3.12)
,
The elements of the vector Cn+ are as given in equation (1.3.9) and the elements
of the sparse Mn+, matrix are given in equation (1.3.10).
The iterative solution may be started at any point in the loop, if Jacobi iteration
is used. Because the number of algebraic variables (X) associated with equations
(1.3.1) or (1.3.11) are small, it is most advantageous to extrapolate these algebraic
variables and commence with a solution using equation (1.3.1 1).
The disadvantage of any multi-step method (k > 2) is that is not self-starting.
After a discontinuity k - 1, steps must be performed by some other self-starting
method. Unfortunately, it is the period immediately after a step which is most
critical as the largest negative eigenvalues are significant. As k - 1 is usually
small, it is not essential to use an A-stable starting method. Accuracy over this
period is of more importance.
where
k,= hf(tn+cih,yn+ x
j= 1
aijkj) for i = 1,2,. . . , u (1.4.2)
Being single-step, these methods are self-starting and the step length need not be
constant. If j is restricted so that j < i, then the method is explicit and c, must be
zero. When j is permitted to exceed i, then the method is implict and an iterative
solution is necessary.
Also of interest are the forms developed by Merson and Scraton. These are
fourth-order methods ( p =4) but use five stages (u = 5). The extra degree of
freedom obtained is used to give an estimate of the local truncation error at that
step. This can be used to automatically control the step length.
Although they are accurate, the explicit Runge-Kutta methods are not A-
stable. Stability is achieved by ensuring that the step length does not become large
compared to any time constant. For a pth order explicit method the characteristic
root is:
where the second summation term exists only when v > p and where ai are
constant and dependent on the method.
For some implicit methods the characteristic root is equivalent to a Pade
approximant to eh'.
The Pade approximant of a function f ( t ) is given by
and i f
then
Morsden
Bunnvthoroe II I
Brunsw~ck
Hoywords
D.C link
Smelter
Tiwai
The data presented here represents the input and output of a transient stability
study of the two primary a.c. systems plus the asynchronous h.v.d.c. in-
terconnection. The a.c. fault is in the South Island (rectifier) system and the
kv.d.c. link is assumed controllable throughout the study.
The loading data in the input is the solution from a load-flow program which
may be an integral part of the transient stability program or, as in this case, be run
separately. Also note that output of the AVR and speed governor information has
been restricted to two synchronous machines each.
Examples of the data as required for a short-circuit analysis program and for
various load-flow programs are given in the relevant chapters.
33333 >033) 13553 >0=30 a355
000J> 30391 03>03 07300 3303
...............
33333 32-35 33333
........................
>>>a> >a>>> 3>333
30031
133,)
0a331
3033
213) ZJ 3
3 3
1 >
-0331
.........................................
-73aU 30-)3
1>>31 33353
=333>
=>a33
3353>
a3333
1 > 3 3 1
3>33>
00
33
- -
C I 0 3 0 OOO=O 01331 00033 00300 00000 00303 00
33a23 ..91d2 0--39 33000 0030
a 3 31
U
d
3a033 30350 00333 030013 -000 Z 0 0 f C O G 0 tO\lU)O OZrOO Nr(?0O OU3-0 O-5-C NO
00300 00000 OOJOO OODOO ooon 4 - o o o o l c o o ~ o m o as o - r o ma-a- no o o o n m m o
33337 03300 0>300 00-oo 0000 9-3 3 o 7 1U3UO ? f a o r y w m n m -mr-v fa--= --==a -3
........................
3311'3
a3372
03a=o a3301
35333
z.aaaa
n033U
0030
3333
1. 3
&,1
...................................
3
o
1
CJ
a
~ ~ - nnr-n
0 3
2J-DU
n c r r a an*=-
17713
o . r a r a aaao3 oo
-3UUa O X 3 a 3 0=
-
53333 33303
-
J
33330 01333 00331 N300, 0003 Y) 0 0 0 03003 00033 0a000 03030 >5331 030- 30
U 0 00 0
Y 3
v n?
30, 3 0 3 3 3 n 3 3 3 3 33333 03533 03=03 00033 0 3 3 0 3 3 3 3 3 3 003-3 33333 0 3 3 0 3 3
onn - - ~ n - 3,-nn tr-17 ~ ~ a n> ~ t q t tOD>-- >>CVI n.-7, 3332- - 3 5 - 1 --u=n 7
U
t3-
>>a
on-
O
--a33 -----
-~ c u n na-n-- --t--
33333 3333-
OCOOO 0 0 0 0 ~ ooOo0
t t n o t U.-UU
0 0 0 ~ 0O O ~
0 0 3 0 0 n ~ q u f ir - - ~ -
...........................................................
171-1 ->3--? 11cc0 3 3 1 1 0 = o n > > ---DO
,-
U U % > ~0 - - 0 0
3 3 2 3 3 3 3 3 > 3 3 3 3 > 3 ---=a
~ o
00000 00000 00000
1177-
0--,3
>:>a-
00303
>ton- nnot-
a ~ n - - nrno?
-7333 33--3
03000 O O ~ O Oo
n
n
3
I
m
I
I
I
I
. I .
IYI
I
.=;:
V) I 0
W 3 r
3 n l r o o ~ ~ o o o o l o o m m ~ o o 3 a O9.-
m
I1
I
A, I
-.
0
3
1
I
1I
0
11
W
5
aaYI
x2
.IC
P-J 2xxxsssXLsauaaauuQ.auaau
DL
"la
32
0
X d
OX
X 3
- 2Z
W
>W
% II f C
",<a
XI.
4x2
a I
3 I
.r
LL I
I
3 :I -.
U
" I P
-a
v
II
~
W
n
N
I Il i u)
I I
I I
J I
-:I
V
rl II
I
>
*
- 1
I 3.
-
t
.
l
I I
Z
I I:
I W
2 1 -
3 1 -
4 1 -
P I " l
I
I
I
I
I .
.
I "0
I
" I *
a
0
:1
'-!
0
t
N
. lI
I
U
a l e
a
'Y
1
1
02 11 -
4L 1I .
-' II "
W I U
-41 II u
U
I I
I I
W.,
J I
II D
- 1 7 .
--
3 1
-&- 1I .
1 3
c 1 z
I
II !,
2 1 -
3 1 0
- 1 -
r l 'A
-
3 c
--. -
3
I- -a
-
0
rc
3.-
0
0
a0
J3
%
111-
Id
-4
a", 3 m
-aC3 u
J
3
J a
3-
1
x-z
X J
*
C
I
c
z
x
J
-
J
-LI --
Z
A!
**J
3.72
0
I+
xc
J)
L
UJ 30
- .. -.
4 4
SZ
23
-0
00
00
An
h.
I
I
w-
>x
J2
4 ..
L
-
-
-J
A 4
J3
5 .
L
00
33
-.-.
f<
U
\23
4-
z
rn 3
-
- n..>
m N
Y . u-
T
3 t N
A! J
C >J
- ..
JU >- JZ
>.
-.
-L > 0 -7
L9
% J)
00 21
3
2 1
-
3
un
m-t
om
dc
JJ 3.5
1- 3.0
1 93.
m
J
a
..
31
,>
01
T
C
..
a-
3.-
03
c-
3, ..
01
-a
00
z J)
un
'3 3 2% 00
7- 1 I
w"3
35 J"l
u
.o .
T2
4-
0-
un
..
OD
nu7
-
-Q 0I %
I
-4
1 I u
.3 .
.7 I 1
.*
9Z W
I.,I LL
"34. 3 r Q
.-
-
3 2 S3 44
.U WT
I I
L 1: I 1
4 n
31d D I fd. d
*PA! l
1114
32
0
a
-
o x
I 3
AJr
Wn
z
Y
I I
'1111
>d
.*
C I
d3 O-Z
0
0-
n
-,-I
C U
-l ax
u
WC
-,
u x
31d at
m
c1Q.O-
rl-o-
I
:!
:
:I
50'0-
on-o-
Subject Index
A-stable, 23 1,397 Characteristic root, 23 1,398
A.C. faults, 180, 184, 189,228, 347 Cornmutating voltage, 46, 129
A.C. load flow, 58.90 Commutation, 45, 129
angle reference, 59 abnormal, 307
balanced, 58 failures, 387
power mismatches, 65, 7 1,94 normal, 307,309
unbalanced, 9 1 Commutation angle, 46
voltage reference, 59 Commutation reactance, 46, 50
A.C. system model Compound turbine, 261
impedance angle, 328 Connection matrices, 6, 13
impedance frequency locus, 327,347 compound matrices, 15
Thevenin equivalent, 282,328,338, three phase transformers, 28
345,347 Constant current control, 53, 277,3 18,
time-variant, 328,338,367 354
A.C.-D.C. conversion, 43,271 Constant extinction angle control, 53,3 18,
A.C.-D.C. load flow, 124, 152 359
control equations, 129, 159 Constant power control, 55,277
converter variables, 126, 156 Contactors, 293
mismatch equations, 125 Convergence, 8 1, 14 1,229
per unit system, 128 Converter
sequential algorithm, 124, 137 commutation (source) voltage, 46
unbalanced, 153 commutation reactance, 46,50
unified (simultaneous) algorithm, 124, controlled rectification, 45
13 1 faults, 344
Admittance matrices, 5 inversion, 48
formation for simple networks, 17 modes of operation, 273,307
formation for three phase networks, 4 1 parallel connection, 5 1
three phase compound, 13.21 phase angle control, 152
Angular momentum, 208 series connection, 50
Applicability of A.C.-D.C. models, 381 six-pulse, 309
Automatic voltage regulator, 95,214,238 transient model, 306
twelve pulse. 309
Backward Euler method, 233 uncontrolled rectification, 44
Boiler, 2 18 Converter control
Branch switching, 228 characteristics, 54,272
constant current, 53,277,3 18,354
Cage factor, 268 constant power, 55,277
firing instant, 3 18, 363 interactive coordination, 336,356, 382
firing limits, 3 18 nodal analysis, 30 1, 345
power modulation, 279 partitioning, 302,304
Converter transformer, 44, 306 state variables, 330, 335
connection, 162, 1 74 topological matrices, 30 1
leakage reactance, 44,307,3 1 1,323 valve voltages and currents, 3 13
magnetization curve, 324 variable topology, 303
magnetizing reactance, 324 voltage and current relations, 302,3 13
on-load tap changer, 44,48, 127, 155, Dynamic stability, 206, 26 1
323
phase-shifting, 43 Eigenvalue, 234,397
saturation, 324 Enforced delay, 274
tertiarv winding, 5 1 Equidistant firing control, 152, 159,3 17
transient model, 323 Exciter, 2 16, 240
Cylindrical rotor, 2 10
Fast decoupled load flow
D.C. faults, 353, 360, 364, 389 assumptions, 80,98
arc extinction and deionization, 357, balanced A.C. system, 78
360,389 balanced A.C.-D.C. system, 13 1
detection, 354, 356, 360 characteristics and performance, 83,
earth resistance, 358,36 1 105, 144, 167
restart, 359, 362, 368 flow diagram, 8 1, 101, 135, 165
D.C. transmission, 2, 53, 124, 154, 168. multiterminal D.C., 139
32 1 sequential solution technique, 137
monopole, 329 simultaneous solution technique, 131
multiterminal, 57, 139 unbalanced A.C.-D.C. system, 152
New Zealand link, 5 1,329, 357 unbalanced A.C. system, 96
n-equivalent, 329 unbalanced system program structure,
quasi-steady state, 277, 337,345, 382, 103, 166
389 unified solution technique, 13 1
steady state, 43 Fault calculations, 184, 189
submarine cable, 329,357 Fault studies, 180
two-terminal, 53 Filters, 46, 5 1, 155,305,326,367
Dead band, 286 Fourier analysis, 47
Decoupled load flow methods, 76 fundamental current, 47, 160
Deep bar rotor, 267 spectral analysis, 370, 372
Delay angle, 45, 53, l52,27 1 spectral leakage, 373
Direct axis, 207 TCS waveforms, 3 7 1
Direct current, 44, 53 Frames of reference, 207,224,267
margin, $5
ripple, 44 Harmonics, 156, 177,353,367,370,374
setting, 55 Hydrogovernor, 2 17
Distance relays, 293 Hydroturbine, 2 18,26 1
Double cage rotor, 267
Dummy busbar, 229 Induction machine, 265
Dynamic A.C.-D.C. model, 300 cage factor, 268
A.C. faults (see A.C. faults) contactors, 293
branch equations, 302 deep bar rotor, 267
computer program, 332, 340 double cage rotor, 267
D.C. faults (see D.C. faults) inertia constant, 266
dependent variables, 303,335 magnetizing reactance, 267
initial conditions, 3 15 mechanical torque, 266
Induction machine (contd.) Nodal formulation
open circuit reactance, 267 dynamic A.C.-D.C. model, 300
slip, 265 load flow, 60
transient reactance, 26 7 transient stability, 223
Inertia constant, 208,266 Nonunit protection, 229
Infinite machine, 214 Norton equivalent, 63,224
Integration, 230,330, 395 Numerical stability, 33 1,396
A-stable, 23 1,397
backward Euler method, 233 Open circuit faults, 192
characteristic root, 23 1, 398 Open circuit reactance, 267
implicit, 33 1 Ordering of sparse matrices
Runge-Kutta methods, 230, 331,400 dynamic, 74
X-stable, 23 1, 397 pre-ordering, 74
step length, 23 1,33 1 Overcurrent relays, 290
trapezoidal method, 23 1
Interactive coordination, 336,356,382 Per unit for D.C. system, 48, 128, 157,33 1
Interceptor governor, 265 Potier reactance, 252,257
Interceptor valve, 26 1 Power factor in A.C.-D.C. systems, 49
Inversion, 48 Power-frequency control, 56
characteristic, 54 Predictor-corrector methods, 398
deionization angle, 48 Primitive networks. 5
extinction angle, 48,53
Iteration schemes Quadrature axis, 207
A.C. load flow, 72,8 1, 101
A.C.-D.C. load flow, 132,165 Reactive power in A.C.-D.C. systems, 50,
transient stability, 247 53
Rectification, 44,27 1
Jacobian matrix elements abnormal modes of operation, 273
A.C. load flow, 68 characteristics, 54,272
A.C.-D.C. load flow, 132, 167 controlled, 45
delay angle, 45,53,271
Lead-lag circuit model, 2 19 dynamic load, 273
Linear transformation, 4, 15,29 enforced delay, 274
Load characteristics, 220 operating mode identification, 276
Loads, 64,220,326 saturable reactors, 27 1
static load, 272
Machine switching, 230 uncontrolled, 45
Magnetizing reactance, 267 Relays, 290
Matrix partitioning, 14,24 distance, 293
Matrix sparsity, 60 induction motor contactors, 293
Mechanical power, 208,218 overcurrent, 290
Mechanical torque, 266 undervoltage, 293
Resonance, 354
Network subdivision, 18 Runge-Kutta methods, 230,33 1,400
New Zealand Electricity power system,
329,402 S.C.R. (see Short circuit ratio)
Newton-Raphson S.V.S., 286
characteristics and performance, 75 deadband, 296
equations, 67 Saliency, 2 10,225
flow diagram, 72 Salient pole rotor, 2 10
general formulation, 66 Saturable reactors, 271
Jacobian elements, 68 Saturation, 251,307,324,353
starting techniques, 75 Saturation curve, 253,257
Saturation factor, 252 Synchronous machine models. 208.21,4.
Sequence components (see Symmetrical 237
components) Synchronous reactance. 2 10
Series element three-phase representation. System damping. 279
11.26
Short circuit faults, 189 TCS. 300
Short circuit ratio, 354 alignment with TS, 377, 383
Shunt element three-phase representation, fault application and removal. 347. 378
25 interaction with TS, 336.369,376
C-stable. 23 1. 397 synchronous machine model. 3 19
Slip. 265 transformer model, 323
Smoothing reactor, 44,308 Thermal turbine, 2 18
Sparse matrix equations, 74, 167, 184, 303, Thevenin equivalent. 223,345
340 Three phase faults. 18 1. 228
Sparsity, 184, 303 Three phase representation, 1 1, 18.26
ordering, 74 Transformer model
programming, 72, 167 saturation, 324
Speed governor, 2 17 single phase, 7
interceptor. 265 three phase, 26, 162
Stability, 206, 25 1 (see also Converter transformer)
numerical, 396 Transient converter simulation (see TCS)
representation of plant in network, 224, Transient reactance, 2 1 1, 267
270,28 1,289,295 Transient stability, 206
Starting techniques for load flow, 75 Transient stability program structure, 240
Static VAR compensation systems (see Transmission line model
S.V.S.) ABCD parameters. 39
Stiffness, 207, 397 line sectionalization. 39
Subtransient reactance, 2 12,269 long lines, 22
Swing curves. 388 mutual coupling, 22
Symmetrical components, 10, 18 1, 188,296 single phase, 6
three phase transformers, 35, 1 76 three phase, 18
Symmetrical firing control (see Equidistant Trapezoidal method, 23 1
firing control) Triangular factorization, 73, 136
Synchronous machine Turbine
angular momentum, 208 compound, 26 1
cylindrical rotor, 209 hydro, 2 18,26 1
damper windings, 208,2 12 thermal. 2 18,26 1
damping coefficient, 209,235
inertia constant, 208
negative sequence braking, 296 Unbalance
negative sequence impedance, 296 effect on converters, 152
Potier reactance, 252,257 source of unbalance, 90
salient pole rotor, 2 10 unsymmetrical A.C. system, 12
saturation, 25 1 Unbalanced faults, 187,295
saturation factor, 252 Undervoltage relays, 293
subtransient reactance, 2 12
transient reactance, 2 1 1 Valve
Synchronous machine controllers, 2 14, extinction, 3 18,332
238,265 switching, 33 1