Ac DC Power System Analysis Arrillaga, - Jos
Ac DC Power System Analysis Arrillaga, - Jos
Ac DC Power System Analysis Arrillaga, - Jos
AC-DC
POWER
SYSTEM
ANALYSIS
Jos Arrillaga
and Bruce Smith
Preface
Although the early decision to generate electric power at 50/60 cycles and
three phases is practically irreversible, the transmission and utilisation of
energy is not necessarily tied to these conditions. The choice between
transmission alternatives is made on the basis of cost and controllability.
The original justification for HVDC transmission was its lower cost for
long electrical distances which, in the case of submarine (or underground)
cable schemes, applies to relatively short geographical distances. At present,
the controllability factor often justifies the DC alternative regardless of cost,
as evidenced by the growing number of back-to-back links in existence.
The merits of HVDC over AC transmission have been explained in
several books by Adamson and Hingorani, Uhlmann, Kimbark, Arrillaga,
and Paddiyar, listed in chronological order.
In earlier days, the dynamic performance of the DC link was assessed with
the help of scaled-down physical simulators. These provided a reasonable
representation of the converter control and protection functions, but were
very restricted in AC-network representation.
With the expansion of HVDC transmission throughout the world, and
particularly the increasing numbers of interconnections between different
countries, few power systems can continue to escape the effect of this
technology in their planning and operation. Such expansion has encouraged the development of analytical models to represent the behaviour of the
AC-DC power system.
An early attempt to describe the HVDC link as a power system component
was made in the book 'Computer modelling of electrical power systems'.
Although the book's main objective was conventional power-system analysis,
it did propose algorithmic modifications for the incorporation of HVDC
transmission. Since then the experience of many years of HVDC operation
has produced more advanced models to represent the behaviour of both the
AC and DC systems.
In particular, the availability of the EMTP (electromagnetic-transient
program) with detailed representation of power-electronic components and,
more recently, its implementation in the RTDS (real-time digital simulator)
xii
Preface
Contents
Preface
Introduction
.1 Basic AC-DC configuration
.2 AC-DC simulation philosophy
.3 Steady-state simulation
.4 Fault analysis
.5 Harmonic analysis
.6 System stability
The AC-DC converter in steady state
2.1 Introduction
2.2 Power frequencysymmetrical operation
2.2.1 Analysis of the commutation circuit
2.2.2 Rectifier operation
2.2.3 Inverter operation
2.2.4 Power factor and reactive power
2.3 Power frequencyunbalanced operation
2.3.1 Terminology and waveforms
2.3.2 Variables and equations
2.4 Characteristic harmonics
2.5 The converter as a frequency modulator
2.5.1 The modulation process
2.6 Harmonic transfer generalisation
2.6.1 From the AC to the DC sides
2.6.2 From the DC to the AC sides
2.6.3 Effect of switching-instant variation
2.6.4 Transfer across the DC link
2.7 Harmonic instabilities
2.8 Generalised harmonic domain converter model
2.8.1 Analysis of the commutation
Star-connected bridge
Delta-connected bridge
2.8.2 Valve-firing control
Current control
Commutation margin control
XI
1
1
2
3
4
4
5
7
7
8
10
12
14
15
17
17
18
23
26
27
29
30
31
32
33
36
41
41
41
43
44
45
46
vi
Contents
2.8.3 Direct voltage
Star-connection voltage samples
Delta-connection voltage samples
Convolution of the samples
2.8.4 Phase currents
Converter side
System side
2.9 Summary
2.10 References
48
48
49
51
54
54
54
59
60
61
61
61
63
66
70
71
72
74
74
77
78
78
78
83
84
86
87
90
92
93
94
95
98
100
102
102
104
106
107
108
The
4.1
4.2
4.3
4.4
4.5
109
109
109
110
113
116
117
117
121
harmonic solution
Introduction
Basic AC-DC system
Functional notation of the converter equations
Mismatch equations
Newton's method
4.5.1 The Jacobian matrix
Numerical differentiation
Analytical derivation
Contents
4.6
4.7
4.8
4.9
flow
flow
vii
122
123
123
124
124
128
130
130
132
133
134
141
141
143
143
143
145
146
153
153
158
161
162
162
165
167
167
169
171
172
175
175
177
178
180
181
182
184
186
189
192
193
199
201
204
209
215
216
viii
Contents
6.7
217
217
220
223
225
228
229
230
230
233
235
237
241
241
258
259
262
263
269
271
Electromechanical stability
7.1 Introduction
7.2 Dynamic model of the synchronous machine
7.2.1 Equations of motion
7.2.2 Electrical equations
7.2.3 Synchronous-machine controllers
7.2.4 Generator representation in the network
7.3 Load representation
7.4 The AC transmission network
7.4.1 System faults and switching
7.5 Static power conversion models
7.5.1 Single-converter loads
Abnormal modes of converter operation
Converter representation in the network
7.5.2 DC links
DC power modulation
DC link representation in the network
7.6 , AC-DC transient stability programs
7.6.1 Program structure
7.6.2 Trapezoidal integration
7.6.3 Initial conditions
7.6.4 Test of operating mode
7.6.5 Converter program interface
7.6.6 Test for commutation failure
7.7 Test system and results
7.7.1 Minor disturbance
7.7.2 Major disturbance
7.8 Summary
7.9 References
275
275
276
276
278
280
282
285
286
287
287
288
290
293
296
298
301
301
301
301
304
306
307
309
312
315
316
320
320
Contents
8
ix
323
323
324
326
326
328
329
331
331
332
335
339
343
346
351
355
356
356
361
362
Appendices
I
Newton-Raphson method
363
1.1 Basic algorithm
363
1.2 Techniques to make the Newton-Raphson solution more efficient 364
Sparsity programming
365
Triangular factorisation
365
Optimal ordering
366
1.3 References
366
II
367
367
368
370
III
Test systems
III. 1 CIGRE HVDC benchmark model
111.2 Simplified test system
111.3 Test systems used in the stability chapters
System A
371
371
371
373
373
System B
375
IV
State-space analysis
379
Numerical integration
383
VI
Curve-fitting algorithm
387
Index
391
Chapter 1
Introduction
k 4
.z
A^ 3
_
ii 5
Figure 1.1
i^ 6
'f 2
Figure 1.2
Converter
transformers
6-Pulse
bridges
D.C. line
|W
2
1
filters
filters
Introduction
Introduction
Chapter 2
2.1 Introduction
The transfers of voltage and current across the AC-DC converter are
completely specified by the switching instants of the bridge valves, being
both the firing and end of commutation instants. On the assumption of a
balanced, undistorted AC-terminal voltage, and infinite smoothing reactance, the converter is readily analysed by Fourier methods.1 Under these
conditions, closed-form expressions can be obtained for the firing angles,
commutation duration, fundamental-frequency voltage and current, characteristic phase-current harmonics and DC-voltage harmonics.
However, under realistic conditions, there is some asymmetry and
distortion, the switching instants of the bridge valves are not equispaced
over one cycle owing to control action, and the transfer function between
the AC and DC system is modified. Even a small modulation of the
switching instants can lead to current components in the AC system at the
modulation frequency sidebands. This effect is equally important to both
commutation duration and to firing-angle variation.
The incorporation of switching-angle modulation in the converter model
permits an accurate derivation of the individual switching instants; their
effect on transfers between the AC and DC systems can then be quantified,
and all causes influencing the modulation accounted for. An early cause of
firing-angle modulation was the use of individual firing control, which was
responsive to harmonic distortions in the terminal voltage.2 The adoption
of equidistant firing control, with its much longer time constant, effectively
eliminated this type of firing-angle modulation. However, the firing angle
is still modulated as a result of harmonics in the DC current, via the
current-control loop, sometimes resulting in harmonic instability.3 The
commutation duration is modulated by terminal-voltage harmonics, DC
ripple and firing-angle variation.4 All these variations can have a significant
effect and must be included in the converter model if good accuracy is
required. Modulation of the commutation-period duration by terminalvoltage harmonics, in particular, has a significant impact on the AC-side
harmonic response of the converter.
(a)
(b)
(c)
- j _
'j>
(d)
T\ r/
Figure 2.1
10
(c)
Figure 2.2
Under these assumptions, Figure 2.1 shows the three-phase system voltages
and (prefiltered) converter currents of the bridge rectifier of Figure 1.1 for
star and delta-connected bridges. The currents in continuous trace correspond to ideal commutations, and the dotted lines show the effect of
commutation overlaps.
Figure 2.2 shows the AC voltages on the secondary side of the converter
transformer, the rectified voltage on the converter side of the smoothing
reactor and the AC current (phase R) on the converter side of the
transformer.
D
e
11
1^1 -
(a)
'3
(c)
Figure 2.3 The commutation process
a Equivalent circuit of the commutation from valve 1 to valve 3
b Voltage waveforms showing early (rectification) and late (inversion)
commutations
c The commutating currents
eR and eY (and must be completed before the lower crossing of these two
voltages).
Since eY > eR, a commutating current ic (=i3) builds up at the expense of
ix so that at all times
h +H = h
(2.2.1)
As the rates of change of i3 and ix are equal (provided that the commutation reactances are balanced), the voltage drops across XCR and XCY are
the same and, thus, during the overlap period, the direct voltage Vd is the
mean value of eY and eR.
12
From the circuit of Figure 2.3a and assuming XcR = XCY = Xc we can
write
eY - eR = 2(Xc/co) d(ic)/dt
(2.2.2)
(2.2.3)
- L ^ aVtermsin(cot)d(cot) = Xc f * d(ic)
(2.2.4)
[cos a - cos(cot)]
(2.2.5)
(2.2.6)
(2.2.7)
13
where Vc0 is the maximum average DC voltage (i.e. at no load and without
firing delay); for the three-phase bridge configuration Vc0 = (5y/2/n)aVteTm
and aVterm is the phase-to-phase r.m.s. commutating voltage.
Eqn. 2.2.7 specifies the DC voltage in terms of aVterm, a and fi. However,
the value of the commutation angle is not normally available and a more
useful expression for the DC voltage, as a function of the DC current, can
be derived from eqns. 2.2.6 and 2.2.7, i.e.
Vd = ^L
aVterm cosa
- ^ 1 ,
(2.2.8)
I J
Irms
w P'P I2d(cot)\ = y/njj*
rms = Jhw
(2.2.9)
Since harmonic filters are assumed to be provided at the converter terminals, the current flowing in the AC system contains only the fundamental
component frequency and its r.m.s. magnitude (obtained from the Fourier
analysis of the rectangular waveform) is
(2.2.10)
If the effect of commutation reactance is taken into account and using eqns.
2.2.5 and 2.2.6, the currents of the incoming and outgoing valve during the
commutation are defined by eqns. 2.2.11 and 2.2.12, respectively
L,(cos a cos cot)
_j*
L for a < G * < a + ix
(2.2.11)
l =
cos a cos(a + fi)
.
cos a cos(cot - 2n/3) _
2%
2n
i = Id-Id
^
- ^ fora + <cot<(x + - + // (2.2.12)
3
3
cos a - cos(a + fi)
In between commutations, the current is
i = Id
(2.2.13)
14
where
k = V{[cos 2a - cos 2(a + fi)]2 + [2/x + sin 2a
- sin 2(a + /*)]2}/{4[cos a - cos(a + //)]}
(2.2.15)
and taking into account the transformer tap position, the current on the
primary side becomes
(2.2.16)
cof
Figure 2.4
15
Moreover, there is a fundamental difference between rectifier and inverter operations which prevents an optimal firing condition in the latter case.
Although the rectifier delay angle, a, can be chosen accurately to satisfy a
particular control constraint, the same is not possible with respect to angle
y because of the uncertainty of the overlap angle fi. Events taking place after
the instant of firing are beyond predictability and, therefore, the minimum
commutation margin angle, y0, must be sufficient to cope with reasonable
uncertainties (values between 15 and 20 are typically used).
The analysis of inverter operation is not different from that of rectification. However, for convenience, the inverter equations are often expressed
in terms of the commutation margin angle y(y = /} n, where /? = n a).
Thus, omitting the negative sign of the inverter DC voltage, the following
expressions apply
(2.2.17)
The expression for the direct current is
l
aV
cos
(2.2.19)
and
cos <t> = VdId/(V3aVtermI)
(2.2.20)
Substituting Yd and ld from eqns. 2.2.7 and 2.2.10 into eqn. 2.2.20, the
following approximate expression results
cos (j) = (|)[cos oc + cos (a + /i)]
(2.2.21)
16
The reactive power is often expressed in terms of the active power, i.e.
Q = P.tan<
(2.2.22)
t a n <f> =
(2.2.23)
(2.2.24)
Figure 2.5
17
2 /2
Primary
Figure 2.6
Secondary
18
////#////
Phase 3
Phase 2
both sides of the converter. The phase-neutral voltage is used as the base
parameter and, therefore
MVAbase = base power per phase
Vbase = phase-neutral voltage base
The current base on the AC and DC sides is also equal. Therefore, the p.u.
system does not change the form of any of the converter equations.
The phase-to-phase source voltages referred to the transformer secondary
are found by a consideration of the transformer connection and off-nominal
turns ratio. For example, consider the star-star transformer of Figure 2.7.
The phase-to-phase source voltages referred to the secondary are
--i-V 1
"~
^ ter
term /^term
tet
(2.3.1)
(2.3.2)
(2.3.3)
which in terms of real and imaginary parts yield six equations.
19
Phase 1
Phase 3 / ^ \
Phase 2
(a)
(b)
wmm
I n
(c)
jj
j i l l j i ill
jj!
iji
I
C3+OL
Figure 2.8
(2.3.4)
(2.3.5)
The three-phase converter transformer is represented by its nodal admittance model, i.e.
* node
Ypp
YpS
YSP
Yss
where p and 5 indicate the primary and secondary sides of the transformer
on the assumption of a lossless transformer (i.e. Ypp =jbpp etc.).
Turning to the AC side, on the assumption of a lossless transformer, the
currents at the converter-side busbar are expressed as follows
(2.3.6)
By subtracting the 6lerm 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
- C,)]
(2.3.7a)
k= l
(2.3.76)
(2.3.7c)
(2.S.7d)
k=l
(2.3.7)
(2.3.7/)
fc=l
21
Three equations are derived from approximate expressions for the fundamental r.m.s. components of the line-current waveforms as shown in Figure
2.8c, i.e.
(2.3.8a)
li|sin(T 2 /2)
(2.3.86)
l3 = - - - k s i n ( T 3 / 2 )
(2.3.8c)
y/y/i
II
y>/s
P
1
-ylj*
y/y/3
~P
-y/y/3
-v
1
-v
2
~P
1
2
3^
1
V
v
terir
V2
V3
v
terir
(2.3.9)
22
The sum of the real powers on the three phases of the transformer
secondary may be equated to the total DC power, i.e.
Vd.Id = X %ih cos(4>t ~ ,)
(2.3.10)
(2.3.11)
Efsin<^ = 0
(2.3.12)
(Xcl
c o s / / + cosa 2 - ld
(Xc
+ Xc3)
=0
(2.3.13a)
ltXc3) = 0
(2.3.13*)
V2U13
V ^
^23
(2.3.13c)
Two further equations result from assuming that the off-nominal tap
position is the same in the three phases, i.e.
a1=a2
(2.3.14a)
a2=a3
(2.3.14*)
23
(2.3.15a)
VA, = 1 7
(2.3.156)
or
11
)
(2.4.1)
+ cos 13co^ - cos 17co^ + cos \9cot
From eqn. 2.4.1 the r.m.s. magnitude of the fundamental frequency is
V^ I ( ) = v^ I(j
(2A2)
24
2v/ 3
/
1
1
1
ld I cos cot + - cos 5cot - cos 7 cot cos llcot
n
\
5
7
11
+ cos \3cot + cos I7cot
cos \9cot
(2.4.3)
TC
(2.4.4)
This series only contains harmonics of order \2n 1. The harmonic currents
of orders &n + 1 (with n odd), i.e. n = 5, 7, 17, 19 etc., circulate between the
two converter transformers but do not penetrate the AC network.
The effect of commutation overlap, represented by eqns. 2.2.11 and
2.2.12, is to reduce the magnitude of the harmonic currents but does not
alter the half-wave symmetry of the waveform and, therefore, the harmonic
orders remain the same.1 Figure 2.9 illustrates the variation in 11th
harmonic with delay and overlap angles.
The DC voltage contains the following three different functions with
reference to voltage crossing Ct in Figure 2.2
vd = y/2aVXmQos\cot
+I
vd = >/2Vtermcos cot + I
= ^ aV term cos [cot]
vd = y/t aVterm
cos II mt
mt -- | |
term cos
for
(2.4.5)
+ -J2aVtermsm[cot]
for
(2.4.6)
(2.4.7)
25
10
p.
>
1'
E
c
3
6\
54-
a 40^vv^X \
or - SO^X
a -
o-
10
20
30
40
Figure 2.9
v2c0
J2(k - 1
(k + I) 2
- I) 2
l)^l cos(2a +
(2.4.8)
Figure 2.10 illustrates the use of eqn. 2.4.8 to derive the variation of the
12th harmonic as a percentage of Vc0, the maximum average rectified
voltage, which for the six-pulse bridge converter is 3y/2aVterm/n. For a = 0
and fi = 0, eqn. 2.4.8 shows that
ir
V c0
/o
/Q
(2.4.9)
k2 - 1 ~ k2
26
30
Angle of overlap JJ.
Figure 2 JO
This means that the higher harmonics increase faster with a. Eqn. 2.4.10
represents the maximum proportion of harmonics in the system, particularly when it is considered that at a = 90, /i is likely to be very small.
27
V = Y #AC .i DC
(2.5.2)
where Y , ^ and Y^c are transfer functions for the voltages and current,
respectively, and \// = 0, 120, 240 for each of the three phases.
28
Y,, =
1 A, c o s n q t
n=l
co
+ $)
(2.5.3)
(2.5.4)
00
Y,, =
A, cos(w,t
n= 1
where
This process is illustrated in Figure 2.12, which shows the modulated output
current on the AC side of the converter in response to a DC current which
contains a ripple frequency.
29
I,, /phase
*t
u+
When the commutation angle is taken into account, the voltage transfer
functions (Figure 2.1 la-continuous line) are shown as in eqn. 2.5.3 but
with the value of A, replaced by
(2.5.5)
sin ( 4 2 )
nd2
(2.5.6)
For an accurate quantitative assessment, the analysis must include the effects
of pulse-position and pulse-duration modulation as affected by the modulating frequencies as well as explicit representation of the converter controller.'
30
AC?
DC
AC,*
AC!
(12^,U
ir,t
{ 1 2 n * 1 ) f , / r , -1
(12m+1)W^2
12n
(12m-1)(12/>'*i)/1 /i
-!
(I2m+1)^
(12m+1)(12ny
(12m-1)(12nU)
Figure 2.13 Harmonic transfers across a 12-pulse HVDC link. The encircled
elements indicate harmonic sources and m,ne(l,
2,3...)
31
forne(0, 1,2,3,...)
(2.6.1)
32
= ((12m
f 1) f k ) f ,
for mG(0, 1, 2, 3 , . . .)
(2.6.2)
15
Figure 2.14
20
25
30
35
40
33
45
50
34
(2.6.3)
35
(n = 0,1,2,...)
(2.6.4)
(2.6.5)
(2.6.6)
One of these frequencies (ft 2/2) will beat with the fundamental-frequency voltage of system 1 at a frequency
/2)
(2.6.7)
which for a 50/60 Hz conversion scheme becomes 20 Hz. This is a flickerproducing frequency. This same frequency will be referred to generator
rotor-shaft torque at 20 Hz, which may excite mechanical resonances.
Again, this type of crossmodulation effect is most likely to happen in
back-to-back schemes due to the stronger coupling between the two converters, although it is also possible with any HVDC scheme in the presence of
a suitable resonance. However, for the McNeill back-to-back scheme in
Alberta, which has no DC smoothing reactor, the calculated and measured
effects were small and, therefore, buried in the continuous random variations found in the power system.
Now consider two AC systems of the same nominal (but slightly different)
frequency.
Substituting m n = 0, and k = 2, for fundamental frequency f0 into eqn.
2.6.5, a current and resultant voltage (through the AC-system impedance)
of frequency
/ A C I = / O 2 ( / O + A/ O )
(2.6.8)
which leads to/ 0 2A/0 is induced on the AC side. This will either beat with
the fundamental frequency f0 at frequency 2A/0 or produce generator/motor
shaft torques at 2Af0. This frequency is generally too low to produce flicker
but may induce mechanical oscillations.
36
convertor
switching
action
positive-sequence
second-harmonic current
distortion I
convertor
switching
transformer
action
core
saturation negative-sequence
dc distortion lacn
AC side
37
fundamentalfrequency voltage
distortion
DC-side
fundamentalfrequency
impedance Zdch
fundamentalfrequency current
distortion l d c h
DC side
positive-sequence
second harmonic
fundamental
frequency
acp
negative-sequence
DC harmonic
Va
Racn
AC side
converter
DC side
Figure 2.16 Equivalent circuit for study of transformer core saturation instability
38
inductance, Lm, and AC-side impedance at the low frequency IACn variation.
Assuming that the AC-side impedance remains fairly constant around 0 Hz,
the impedance can be simplified to the AC-system resistance at 0 Hz of
RAOl. On the DC side, the presence of fundamental-frequency voltage
distortion causes an equivalent current distortion to flow through the
DC-side impedance.
The converter block accounts for the transfer of AC-side voltage distortion to the DC side and the DC-side current distortion to the AC side, as
described in section 2.6. Considering only the most significant low-order
harmonics, these interactions can be represented by eqn. 2.7.1. The elements within the matrix describe the amplitude change and phase shift
introduced during the transfer of distortions from either side of a converter.
The response of the converter controller and the signal transducers is
included in the matrix
^DCh
^ACp
= d
h i
J^ACn __
vAC;
b c
*ACn
(2.7.1)
'DCh
Implicit in this expression is that the VACp and IACp are the higher-frequency
positive-sequence components on the AC side at one harmonic higher than
the DC-side components of VDOl and IDC/I, which are in turn one harmonic
higher than the negative-sequence components of VACn and IACn on the AC
side. In the analysis of core saturation instability, VACp and IACp are at the
second harmonic, Vj^h a n d IDC/I a t the fundamental frequency and VACw and
IACll at DC. The elements within the matrix are calculated from average
values of the converter operation in the steady state. The responses of the
converter controller and the signal transducers are incorporated. The
response of a constant-current controller is included in elements c, f and i
which accounts for the transfer to the three left-hand-side harmonics
through the controller loop.
In addition to the linearised converter transfer function, eqn. 2.7.1
describing the converter operation, a set of simultaneous complex equations
is written to describe this equivalent circuit
YACp I A cpZ ACp + 1 2 + Z AC|J
(2.7.2b)
I,+ = - X L ,
X
(2.7.2fl)
ACn
ACn
dt
(2.7.2c)
(2.7.2d)
(2.7.2e)
39
V* ' ' 3 a /
RAC
1 4- AZACp + XCZ AQ ,
Lm (1 +R A C w D)(l + AZ A C p )~R A O l BCZ A C p
**
A= d+ ^
(2.7.3c)
= *+ ^
(2.7.M)
^
D= A + ^
(2.7.3.)
(2.7.3/)
The term a in eqn. 2.7.3a is defined as the saturation stability factor (SSF)
which basically defines the susceptibility of the system to the development
of this type of instability. A positive SSF value indicates that the harmonic
sequence will decay over time and, hence, the system is stable. On the other
hand, a negative SSF suggests negative damping and the development of
instability as the distorting harmonic sequence increases over time. The
value and sign of the p term determine the speed and the direction of the
variation of these harmonic sequences.
This saturation stability-factor technique has been verified against dynamic simulation,10 and is included as an example of electromagnetictransient simulation of HVDC schemes in Chapter 6, section 6.8. In that
section, electromagnetic simulation is used to verify several solutions to a
core saturation instability, derived by consideration of the SSF.
From Figure 2.15, the contribution from the transformer saturation to the
instability feedback loop comes in the form of an additional positivesequence second-harmonic current resulting from the saturation. This effect
is modelled as a positive-sequence second-harmonic current injection I 2 +
into the second-harmonic part of the equivalent circuit. The depth of the
saturation is calculated according to the amount of negative-sequence DC
from the converter that is flowing into the transformer magnetising inductance, hm (i.e. I o _).
The transformer saturation related to this instability can be regarded as
asymmetrical saturation since the transformer is only saturated in one half
of a fundamental cycle. As the transformer is coming in and out of
saturation, the magnetising inductance, LM, is nonlinear, but it is only in
saturation for a short period of each cycle as indicated by the width of the
distorted magnetising current pulse. Therefore, it is reasonable to assume
40
0<X<l
(2.7.4)
a low and predominantly capacitive DC-side impedance at the fundamental frequency with the presence of a series resonance near to but
higher than the fundamental frequency;
a high and predominantly inductive AC-side second-harmonic impedance with the presence of a parallel resonance near to but higher than
the second harmonic frequency;
a high AC-side resistance near 0 Hz.
41
In the commutation circuit shown in Figure 2.17, Vfl, Vb, lc and Id are sums
of harmonic phasors. In this diagram phase a is commutating off, phase b
is commutating on and the commutation ends when Ic = ld. Considering a
Figure 2.17
periodic steady state, and summing voltage drops around the commutation
loop at harmonic order k
Vak +jkXaIck -jkXaldk
+jkXblck
- Vbk = 0
(2.8.1)
'
- vabk
jk(xa + xb)
(2 8 2)
- -
(2.8.3)
where
(2.8.4)
.k = l
43
loop. Assigning this value to D ensures that at the moment of firing a valve,
6i9 the current in it is zero.
The resulting commutation current is the steady-state solution of the
commutation circuit and, since there is no resistance in the circuit, the
steady state is achieved instantaneously when the appropriate valve is fired.
The commutation ends when the instantaneous commutation current is
equal to the instantaneous DC current. This angle, the end of commutation
0t., cannot be solved directly, as eqn. 2.8.3 is transcendental. Instead, the
end of each commutation is determined by the zero crossing of a differentiate mismatch equation, solvable by Newton's method. The mismatch
equation is easily constructed by substituting cot = ${ into the Fourier series
for the DC and commutation currents, and taking the difference
(I* - IoJ^j
(2-8.5)
I** + hk - hk = 0
hk ~~ lck ~~ ^dk ~ 0
=0
vck -jiCK*ck + y* = o
V^-AA-V^O
(2.8.6)
44
VckKak
-h]IdkXckXak
(2.8.7)
(2.8.8)
45
1 +jk<oT
With reference to Figure 2.20, it can be seen that firing occurs when the
elapsed angle from the equidistant timing reference is equal to the instantaneous value of the alpha order i.e. a = 0i pt. The equidistant timing
references are represented by /?,. = (i 1)TT/3. The firing-mismatch equation
is, therefore
= 0
(2.8.11)
k=l
This analysis of the firing process is also valid for a bridge connected via a
star-g/delta bridge to the AC system, in which case the equidistant timing
references should be advanced by 30.
controller
current transducer
p
Id
current order
Figure 2.19
Current controller
46
A
a
Figure 2.20
Method of finding the firing instants; the timing instants are assumed
perfectly equidistant (7 3 )
The constant component of the alpha order, oc0, cannot be solved directly
as the PI control has a pole at zero frequency. However, in the steady state
the average delay angle, <x0, takes on a value which causes the DC component of the DC current to be equal to the current order. This requirement
is easily expressed as a mismatch equation with a zero crossing at the
current order
jr
ry
4V
E ^Y
(2 H 19^
where Vfwd is the constant forward voltage drop through a group. This
equation states that the DC voltage, when applied to the DC system, causes
the current order to flow in the DC system. The DC voltage is obtained from
eqn. 2.8.35. Note that eqn. 2.8.12 is not a function of a0. The average delay
angle does, however, feature in the firing-angle mismatch eqns. 2.8.11.
When the converter is solved in Chapter 4, the average alpha order emerges
from the Newton solution.
A similar equation to eqn. 2.8.12 can also be written to describe a
constant-power control
= Pspeoified
- 4Vfwd)I,dO
(2.8.13)
47
yn = dn-4>n
where
y = angle of extinction
3 = voltage-crossover angle
(f> = angle of end of commutation
A common inverter-control mode is extinction angle (or gamma) control, to
minimise reactive-power demand with little chance of commutation failure.
Under normal levels of distortion, the twelve gammas will differ, and it is
the smallest of these which is controlled to a specified minimum, typically
15 to 20 degrees. As the switching with the minimum gamma does not
change from cycle to cycle when in the steady state, the firing order remains
constant and unresponsive to harmonics.
Since the commutation voltage is assumed to be distorted with harmonics,
an iterative solution is required for the crossover points, 3n. A suitable
equation for solution by Newton's method is just the statement that the
commutating voltage is zero at 3
i
dek)e
y^.O.lJ)
where Vdbk and Vdek are the relevant DC voltage samples derived in section
2.8.3 corresponding to the noncommutating states before and after the
commutation, respectively. The zero crossing being solved for S is illustrated
Vdc
Figure 2.21
48
(2.8.16)
(2.8.17)
(2.8.18)
Sample (p)
2
3
4
5
6
7
8
9
10
11
12
DC voltage
Icl
-Id
Id Icl
Id
-Id
Id
el-h
hi
0
I c3
Id
h-h
C
B
Id
*d
Ic4
Id
~Ic4-Id
Id
c5
"Id
-he-
Id
Id
0
I c5
ld
B
,
C
A
.
B
A
C
-Id
A
C
Id
Id
-Id
-Id
b
C
eqn.
2.8.19
2.8.18
2.8.20
2.8.18
2.8.19
2.8.18
2.8.20
2.8.18
2.8.19
2.8.18
2.8.20
2.8.18
Vw-L-V' + L ' ( V r - i t o W - V . - ; t o L , U
(2.8.19)
and
i ^ + y-jfca.w
(2 8 20)
(2.8.21)
50
= 0
*V
(2.8.24)
,x
(2.8.^5)
^ck
*dk
hk - hk ~hk =
z
ck = hk
Vak jiakXak V dk = 0
V^-AA+V^ + V^-iuX^O
(2.8.26)
ak
As for the star-connected source, the solution for the DC-voltage samples is
generalised over all 12 conduction periods into a matrix of coefficients of
the DC and AC sources, i.e.
Vm = VaikVak + PwkV* + PclkVck + VmIA
where
I = 2, 4, 6, 8, 10, 12
(2.8.28)
51
Figure 2.22
[ # P - ap)/2rt
bp>ap
otherwise
(2.8.29)
(2.8.30)
Since the end of one conduction interval is the beginning of the next, all
of the trigonometric evaluations are used in two consecutive sampling
functions, thus halving the number of calculations. The DC voltage can now
be written as
vd=
12
(2.8.31)
52
sampling function p
2K
Figure 2.23
(2.8.32)
Table 2.2
Sample (p)
1
2
3
4
5
6
7
8
9
10
11
12
ap
02
4>2
02
4>2
03
<t>3
05
4>5
06
4>6
04
04
05
05
06
06
01
03
53
-1 cycle-
Figure 2.24
phasors
(2.8.33)
fe=l
1= 0
54
1=1
=k
1=0
(2.8.34)
V.0 = kj I j E O V * ) - 2/(V<Jp0p0)}
P=I
u=o
(2.8.35)
(2.8.36)
and similarly for one of the other two phases. The third phase must always
be the negative sum of the first two, since there is no path for zero sequence
into a bridge. This leads to a total of eight convolutions to calculate the
three-phase currents. As evident in eqn. 2.8.36, the periodic samples for the
phase-current calculation are just the DC-side current, and the commutation
currents derived in section 2.8.1. The calculation of the phase current
flowing into the transformer primary is addressed in the next section. The
derivation of a phase current is illustrated graphically in Figure 2.25.
System side
Unbalance in the tap-changer setting between the two six-pulse groups will
lead to imperfect cancellation of six-pulse harmonics on the AC and DC
sides of the converter. If the impedances in each phase of a three-phase
bank are not all equal, the converter will generate positive and negative-
55
convolved dc current
Figure 2.25
56
Figure 2.26
saturation and hysteresis are not modelled. The resistance and reactance
may be unbalanced, but the tap settings are assumed to be the same on all
phases. The tap-change controller is not modelled as it does not respond to
harmonics. The magnetising current injection is approximated by a shunt
to ground at the primary terminal.
The outcome of this analysis is a transfer model of the transformer; the
primary currents are related directly to the secondary currents and the
secondary voltage to the primary voltage. This is easily achieved for the
star-g/star connection, shown as a single-line diagram in Figure 2.26. The
transformer and thyristor resistances have been referred to an equivalent
AC-side resistance, RAC
(2.8.37)
(2.8.38)
Since all impedance has be^n removed from the transformer, the secondary
voltage is now independent of the current through the transformer
v
(2.8.39)
I1 s
(2.8.40)
These equations are repeated over all harmonics, and all three phases.
57
star-g-delta
primary to secondary
Figure 2.27
58
The J3 scaling for the delta winding does not affect the transfer of thyristor
resistance through the transformer, since it is not connected in delta. Thus,
the referred AC-system resistance is the same as that given by eqn. 2.8.37.
Calculation of the primary-side phase currents in terms of the secondary
currents is complicated by the circulating current in the delta winding. If
the transformer is unbalanced, some of this appears as a positive or
negative-sequence current on the primary side.
The admittance matrix for an unbalanced star-g/delta transformer is
readily obtained
I
where
(2.8.44)
(2.8.45)
(2.8.46)
IPb
IPC
h a
A B
=
[c
D]
(2.8.47)
59
where
a2Ya
A=
c=
a Yfc
0
0
B=
0
2
0
a Y,
-ajSY a
ajSY,
-ap\
a^Yc
-a.pYa
aj8Yfl
(2.8.48)
-a/T
D=
(2.8.49)
(2.8.50)
(2.8.51)
(2.8.52)
(2.8.53)
2.9 Summary
Following an introductory discussion on model accuracy, this Chapter has
described a variety of steady-state AC-DC converter models with varying
degrees of complexity.
The starting case is the fundamental-frequency model based on symmetrical operation, with perfect AC-current filters and DC-voltage smoothing.
This model is used in the power flow solution considered in Chapter 3. The
removal of the symmetry constraint is considered next, and a three-phase
model is developed for use in more detailed power-flow analysis, a subject
60
2.10 References
1 ARRILLAGA, J.: 'High voltage direct current transmission', IEE Power Eng. Ser.,
1988, pp. 246
2 AINSWORTH, J.D.: The phase locked oscillator-a new control system for
controlled static converters', IEEE Trans. Power Appar. Syst., May 1968, 87, (3),
pp. 859-865
3 CHEN, S., WOOD, A. R., and ARRILLAGA, J.: 'HVDC converter transformer
core saturation instability: a frequency domain analysis', IEE Proc, Gener. Transm.
Distrib., January 1996, 143, (1), pp. 75-81
4 WOOD, A.R.: 'An analysis of non-ideal HVdc converter behaviour in the
frequency domain, and a new control proposal'. PhD thesis, University of
Canterbury, New Zealand, 1993
5 GIESNER, D.B., and ARRILLAGA, J.: 'Behaviour of HVDC links under unbalanced ac faults', Proc. IEE, 1972, 119, (2), pp. 209-215
6 SHORE, N.L., ANDERSSON, G., CANELHAS, A, and ASPLUND, G.: 'A
three-pulse model of dc side harmonic flow in hvdc systems', IEEE Trans. Power
Deliv., July 1989, 4, (3), pp. 1945-1954
7 HARKER, B.J.: 'Steady state analysis of integrated ac and dc systems'. PhD thesis,
University of Canterbury, New Zealand, 1980
8 SWARTZ, M., BENNETT, W.R., and STEIN, S.: 'Communication systems and
techniques' (McGraw-Hill, 1966)
9 PERSSON, E.V.: 'Calculation of transfer functions in grid controlled converter
system', IEE Proc, May 1970, 117, (5), pp. 989-997
10 HU, L., and YACAMINI, R.: 'Harmonic transfer through converters and HVdc
links', IEEE Trans. Power Electron., July 1992, 7, (4), pp. 514-525
11 CIGRE SC14, WG 14.25, 'Harmonic cross-modulation in HVDC transmission'.
HVDC colloquium, Johannesburg, South Africa, September 1997
12 SMITH, B.C.: 'A harmonic domain model for the interaction of the HVDC
converter with ac and dc systems'. PhD thesis, University of Canterbury, New
Zealand, 1996
Chapter 3
3.1 Introduction
Power-flow analysis is used to determine the steady-state operating characteristics of the power-generation/transmission system for a given set of
busbar loads. In this context, steady state means a time-independent
condition which implicitly takes into account the final adjustment of the
parameters involved in maintaining the specified operating condition at
each bus, and the component-related capability constraints.
The solution provides information on 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 Newton-Raphson method,1 complemented with sparsity programming and optimally-ordered Gaussian elimination, has become the cornerstone of modern power-flow programmes. The method is described in
Appendix I.
Most of this Chapter describes the operating state of the combined AC
and DC systems under the specified conditions of load, generation and
DC-system control strategies. The superiority of the fast decoupled AC-load
flow2 is now generally accepted and, therefore, the incorporation of HVDC
transmission is described with reference to this algorithm.
In short-term studies, such as the system response to a change in load
specification, the power-flow assessment is made without the assistance of
generator-control adjustment. The term quasisteady state is used here when
referring to such a condition. An assessment of the power-transfer capability
of the DC link also considered in this Chapter comes under this category.
62
Vk - voltage magnitude
6k- voltage phase angle
Two of these variables are specified at the outset and the aim of the power
flow is to find the remaining two variables at each bus.
Three different bus conditions exist:
(i)
63
input
read system data and
node specifications
use system data to form
the system admittance
matrix
initialise voltages at all
system busbars
update voltages and angles
to satisfy the specified
conditions of load and
generation
are
specified
conditions
satisfied
9
the voltage errors less than rj2. Typical figures for r\x and Y\2 are 0.01 p.u.
and 0.001 p.u., respectively.
When the solution is reached, the terminal conditions at each bus are
calculated and with that information the line power flows and system losses
are determined.
lk =
for all k
mek
(3.3.1)
64
where Ik is the current injected into bus k. The power at the k bus is then
given by
S* = P* + jQ* = E*I
= Ek I J&.E*
(3.3.2)
mek
or
P(e,f)
or
QLe,f)
In polar co-ordinates, the real and imaginary parts of eqn. 3.3.2 become
P* = I V,V m (G tm cos 6km + Bkm sin 0km)
(3.3.3)
mek
(3.3.4)
mek
where
AP^lg^+Xj^AV,,,
mek oum
oy
mzk
(3.3.5)
and
^
me* ^ m
V
mek
Oy
(3.3.6)
for a PV busbar
no equations
The voltage magnitudes appearing in eqns. 3.3.5 and 3.3.6 for PV and slack
65
busbars are not variables, but are fixed at their specified values. Similarly,
9 at the slack busbar is fixed.
The set of linearised equations consists of two for each PQ busbar and
one for each PV busbar. The problem variables are V and 9 for each PQ
busbar and 9 for each PV busbar. The basic Newton solution is then
expressed by the equation
P mismatches
for all PQ and
PV busbars
Q mismatches
for all PQ and
PV busbars
AP""1
AQ* *
HP-1
A9P
Np-1
Lp-1
AVP
yp-i
^ 9 corrections
I for all PQ and
J PV busbars
1 V corrections
> for all PQ
J busbars
(3.3.7)
= Vm
and for m k
_ P _ r v2
66
yes
Figure 3.2
67
(3.4.4)
mcok
(3.4.5)
(3.4.6)
u tt = - I u ta
where Zkm and Xkw are the branch impedance and reactance, respectively.
Eqns. 3.4.1 and 3.4.2 can be solved using Newton's method, by expressing
the Jacobian equations as
AP
A0
AQ/V
AV
(3.4.7)
or
(3.4.8)
[AP] = [T][A0]
[AQ/V] = [U][AV]
(3.4.9)
where
[AP] = [AP]
and
[AQ/V] = [AQ]
The decoupled power flow based on the Jacobian matrix equation of the
formal Newton method is
AP
Ad
AV
(3.4.10)
If the submatrices N and J are neglected, since they represent the weak
coupling between P 6 and Q - V, the following decoupled equations
(3.4.11)
[AQ] = [L][AV]
(3.4.12)
The latter equation can be unstable at some distance from the exact solution
and improved convergence is obtained by dividing the right-hand side of
both eqns. 3.4.11 and 3.4.12 by the voltage magnitude V, i.e.
[AP/V] = [A][A0]
(3.4.13)
= [C][AV]
(3.4.14)
These equations are solved successively using the latest values of V and 6
available. [A] and [C] are sparse, nonsymmetric in value and are both
functions of V and d.
The efficacy of the decoupling solution is further improved by making the
following assumptions:
(i) V t > V m =1.0p.u.;
(ii) Gkm BkOT, and hence can be ignored;
(iii) cos(6k-em) = h0
since angle differences across transmission lines are small under normal
loading conditions.
This leads to the following decoupled equations
[AP/V]==[B*][A0]
(3.4.15)
[AQ/V] = [B*][AV]
(3.4.16)
(3.4.18)
no
no
solve eqn. 3.4.17
and update (0)
calculate (AQ/V)
/reactive
powers \
converged
yes
output results
Figure 33
where
Bklfl --
mmk
Bkm = Bkm
m<ak
for m # k
for m # k
69
70
The equations are solved alternatively using the most recent values of V and
Q available as shown in Figure 3.3.6
The matrices B' 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 (between two and 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 AP/V and AQ/V, and are not
susceptible to any humps in the defining functions.
(3.5.1)
(3.5.2)
71
The injected powers Qterm(DC) and Pterm(DC) are functions of the converter AC terminal busbar voltage and of the DC-system variables, i.e.
Pterm(DC) =/(V term) x)
(3.5.3)
(3.5.4)
=0
(3.5.5)
AQ term (V, 6, x)
R(v t e r m , x)k = o
(3.5.6)
=0
(3.5.7)
AQ term (V, 9, x)
R(V term , x)
72
(i)
EUL
a.c. terminal
busbar
QBf
(a 2 -a) Bt
Id
JS
(\-Q)Bf
(b)
Figure 3.4
73
74
(3.6.1)
(3.6.2)
(ii)
(3.6.4)
(iv)
75
(3.6.5)
(3.6.6)
(v)
(vi)
(vii)
(3.6.8)
76
Thus, eqns. 3.6.3, 3.6.4, 3.6.7 and 3.6.8 can be combined into
V^A1.a.Vtcrmcos^ = 0
(3.6.10)
where
The final two independent equations required are derived from the specified control mode.
The DC model may thus be summarised as follows
R(x, V term ) t = 0
(3.6.11)
where
R(l) = Vd-k1.a.VteTm.cos<t>
R(2) = Vd - k , .a .V t e r m .cos + ^ . X ,
(3.6.12)
= Vterni.Ip.cos4>
(3.6.13)
and
or
V (( .I d
(3.6.14)
77
(ii)
sp
a
= 0
Specified DC voltage
V, - Yf = 0
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 AC terminal voltage may also be included as
converter-control equations.
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 freed.
78
(3.6.15)
In this case, the above equation becomes one of the two control
equations and some other variable (e.g. tap ratio) must be specified
instead.
(3.6.16)
This equation is valid for rectification or inversion. However, under inversion the value of Vd, calculated from eqn. 3.6.16, will be negative.
To specify operation with constant extinction angle, the following equation is used
cos (n~y)~
cos (n - ysp) = 0
(3.6.17)
79
Ad '
APterm(V, 6, x)
AQ(V, S)
AQterm(V, 8, x)
A0 t e r m
AV
(3.7.1)
AV t e r m
Ax
R ( V t e r m , *)
(3.7.2)
(3.7.3)
and
P l m n (DC) =/(V t e r m , x)
(3.7.4)
(3.7.5)
AP/V
Ad
B'
DD AA7
AP term /V tern
AQ/V
AV
AQ erm /V tem
BH
AA!'
AV te
BB"
Ax
(3.7.6)
where all matrix elements are zero unless otherwise indicated. The matrices
[B'] and [B"] are the usual single-phase fast decoupled Jacobians and are
kept constant. The remaining matrices vary at each iteration in the solution
process.
80
'term
D D = 0 + - L - {5Pterm(DC)/<5Vterm}
* term
and since
P term (DC) = V(
3P term (DC)/dV lerm = 0
Therefore
DD = 0
Similarly
-^[MPter
"term
i-[5P u
Merm
The powerflowsolution 81
[AA"]=-[8AQterJdx]
"term
1
Mem
BB" = dR/dVteil
[A] = dR/dx
1
5Q,erm(AC)/5Vter
V ter
[aQUrm(DC)/avterm]
= B;;(AC)
In the above formulation the DC variables x are coupled to both the real
and reactive power AC mismatches. However, eqn. 3.7.6 may be separated
to enable a block-successive iteration scheme to be used.
The DC mismatches and variables can be appended to the two fast
decoupled AC equations in which case the following two equations result
AP/V
APtCrm/Vterm
R
AQ/V
AQ term /V ter ,
R
Ad
B'
AA'
A
Ax
(3.7.7)
"A?"
B"
BB"
AA"
AVterm
Ax
(3.7.8)
82
IT
Figure 3.5
83
These features justify the removal of the DC equations from eqn. 3.7.7 to
yield a -P,QDC- block successive iteration scheme represented by the
following two equations
(3.7.9)
[AP/V] = [B'][A0]
AQ/V
AQ
IV
term'
term
AV
B'
AA"
AVterm
BB
Ax
(3.7.10)
AQ/V
AV
AA';t
AQ t e r m i/V t e r m l
AV t e r m l
AV t e r m 2
AQterm2/ M e r m 2
AQ,e3/V ter m3
AA'33
AQ,erm4/Vter m4
AA44
^MOD
A Q t e r m 5 / V t e r m5
AA'5 5
AR2
AR3
AR 4
AR5
AR6
AV t e r m 4
AV t e r m 5
AV t e r m 6
AQ t e r m 6/V t e r m6
ARX
AV t e r m 3
Bn
AXl
K2
AX2
A
(30 x 30)
B33
n"
B55
B'e 6
AX3
AX 4
AX5
AX6
(3.7.11)
84
Figure 3.6
where B^OD iS ^ e part of B'' which becomes modified. Only the diagonal
elements become modified by the presence of the converters.
Off-diagonal elements will be present in B^OD *f there is any AC connection between converter terminal busbars. All off-diagonal elements of BB"
and AA" are zero.
In this example, matrix A is block diagonal in 5 x 5 blocks with the
exception of the DC interconnection.
The DC-side configuration is used to derive the set of DC equations,
corresponding with the term R(3) of the point-to-point interconnection
(eqn. 3.6.11). For the test system of Figure 3.6 these are
- I<)3 R<(3 = 0
hi ~ hi = 0
(3.7.12)
85
v term J
(3.7.13)
where
is the mismatch calculated in the absence of the DC converter, and B" is the
constant AC fast decoupled Jacobian.
After triangulation down to but excluding the busbars to which DC
converters are attached, eqn. 3.7.13 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.
The DC-converter equations may then be combined with eqn. 3.7.14 as
follows
(AQ/V)"
ermy
V
T
, AQ terro (DC)
V
ter
AV
0
B"'
AVtern
B;r + B;;(DC)
AA"
BB"
Ax
(3.7.15)
where
'ter
86
AQ term Y , AQterm(DC)"
vten
Merir
B;;' + B;;(DC)
AA"
BB"
AV tern
(3.7.16)
may then be solved by any method suitable for nonsymmetric matrices.
The values of Ax and AVterm are obtained from this equation and AVterm
is then used in a back substitution process for the remaining AV to be
completed, i.e. eqn. 3.7.14 is solved for AV.
The most efficient technique for solving eqn. 3.7.16 depends on the
number of converters. For six converters or more the use of sparsity storage
and solution techniques is justified; otherwise, all elements should be stored.
The method suggested here is a modified form of Gaussian elimination
where all elements are stored but only nonzero elements processed.8
am P d m
am * dm
am * dm
a m Prfm
am ?dm
am Pdm
am h
vdm
7n
an
an
dn
dn
In
In
n
7n
In
' dn
dn
The DC-link data and specified controls for case 1 are given in Table 3.1.
Initial values for the DC variables x are assigned from estimates for the
DC power and DC voltage and assuming a power factor of 0.9 at the
converter-terminal busbar. The initial terminal busbar voltage is set at
1.0 p.u.
87
Characteristics of DC link
AC busbar
DC voltage base
Transformer reactance
Commutation reactance
Filter admittance B*
DC link resistance
Converter 1
Converter 2
bus 5
100kV
0.1260
0.1260
0.4780
bus 4
100kV
0.0728
0.0728
0.6290
0.3340 ohms
58.6 MW
7
10
-128.87 kV
3.9 Modification of the power flow for use with the unit
connection
The use of direct-connected generators to HVDC converters is under
serious consideration.10 The two alternatives, termed unit and group
88
Figure 3.7
Unit-connected HVDC power station, in which generators and converters are integrated into single units, series-parallel and/or parallel
combination of units at the DC side is also allowed
connection and illustrated in Figures 3.7 and 3.8, respectively, are equivalent from the power-flow modelling viewpoint.
In the absence of filters, the voltage at the converter terminal is not
sinusoidal and cannot be used as the commutating voltage in the conventional formulation. Instead, the use of the generator internal e.m.f. behind
subtransient reactance (E") is suitable.
The fast-decoupled power-flow formulation described in previous sections
should be applicable by shifting the converter interface to a fictitious
terminal (the internal e.m.f. behind subtransient reactance) where the
waveform is assumed sinusoidal.
Although E" is neither accessible nor directly controllable and varies with
the load, its magnitude can be derived from the generator's subtransient
Y Y
Figure 3.8 Group-connected HVDC power station, in which machines and converters can be combined in groups via a transfer bus
89
xlc
Figure 3.9
phasor diagram for the nominal operating point and then kept constant (as
shown in Figure 3.9). This is done to satisfy the generator's conventional
power-flow specification (i.e. constant voltage) but at the fictitious internal
E" bus. The conventional power flow is then carried out on the assumption
of perfect filtering at that point.
However, the internal e.m.f. behind subtransient reactance is not directly
controllable. Instead, the generator excitation should be controlled to
provide the specified firing angle (amin).
Therefore, the commutating voltage (E") will not be shown in advance
and its magnitude and phase angle (/?)> as shown in Figure 3.9, must be
derived as part of the iterative solution.
The five-variable model of section 3.6 at the rectifier end is replaced by
a set of seven variables
(3.9.1)
Therefore, seven residual equations are needed to formulate the power flow
problem at the rectifier end. 11
The first two equations are common to the conventional model (with the
tap-ratio variable removed by making ar = 1), i.e.
Vd - k^'cosot - ~xcld = 0
n
where
3V2
(3.9.2)
90
Two more equations are derived from the generator subtransient phasor
diagram of Figure 3.9. That is
E - E"cos 0 - (x - x")\Ip\sin(fi + <j>) = 0
E" sin P-(x-
(3.9.3)
(3.9.4)
(3.9.5)
E- W = 0
(3.9.6)
(3.9.7)
(3.9.9)
which contains the five variables of the conventional (inverter) set and an
extra variable (cos a) representing the rectifier end of the link.
91
the specified current level derived from a constant-power setting (P^p) at the
inverter end.
Besides Psf and ym, a third control specification must be made to match
the six-variable formulation. It is suggested that am is used and the inverter
transformer tap position made free. This will provide the highest transmission voltage. Should the tap changer attempt to violate one of its limits,
this limit becomes the new control specification while freeing the value of a.
The vector of independent variables, x, which describes the state of the
DC link at the inverter end, can still be obtained by the application of the
Newton-Raphson algorithm. The set of residual equations, |R|, to be
satisfied upon convergence is
(3.9.10)
3
Vd - &iaVterm cos {n - y) - -xcld = 0
(3.9.11)
(3.9.12)
(3.9.13)
(3.9.14)
(3.9.15)
The incorporation of the equivalent inverter into the fast-decoupled ACDC power-flow algorithm requires modification of the B", AA', AA" and BB"
submatrices of the Jacobian matrix in eqns. 3.7.7 and 3.7.8. These submatrices are greatly simplified in the equivalent inverter model as they contain
only a single element each. The new elements are
[AA'] = [<5Pterm(DC)/&]/Vterm
or
[AA']T =
[AA"] = [<5Qterm(DC)/<5x]/Vtern
92
or
[AA"] T =
_2a i B i V t e r m i
[BB"] = <5R/<5Vtern
or
[BB"]T =
Figure 3.10
real power
'practical
stability
limit
93
Benmore
16kV
Figure 3.11
Haywards
110kV
94
reactance of 2% (on the system MVA base). Since, for the purposes of the
model, it is fed by a sinusoidal voltage source, the transformer reactance can
also be regarded as the commutating reactance of the AC-DC converter.
The converter is modelled as a three-phase 12-pulse bridge rated at 1.2 kA
DC-side current and 525 kV DC-side voltage. The converter station actually
consists of four six-pulse bridges in series.
An overhead transmission line and an undersea cable transmit the power
to an identical converter at Haywards. The large smoothing reactors at each
end of the line are assumed to maintain a constant DC-side current. The
total resistance of the converters, smoothing reactors, transmission line and
cable is 25.56 Q.
Another 750 MVA transformer connects the Haywards converter to the
HOkV Haywards busbar. As at Benmore, the converter transformer reactance of 2% is regarded as being the commutating reactance because a bank
of harmonic filters maintains a sinusoidal voltage on the HOkV busbar.
The voltage level of the HOkV busbar is maintained by the voltage
regulators on synchronous compensators which can provide up to 260 MVAr
of reactive power. This is supplemented by 110 MVAr from the harmonic
filters. The busbar feeds the North Island AC network via a transmission
line.
Similarly to the capability chart of a generator, which represents the
power available from the machine terminals, the most useful information to
be derived from an HVDC-link capability chart is the power received at the
inverter end. In the test system, this is the power available to the Haywards
110 kV busbar from the converter transformer.
L
M
N
O
P
Q
R
S
95
(3.10.1)
~~ M>C*D
(3.10.2)
(3.10.3)
(3.10.4)
VDC =
(3.10.5)
(3.10.6)
P+jQ
ysu
A
Xc
zx
h
dc
dc
96
Eqn. 3.10.6 is used to draw the loci for maximum and minimum valve
currents. These are circles centred on the origin with radii of Vp#K(3^/2/
^Iocmax and VpaK(3V2/7t)IDCmin.
The maximum r.m.s. converter transformer current It is equal to the
per-unit MVA rating of the transformer. If the effects of commutation are
ignored, the converter transformer current can be related to the DC-side
current in per unit by
The DC-side current from eqn. 3.10.7 is substituted into eqn. 3.10.6 to
produce a circular locus centred on the origin with radius VpaK(3^j2/n)It/
yj2 describing the maximum converter transformer current. If the effects of
commutation are considered, the locus would be an approximate circle with
a slightly larger radius. Consequently, there is a small margin of safety
inherent in the circular locus described here.
The loci representing maximum DC-side voltage VDCmax are obtained by
combining eqns. 3.10.1, 3.10.2 and 3.10.3 to give
Q = - |P| Vt(V,aK(3%/2/7r)^rDCmax)2 - 1]
(3-10.8)
This is the equation of two straight lines extending from the origin of the
complex power plane. The gradients of the lines are
V[(VpaK(3 V2/7t)/VDCmax)2 - 1]
Two parametric equations are used to describe the locus representing the
minimum firing angle, amin. The variable parameter used is the DC-side
current 1 ^ . The real power co-ordinate of the locus is given by substituting
eqn. 3.10.4 into eqn. 3.10.1 to give
P = - [(3V2/7t)Vpcos(amin)IDC - (3/7r)XcIjy
(3.10.9)
The reactive power co-ordinate is given by substituting eqn. 3.10.3 into eqn.
3.10.2 and using the real power from eqn. 3.10.9 to give
Q = - V[V> 2 K 2 (3 V2M) 2 I^ - P2]
(3.10.10)
The parametric equations are used to draw the locus by gradually increasing
the variable parameter IDC starting from zero. For each new value of IDC, a
corresponding locus-point co-ordinate is given by eqns. 3.10.9 and
3.10.10.
The locus representing minimum extinction angle, dmin, is also parametrically described and is obtained by changing the sign of the power flow in
97
(3.10.11)
(3.10.12)
This value of current is then substituted into eqns. 3.10.1 and 3.10.4 to
give the real-power co-ordinate
P = - [(3^2/7c)aVpcosaIDC - (3/TC)XCI^C]
(3.10.13)
(3.10.14)
Qr2)
(3.10.15)
2
- P ]
(3.10.16)
(3.10.17)
98
= VrI?*
max
-1000
real power, MW
Figure 3.13
1000
99
This combination of converter constraints restricts the Haywards converter to operating within the area bounded by loci D-F-H-C-E of Figure
3.13. However, it is not sufficient to consider only the constraints of the
Haywards converter. The Benmore converter and Benmore AC system must
also affect the power available from the link. The operating constraints of
the Benmore circuit can easily be referred to the Haywards end of the link
by using the DC-link power-transfer mapping. This is a point-to-point
mapping which relates the complex power entering the link from the
Benmore 16kV busbar to the complex power leaving the link at the
Haywards HOkV busbar. The mapping is carried out using eqns. 3.10.15
and 3.10.16.
The constraint-locus equations of the Benmore converter can now be first
formulated in the same manner as those of the Haywards converter. This
produces a set of Benmore-converter loci which refer to the power entering
the link at the Benmore 16kV busbar. These loci are identical to the
Haywards-converter loci shown in Figure 3.13. The DC-link power-transfer
mapping is then used to transform these loci to produce a new set of
Benmore-converter constraint loci which are referred to the Haywards
110 kV busbar.
The transformed Benmore-converter loci are shown in Figure 3.14. All
the constraints involving current flow are still represented by arcs of circles.
However, the shapes of all loci have been altered by the transformation.
The complex power supplied from the South Island 220 kV busbar to the
Benmore 16kV busbar is dependent on the voltage angle between the
busbars. This power can be represented by a circular locus. Locus P in
Figure 3.15 is a portion of this circle. The power is restricted by the current
rating of the 400 MVA interconnecting transformer. The circular locus O
represents maximum current flow in the transformer. The intersection of
loci P and Q therefore represents the maximum possible power delivery to
Benmore. The arc of locus P that lies within locus Q represents the full
range of power delivery.
-1000
real power, MW
-1000L
reactive power, MVAr
1000
100
-500L
Figure 3.15
50
Q
real power. MW
1000
Q
1000
5
-500
Figure 3.16
-1000
101
1000
-1000L
reactive power, MVAr
Figure 3.17
posing Figure 3.13 onto Figure 3.14. This chart describes the power that
may be delivered to the North Island from the link assuming that an
infinitely strong 16kV busbar exists at Benmore. The Benmore currentconstraint loci overlay the corresponding Haywards current-constraint loci
because the two converters are identical.
Manipulation of the control angles at Benmore and Haywards allows the
operating point of the link to vary within the area bounded by loci
KD-L-H-CJ-O-E. When maximum power is being delivered to the North
Island, the operating point is positioned by the intersection of loci CJ and
H. This closely matches the practical operation point.
If the Benmore generators are run as synchronous compensators and the
power-generation capability of the South Island is considered, the loci of
Figure 3.15 must be incorporated into the capability chart. This is achieved
by applying the DC-link power-transfer mapping to the loci of Figure 3.15
and superimposing the transformed loci onto Figure 3.17.
real power, MW
-1000
1000
iliillll
y0
\
CJ
XX
H
_mnr>
reactive power, MVAr
Figure 3.18
102
ESCR
SCR<
Pd
-GD-S
Xc
Id
Qc
UL
Ud
103
1.8
.^
1.6
1.4
1.2 m
\
Xc = 0.15 pu.y= 18
2.0(1.46) N B , Qc = Qd = 0.54 Pdn at UL = 1.0 pu
^
\
DC Power
^
^
AC Volts
3.0(2.46ir\
^ 1.0
o
0.8
3.0
0.6
0.4
0.2
u
1.0 1.2
1.4
1.6
1.8 2.0
DC Current (pu)
Figure 3.20
104
0.2
0.4
0.6
0.8
1.0
1.2
105
1.4
Xc = 0.15pu,Y=18
Qc = Qd = 0.54Pdn at U L = 1.0 pu
1.6
1.4
SCR = 2.0
MAP-1
1.2
B
Pd DC Power
UL - - AC Volts
1.0
o
o 0.8
MPC-1
SCR = 3.0
CL
0.6
MPC - 2 \
SCR = 2.0
0.4
0.2
0.2
0.4
0.6
0.8
1.0
1.2
1.4
1.6
1.8
2.0
DC Current (pu)
Figure 3.22 AC-DC systemlow SCR power and AC-voltage curves sudden
change of SCR from 3 to 2
CIGRE, 1992 Reproduced by permission
106
The required power cannot be maintained but the current limit imposed
by the rectifier prevents the collapse of the system despite the demand for
higher power from the master controller.
The operating point B on MPC-2 will be in the region where dd/dld
is negative and stable operation would not be achieved in power-control
mode but stable operation continues in constant-current control mode.
If the DC system is used to control AC frequency or to provide system
damping, operation in the constant-current mode would not be acceptable
except for a very brief period of time; e.g. while the fault is being cleared.
ESCR
K
a
e.m.f.
Example data
Rectifier
Inverter
2.75
15%
12
1.05
1.00
14%
25
1.17
107
Case
Power
p.u.
AC voltage on the
remote busbar
1
2
.970
1
.990
1
udc
.970
1
ld
a
rectifier
7
inverter
1
1
3
12
26.7
18.4
case). For example, an increase in current limit by 0.05 p.u. for case 1 would
cause the power to be maintained at 1 p.u. (instead of 0.97 p.u.) but with a
larger reduction of remote AC voltage, by 3.2 per cent instead of one per
cent. In other cases an increase in current limit would actually reduce power,
with a further voltage depression, hence would be pointless.
In the final steady state, the capacitor switching and tap-changer control
would restore the specified conditions.
3.12 Summary
Using the basic converter steady-state formulation developed in Chapter 2
for power frequency and symmetrical operation, this Chapter has described
the incorporation of AC-DC converters, point-to-point DC links and
multi-invertor schemes in the conventional power-flow solution.
The overall convergence rate of the AC-DC algorithms depends on the
successful interaction of the two distinct parts. The AC-system equations are
solved using the well-behaved constant tangent fast-decoupled algorithm,
whereas the DC-system equations are solved using the more powerful, but
somewhat more erratic, full Newton-Raphson approach.
In general, those schemes which acknowledge the fact that the DC
variables are strongly related to the terminal voltage give the fastest and
most reliable performance.
In practice, converter operation has been considered down to an SCR of
3. A survey of existing 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 power flow.
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 power-flow sense.
An alternative model has been proposed for the integration of unitconnected schemes.
Constraint loci have been developed to take into account the limited
capability of the converter plant. The question of system-related DC-link
108
3.13 References
1 VAN NESS, J.E., and GRIFFIN, J.H.: 'Elimination methods for load-flow studies',
Trans. AIEE, 1961, PAS-80, pp. 299-304
2 STOTT, B., and ALSAC, O.: 'Fast decoupled load flow', IEEE Trans., 1974,
PAS-91, pp. 1955-1059
3 BROWN, H.E., CARTER, G.K., HAPP, H.H., and PERSON, C.E.: Tower flow
solution by impedance matrix iterative method', IEEE Trans.,, 1963, PAS-82, pp.
1-10
4 TINNEY, W.F., and HART, C.E.: Tower flow solution by Newtoh's method',
IEEE Trans., 1967, PAS-86, (11), pp. 1449-60
5 STOTT, B.: 'Decoupled Newton load flow', IEEE Trans., 1972, PAS-91, pp.
1955-1959
6 TINNEY, W.F., and WALKER, J.W.: 'Direct solutions of sparse network equations by optimally ordered triangular factorization', Proc. IEEE, 1967, 55, (11),
pp. 1801-1809
7 ARRILLAGA, J., ARNOLD, C.P., and HARKER, B.J.: 'Computer modelling of
electrical power systems' (John Wiley & Sons, 1983)
8 BODGER, B.P.: 'Fast decoupled ac and ac-dc load flows'. PhD thesis, University
of Canterbury, New Zealand, 1977
9 IEEE Computer Applications Sub-Committee Standard Test System, American
Electric Power Service Corporation, 1962
10 CAMPOS BARROS, J.P., ARRILLAGA, J., BOWLES, J., KANNGIESSER, K.W.,
and INGRAM, I.: 'Direct connection of generators to HVDC converters: main
characteristics and comparative advantages', ELECTRA, August 1993, 149, pp.
19_39
11 ARILLAGA, J., ARNOLD, C.P., CAMACHO, J.R., and SANKAR, S.: 'AC-DC
load flow with unit-connected generator-converter infeeds', IEEE Trans., May
1993, PS-8, pp. 701-6
12 DE SILVA, R., ARNOLD, C.P., and ARRILLAGA, J.: 'Capability charts for an
HVDC link', Proc. IEE, May 1987, 134, (3), pp. 181-6
13 CIGRE Working Group 14.07: 'Guide for planning dc links terminating at ac
systems locations having low short-circuit capacities. Part I: AC/DC interaction
phenomena'. Report 68, June 1992
Chapter 4
4.1 Introduction
Chapter 2 (section 2.8) has described a model of the AC-DC converter in
the harmonic domain which, for a given set of switching instants, converter
terminal voltage and DC-current waveforms, can calculate the returned
currents.
currents.
In practice, however, the terminal-voltage waveforms are part ot the
overall solution which involves not only converter variables, such as controller characteristics and firing-angle constraints, but also the system voltage
sources, current sources and impedances. Iterative techniques are thus
necessary to solve all these variables together to reach a final correct
solution.
. .
.
.
The simplest iterative scheme used to simulate harmonic interactions is
based on fixed-point iterations. At each iteration the latest values of the
distorted converter terminal voltages are used as the commutating voltages
for the converter solution.
The calculated DC-voltage waveform is then impressed upon the DC-side
impedance to derive the DC-current waveform, which in turn, together with
the switching instants and commutation process, provides the AC-side
current harmonic injections. The latter are then used to derive the ACvoltage harmonics in the frequency domain for the next iteration.
However, the fixed-point iteration solution diverges when the AC-system
harmonic impedance is large and the commutating reactance small. This
Chapter describes the use of the Newton method to model the interaction
between the converter and the AC and DC systems.
110
J J
\
controller
YgD
uc
alpha order
\\ filter
VD
Figure 4.1
l'-a
111
Figure 4.1. The starting point will be a relationship between the converter
phase currents and the converter terminal voltage.
Given the primary phase currents, the terminal voltage can be found by
(4.3.1)
m * = [ Y j ^ - t i p + i & + [Ytc]k\yth]k)
V=f1(I%I$
(4.3.3)
= [V]k - [R ][I^
(4.3.4)
(4.3.5)
V? =/ 3 (V)
(4.3.6)
tf)
tf,
tf)
(4.3.7)
(4.3.8)
112
(4.3.9)
(4.3.10)
(4-3.12)
(4.3.14)
The transformer analysis of section 2.8.4 then describes how the primary
currents are obtained
IP=/IO(IS)
I?=/u(V?,l)
(4.3.15)
(4.3.16)
K^fufyUa, 0?,$)
(4.3.17)
1%=fi3(yZ,hk,tf,4>?)
*Wi4fl.I>. 0?>o)
tf=fis(hk>ho,0?>*o)
F s =/ 16 (V d0( I d0 )
(4-3.18)
(4.3.19)
(4.3.20)
(4.3.21)
Fo=/i7(VdO.I,/o)
(4.3.22)
113
vPD
Vdk
OK
c
CO CO
I*
's
Ip
Ip
F^
?
r
e
F?
Fs
Fa,
No. of vars
300
300
300
100
1
100
1
300
300
300
300
6
6
6
6
1
1
Function of
'i(lp,lp)
400 '
400s
f (\i
\#D 1 jk>
r\S r\O
iS
Y^/
JD\
Wi )
f6(Vdk)
^(vdo)
fe(V% [dk, \6
f g (Vp, \dk, \0f0,
f10(\l)
fu (Vp, IQ)
f12(Vp, 1 ^ t?f,
'\4\'dk>
'd0' ^
f,5(\dk, l d 0 , t >?,
0f, 0,D)
0f)
oc0)
f^dO, ld0)
, ho, 01 <^)),
114
how the DC current and terminal voltage, together with the switching
angles, are used to calculate the primary phase currents. The primary phase
currents are then injected into the AC-system impedance to yield the
terminal voltage. The relationship
V =/ 18 (V, I*, I, o , 9f, 4> $, 0?)
(4.4.2)
(4.4.3)
=f19(y,idk,o?,<p?,es,ti)
(4.4.4)
(4.4.5)
Eqn. 4.4.5 yields 100 equations when decomposed into harmonic real and
imaginary components and is a composition of functions which describe the
calculation of the DC voltage and its application to the DC-system model,
to yield the DC current. A further 26 mismatch equations are obtained in a
similar manner, related to the converter controller and the switching
instants
=fzo<y,l, $.<!>?)
(4-4.6)
(4-4.7)
( 4 - 4 -8)
(4-4.9)
115
(4-4.10)
V?, Idk, 0?, 0?, $ , #
3(V),
I*, fl?,
tf)
(4-4.11)
(4.4.12)
rv irv 1 f (f (\r\ i
#s rfi\
f (f (\r\ T
PP
(4.4.13)
Table 4.2 Mismatches and variables for the 12-pulse rectifier
Variable
No. of vars
Fv
300
lld
100
Function of
y v . i* u fl?, 0, ^f, 0)
o)
116
Note that the current mismatch is still expressed in terms of the same
variables as the voltage mismatch. The current mismatch has the advantage
that it does not require the system admittance to be inverted, a possible
difficulty if it has a high condition number. The current mismatch is also
the preferred method of modelling the interaction with a purely inductive
AC system, such as the unit connection, as the system admittance will
decrease with increasing harmonic order. Only the voltage mismatch will be
implemented here, as the AC-system admittance is usually invertible, and
the AC-system impedance is typically much less than one per unit.
The DC-current mismatch, F w , defining the interaction with the DC
system, can also be written as a DC-voltage mismatch, F yd . This mismatch
is obtained by injecting the estimated DC current into the DC-system
impedance and comparing the resulting voltage with the calculated DC
voltage
(4.4.14)
(4.5.1)
(4.5.2)
(4.5.3)
117
i) =
(4.5.4)
-Y
(4.5.5)
with convergence deemed to have occurred when some norm of the residual
vector F(XN) is less than a preset tolerance. Newton's method is not
guaranteed to converge, but convergence is likely if the starting point is
close to the solution.
SiJ
AF,AX,.
F(* + A * ) - F ( x )
Ax
5bJ
Provided that AXj is small enough, this gives a good approximation to the
Jacobian.
The numerically-derived Jacobian for the test system has been plotted in
Figure 4.3, for a solution up to the 13th harmonic. This Jacobian was
calculated at the solution, and so represents a linearisation of the system of
equations in Table 4.2 around the converter operating point.
Referring to Figure 4.3, the elements of the Jacobian have been ordered
in blocks corresponding to the three phases of terminal voltage, the DC
current, the end-of-commutation angles, the firing angles and the average
delay angle. The blocks associated with interactions between the DC-current
harmonics and the AC-voltage harmonics comprise the AC-DC
118
--
co
^T
c2
a.*
S~2
s
fid
CO
-U
""!
J5
>
A.*
co*
fa
HH3
CO
s
s
-S 2
d fit
fit
c^ g
s-2
CO
v-'
s
c2
fa^ N J fa^
co
^ fit
NJ
CO
CO
*m*r
CO
CO
CO
fe
cS
co
1
J
^ t fa
od ^T
co
CO
c^
% &
CO
HH3
cS
fa
fa
N? k j
CO
co
s
co
CO
CO
CO
I co
'fid
c^
^ r ftj
co I co
rig 3|S :
'Is S'
S s"
9 fe*
CO
8 28 ^ Es* 83 V j
w
& r, fa
CO
h-< fit
c^
p P 2 P P ^ p 9 P fit
si s s" s*i s s
s" S S s"
co
8 2P JP 3 P 5
P JS
CO
CO
2
s
CO
co
CO
82
8 2p J
2P
SN
s" s"
8 28 2
119
Numerical Jacobian
mismatches
Figure 43
variables
partition, which is 104 elements square. All other parts of the Jacobian are
called the switching terms. Within the AC-DC partition, elements have
been arranged within each block in ascending harmonic order, with the real
and imaginary parts of each harmonic alternating. Each block in the
AC-DC partition is, therefore, 26 elements square in Figure 4.3 but 100
elements square for a solution to the 50th harmonic.
The Jacobian displays several important structural features:
A The test system contains a parallel resonance in the AC system at the
second harmonic. This leads to rows of large terms in the Jacobian
aligned with the second-harmonic terminal-voltage mismatch (the resonance terms).
B A change in harmonic k on one side of the converter affects harmonics
k + 1 and k 1 on the other side of the converter, causing the double
diagonal structures in the AC to DC and DC to AC blocks. These are
the three-port terms.4
C The end-of-commutation mismatch is very sensitive to harmonics in the
terminal voltage and DC current.
D The individual firing instants are sensitive to harmonics in the DC
current.
120
AC-DC
Va
Vb
Id
100
120
20
40
60
80
100
I,
Figure 4.4
120
Va
Vb
121
Vc
Id
120
r
20
40
60
80
100
120
Get
Figure 4.5
122
as that required to calculate the complete set of mismatches just once. For
the converter system of 426 mismatch equations, the analytic derivation of
the Jacobian is 20 times faster than the numerical derivation. A full
derivation of the Jacobian elements can be found in Reference 5.
As an illustration of the analytical derivation, the elements corresponding
to the end-of-commutation mismatch equation are described here.
The current in the phase that is commutating off at the end of the
commutation should be zero, and is given by
= 0
(4.5.7)
jkcoL*eb
t
(4.5.10)
variations
dR{V3m}
dR{V6m}
dI{Vdm}
dl{\dm}
(4.5.11)
and
dI{Vdm}
(4.5.12)
The partial derivatives in these equations are obtained from eqns. 4.5.8 to
4.5.10 as
(4.5.13)
(4.5.14)
dI{V*m}
123
ejm
(4.5.15)
(4.5.16)
VIX
\Kdmi
UMX K
\ dml
urL L
\ dm}
(4.5.17)
)
and
The partial derivatives in these equations are obtained from eqns. 4.5.8 to
4.5.10 as
(4.5.19)
L. + L,
Jk0i
dl{ldm}
Le + I ,
T = v^h-
(4.5.20)
(4.5.21)
(4-5.22)
This analysis assumes that the commutation is on the positive rail. A similar
analysis holds for a commutation on the negative rail but with \dm
substituted into eqn. 4.5.7.
124
(4-5.23)
J
125
Solving power-system elements in sequence components affords a common frame of reference within which to model the system as a whole. In
such a case, it is necessary to integrate machine models, the load flow,
HVDC, FACTS and transmission-system components into a single iterative
solution. In addition to being more computationally efficient, the sequencecomponents frame of reference leads to greater insight into the interactions
between component nonlinearities, communications interference and power
quality.
In the existing model, interaction with the DC system is specified by
summing the DC-voltage harmonics across each bridge and applying the
resulting voltage harmonics to the DC-system linear model. The resulting
harmonic current
(4.5.25)
(4.5.26)
(4.5.27)
(4.5.28)
A sequence-components solution is implemented by using the sequencetransform matrix, T, to interface complex phase-component calculations to
the Newton solution in sequence components. For complex phasors, the
complex 3 x 3 sequence transform is
1 a2
1
JV
a*
(4.5.29)
126
v
V"
as the sequence components yields W = TV and V = T 1W. If a threephase quantity has been decomposed in real rectangular components, as in
the Jacobian matrix, a real components-sequence transform matrix, T, can
be constructed
W = TV
which in full is
"l
v? 0
VR" 1
V," 0
V
1
R
+ 0
v,
v?
Ra
-la
VR
0 Ra
-la
VR
la
Ra
la
Ra
vf
Ra
-la
Ra2
-la2
VR
la
Ra
la
Ra
(4.5.30)
VJ
where R, I denote the real and imaginary parts of a quantity. The existing
phase-components mismatch equations summarised in Table 4.2 can now
be written as sequence-component mismatches, in terms of the sequencecomponents terminal voltage, using the sequence-transform matrix. Sequence-component mismatches for the 12-pulse controlled rectifier are
listed in Table 4.3.
Table 4.3 Mismatches and variables for the 12-pulse rectifier with the terminal
voltage in sequence components
Eqn.
\dk
Y
F?
r
Functional notation
127
128
Table 4.4 Convergence and performance of the converter model, (b) constant unified
Jacobian, (d) constant sequence components Jacobian
Test no.
Main iterations
1fc
11.1
10.9
14.8
15.0
12.3
11.5
6
5
11
12
7
6
0.1712
0.0851
0.2562
0.2973
0.1130
0.1056
Id
2b
2d
3b
3d
The resulting Jacobian matrix is plotted at the bottom of Figure 4.6, and
shows a greater degree of sparsity than that for the phase components,
shown in Figure 4.3. This is primarily due to the absence of any coupling
between the zero-sequence AC-voltage harmonics and any of the other
converter variables.
Convergence of the sequence-components solution is compared with that
of the constant Jacobian phase-components solution in Table 4.4, for three
of the test cases. They are very similar, indicating that the converter model
can be directly combined with other components (such as the synchronous
machine) modelled in sequence components.
129
update switching
system variables
(a)
130
4.6.1 Initialisation
An initial estimate for the converter delay angle is obtained from the
equation
V,o = |V ( f c l |cosa- I, o
n
n
(4.6.1)
(4.6.3)
These angles are then used to assemble the individual firing and end-ofcommutation angles
6t = pi + a
(4.6.4)
& = 0 | + A
(4.6.5)
FaO(0?,6l?,V,I<l) = O
131
(4.6.6)
A flow diagram for the switching system is shown in Figure 4.7b, and the
structure of the switching Jacobian can be seen in Figure 4.8. Those
elements corresponding to the partial derivatives of voltage mismatch with
respect to end-of-commutation angle have been set equal to zero, since they
are always insignificant.
'dO
5-
10-
520-
25--
(
30--
10
Figure 4.8
15
20
Variable number
25
30
132
l < 0.001
2} < o.ooi
15ff!
0.001
(4.6.7)
133
< 0.001
<
Haoi
< 0.001
J ^ U 0.001
|Fgi|<5 x 10~ 8
|F^| < 5 x 10" 8
(4.6.8)
134
to a tolerance of |X|X < 10 5 is suitable for general purpose use. This type
of convergence test is fast and easy to apply but does not imply that all
harmonics have converged to a satisfactory accuracy.
135
EMTDC
n. Harmonic Domain
1O
16
2O
2S
-4O
3O
35
3O
35
AO
3O
35
AO
3O
35
-*O
15
2O
25
harmonic order
Phase A current, Test One
53 O . 1 2
EMTDC
3*
O.1
O.O8
3 O.O6
Harmonic Domain
2) O.O4
CO
0.02
1O
15
2O
harmonic order
25
1O
Figure 4.9
15
2O
25
The results of test one, the base case, are shown in Figures 4.9 and 4.10. It
can be seen from the spectra that only the characteristic harmonics are
present and that there is a close match between the time and harmonicdomain solutions. The time-domain solution yielded small residual noncharacteristic harmonics in the spectra, which were suppressed in the phase
graphs. The DC-voltage waveform was generated from the harmonicdomain solution by plotting the Fourier series for each DC-voltage sample
during the appropriate interval, rather than inverse transforming the
136
time (sec)
0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018 0.02
time (sec)
137
0.01
0.02
time (sec)
400
Figure 4.11
0.02
harmonics are present on the DC side and even harmonics on the AC side.
In particular, the fundamental resonance on the DC side is excited and
there is a large fundamental component in the DC voltage.
Unbalancing the star-g/delta transformer leakage reactance (test 3),
shown in Figures 4.13 and 4.14, causes the generation of many harmonics
and a large zero sequence on the AC primary side. This is due to sequence
transformation by the transformer, irrespective of its connection to a
converter. The zero sequence, also plotted in Figure 4.13, shows a close
agreement between the harmonic and time-domain solutions.
The main purpose of test 4 was to show that convergence is still acceptable
when the resistance present in the commutation circuit is not represented
in the Jacobian matrix. Also, setting the tap changer on one transformer but
not the other introduces six-pulse unbalance. It was necessary to increase
138
O.3
O.2
EMTDC
Harmonic Domain
O.1
15
2O
25
3O
35
4O
harmonic order
DC voltage, Test Two
1O
15
2O
26
3O
35
4O
3O
35
4O
3O
35
4O
harmonic order
Phase A current, Test Two
SS O.12 -
EMTDC
Harmonic Domain
O.O2
1O
15
2O
harmonic order
25
nil! !L
25
Figure 4.12
139
0.02
400
0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018 0.02
time (sec)
Zero sequence ac current, Test Three
I
-0.1
-0.2
Figure 4.13
0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018 0.02
time (sec)
140
EMTDC
ih Harmonic Domain
1O
15
2O
25
3O
35
3O
35
4O
harmonic order
DC voltage, Test Three
15
2O
25
4O
harmonic order
Phase A current, Test Three
d
d.
0.12
EMTDC
0.1
Harmonic Domain
0.08
o.oe
O.O4
0.02
1O
15
2O
IH
25
m.
harmonic order
3O
Bl
35
til
...
4O
HI w
f 1
III 111
m ijj
Hi
-100
-200'
Figure 4.14
1O
15
2O
25
3O
35
4O
141
Table 4.5 Convergence and performance of the solution, (a) updating switching
terms, (b) constant Jacobian
Test no.
CPU time
(seconds)
1a
1*>
2a
2b
3a
3b
4a
45
5a
5b
15.5
11.1
24.9
14.8
17.7
12.3
22.7
14.1
31.5
20.8
Switch
iterations
Main
iterations
-
6
6
6
6
6
9
9
6
6
6
8
11
5
7
7
9
11
21
X
0.1314
0.1712
0.1440
0.2562
0.1069
0.1130
0.2498
0.3129
0.3068
0.5382
Jacobian might be expected. The execution times listed in Table 4.5 are for
a Sun Sparcstation IPX.
4.8 Summary
This Chapter describes the use of iterative harmonic analysis for the
solution of static converters embedded in AC systems. The electric-circuit
constraints, control and switching functions are incorporated in a unified
way, and solved by the Newton method.
The steady-state interaction of a 12-pulse converter with weak AC and DC
systems has been analysed. In the test system, the AC side has been
represented by a three-phase frequency-dependent Thevenin equivalent,
and the DC system by a harmonic Norton equivalent.
The results of the Newton solution have been compared with a timedomain simulation of the same test system. The close agreement between
two such different methods of solution for the test system validates both of
them. Since the unified method converged quickly for such a difficult test
system, convergence is highly likely for any realistic system.
4.9 References
1 CALLAGHAN, CD., and ARRILLAGA, J.: 'Convergence criteria for iterative
harmonic analysis and its application to static converters'. ICHPS-IV, Budapest,
Hungary, October 1990, pp.38-43
142
2 SZECHTMAN, M., WEISS, T., and THIO, C.V.: 'First benchmark model for HVdc
control studies', Electra, April 1991, (135), pp. 55-75
3 SMITH, B.C., WATSON, N.R., WOOD, A.R., and ARRILLAGA, J.: 'A solution for
the steady state interaction of the ac/dc converter with weak ac and dc systems'.
Trans. IEEE, Paper 96-SM 447-3 PWRD
4 LARSON, E.V., BAKER, D.H., and McIVER, J.C.: 'Low order harmonic interaction on ac-dc systems', IEEE Trans. Power Deliv., January 1989, 4, (1), pp. 493-501
5 SMITH, B.C.: 'A harmonic domain model for the interaction of the HVDC
converter with ac and dc systems'. PhD thesis, University of Canterbury, New
Zealand, 1996
6 ZOLLENKOPF, K.: 'Bifactorisationbasic computational algorithm and programming techniques'. Conference on Large sets of sparse linear equations, Oxford,
England, 1970
7 PRESS, W., TEUKOLSY, S., VETTERLING, W., and FLANNERY, B.: 'Numerical
recipes in FORTRAN, the art of scientific computing' (Cambridge University
Press, 1992, 2nd edn.)
Chapter 5
5.1 Introduction
The power-flow solution, described in Chapter 3, is based on the assumption of symmetrical and undistorted operation, the latter achieved by
perfect filtering at the converter terminals. In practice, there is always some
asymmetry in the AC system which, as described in Chapter 2, generates
noncharacteristic harmonics. These harmonics can excite resonances and
jeopardise the operation of control and protection equipment, as well as
involve extra power loss which should be accounted for in the power-flow
solution.
Thus, a three-phase fundamental and harmonic-frequencies solution of
the AC system and converter is required. The elements for such a solution
are the three-phase power flow,1'2 three-phase harmonic representation of
the AC system3 and the three-phase harmonic converter model. 4 " 6 In the
models proposed so far,7'8 the problem is partitioned into several modules
which are then solved separately. For example, a separate Newton-Raphson
module is used at each iteration to find the switching instants. Interaction
between these modules is then solved iteratively.
Convergence of the sequential method depends directly on the degree
of decoupling between modules; even a relatively small degree of coupling
can cause divergence. For a converter operating under conditions of
unbalance, there is no a priori justification as to why a sequential method should
work. This being the case, a Newton solution of all the unknowns is preferred.
In this Chapter, a three-phase power-flow and harmonic-converter model
are combined and solved using both the sequential and unified Newton's
methods.
144
145
At the slack bus a similar argument holds, except that in this case both
Cartesian components of the positive-sequence voltage are specified, rather
than just the magnitude. The presence of negative and zero-sequence
currents in the network is primarily due to sequence transformation in
transmission lines and unbalanced loads. Load buses are, therefore, best
specified by independent real and reactive powers in each phase.
I*
+ sw}
** =
F =
-/{v e ,i5 -
(5.2.1)
where the subscripts a, b, and c refer to the phase, I to the bus and V and I
are the nodal-phase voltage and current r.m.s. phasors. V and I are related
by the nodal admittance matrix for the system
I = YV
(5.2.2)
F 2/ = | V S / - V + / V * ;
(5.2.3)
Four more equations are required. These are obtained with reference
to Figure 5.1, by writing KirchofFs current law for the negative and
zero-sequence currents at the generator terminal and taking Cartesian
146
positive sequence
V
negative sequence
zero sequence
sP
VI
Figure 5.1
components
F 5 , = R{V0lY0l
F 6 i = /{V 0 ,Y 0 ,
(5.2.4)
At the slack bus, the positive-sequence voltage magnitude and phase are
specified, which is equivalent to specifying both the real and imaginary
parts
2(
= /{Vslack-V+i}
(5.2.5)
(5.2.6)
Y = L,
Y,~l
(5-2-7)
In the following, the bar convention is used to indicate a real tensor matrix,
the subscripts R and I to refer to the real and imaginary parts of a complex
number and the subscripts / and m to refer to the number of a bus. A
submatrix of a larger matrix will be referred to by specifying the range of
row and column indices. For example, H = ](a:btC:d) is a matrix consisting of
rows a to b, and columns c to d of the matrix J.
The advantage in using real matrices as opposed to complex numbers in
this analysis, is that the conjugation operator is replaced by a left-multiplicative constant,
a*o[*]a
(5.2.8)
where
(5.2.9)
This makes conjugation differentiable in terms of real Cartesian components.
With nodal voltages in phase components, partial derivatives of mismatch
equations at PQ buses are easiest to derive. Assume that / is the bus number
of a PQ bus. There are six power mismatches at this bus, labelled Fkl, where
k = 1,2,3,4,5,6. The power mismatches at bus / depend on the nodal
voltages, in all three phases at every bus, m, to which it is attached (including
itself). Let Yb = YR +jYY be the branch admittance connecting phase p of
bus / to phase q of bus m. Then assuming that p " q, or that I ^m, the
following partial derivatives hold
3Vi
-R-K M - l
(5.2.10)
148
RYR
(2p-l)l _ y y
vV
(5.2.11)
V.Y.-IR
<5V,
(5.2.12)
oi
1 a
AV + 1 I
AV 0IR
AV,
on
AV_ IR
AV_,
2 0
0 2
2 0
0 2
2 0
0 2
-1
73
-73
-1
-1
-73
73"
"AV O I R
-1
AV a ( I
AVWR
IsO
AV + IR
AV W I
-1
73
-1
-73
AV C / R
-1
AVd,
-73
-1
(5.2.13)
+,
AV,
0(
=T
(5.2.14)
dV,6(1
0
-2
0
2
0
2
-V3 - 1
-1
1
0
2
2
0
2
0
-1
-1
-V3
V3
-1
/Q
1
0
2
/SI
-l
AI a i R
AI a i l
AIWR
AIM,
(5.2.16)
AI C(R
AI,n
= T
(5.2.17)
AIC,
AI-,
T h e chain rule can now be applied to relate changes in the phase voltage at
other buses to the right-hand side of eqn. 5.2.16. If Ylm is the 3 x 3
admittance connecting bus / to bus m, given by
(5.2.18)
3>21
}>33
150
"Jill
Jl2R
-Jill
Jl3R
Jill
JllR
Jl2I
Jl2R
3*131
Jl3R
J21R
"Jill
J22R
-J22I
J23R
-J23I
Jill
J21R
J22I
J22R
J23I
J23R
(5.2.19)
"AVflm"
(5.2.20)
= T'Y AV hm
AIOi
AV cm
(5.2.21)
If / 7^ m, t h e n t h e last four rows of T ' Y a r e t h e partial derivatives
= T'Y( 3:6>1:6 )
(5.2.22)
= (
0
0
0
0
0
0
0
0
0
0
0
0
Y0R
-Y0I
Y0R
0
Yo,
0
Y-R
Y_,
Y_
(5.2.23)
(3:6,1:6)
(5.2.24)
m=l
(5.2.25)
"T
* / \(\ 1 fi'i
ri^/^
, riVx
' ' ^*
1
...
(5.4.26)
(5.2.27)
F 2( = V s 2 p -[fV,T'V ( ] (1>1)
(5.2.28)
*m
(5 2 29)
F 2( = /{V s l a c k -V + i }
(5.2.30)
(5.2.31)
The partial derivatives with respect to terminal voltages are, therefore, just
^ < - = {f:V'OTO:O'O] _ v * f
\ * l (5.2.32)
152
output results
Figure 5.2
Flow diagram for the three-phase power flow; FNIT is the number of
full Newton iterations, MAXIT is the maximum number of iterations
(5-2.33)
In the above derivation of 6 x 6 Jacobian blocks, many redundant calculations have been specified which are not carried out in practice. For example,
in eqn. 5.2.26 only the first row of [V+lT'Ylm] is calculated.
(5.2.34)
154
ROXBURGH-011
(Slack)
ROXBURGH--220
P+JQ
INVERCARG011
P+JQ
INVERCARG220
TIWAI220
Figure 5.3
The Jacobian for the test system is illustrated in Figure 5.4 by a meshplot
of the magnitude and a sparsity plot of the location of elements in the
matrix. Structural sparsity would be much greater for a realistically-sized
system and is clearly not symmetric. The largest elements are associated
with PQ mismatches at buses with the most connected lines, or the highest
AC-system admittance.
In order to assess the performance of the power flow, a variety of test cases
are solved. The base case is that described in Table 5.1, which contains
severely unbalanced loads. After converging in four iterations, the sum of
the magnitudes of all the mismatches, the 1-norm residual, was
2.2 x 10~ 10 p.u. This was for a solution in which the Jacobian was updated at
every iteration and with a starting point of 1 p.u. positive-sequence voltage at
every PQ bus. Updating the Jacobian only twice increased the number of
iterations to seven, with a final 1-norm residual of 3.9 x 10 ~7. The nodal
voltages at the solution are given in sequence components in Table 5.2. Note
that the voltage constraints at generator buses have been satisfied.
Bus name
Bus type
INVERGARG220
PQ
INVERGAR011
PQ
ROXBURGH-220
PQ
MANAPOURI220
PQ
MANAPOURI014
PV
TIWAI-220
PQ
NTH-MAKAR220
PQ
INV-MAN--TMP
PQ
ROXBURGH-011
Slack
Constraints
Sa =
Sb =
Sc =
Sa =
Sb =
Sc =
Sa =
Sb =
Sc =
Sa =
0+y0MVA
0+y0MVA
0+y0MVA
50+y15MVA
45 + y'20 MVA
52+y18MVA
48+y20MVA
48+y12MVA
51.3 + y'28.3 MVA
0+y'0MVA
s c = o + yo MVA
P+ =500MW
|V+| = 1.05 p.u.
S a =170+y85MVA
Sb=160+y85MVA
Sc=170+y70MVA
Sa = 0+y0MVA
s b = o + yo MVA
s c = o + yo MVA
Sa = 0-fy0MVA
sb = o + yo MVA
Sc = 0+y0MVA
V + = 1.05^0p.u.
Placing a star-g/delta at the slack bus effectively shifts the angle reference
by 30, so that the starting voltages are offset by 30. In this case, convergence required 23 iterations to a final residual of 2.4 x 10" 6 if the Jacobian
was updated twice, five iterations for a full Newton solution and seven
iterations if updated four times. The large number of iterations required in
the first case was due to the large initial mismatch, implying that the
Jacobian calculated after the second iteration was significantly different to
that at the solution. A more versatile approach is to update the Jacobian at
every iteration until the largest mismatch is smaller than a preset tolerance,
at which point the Jacobian is held constant.
If all transformers are star-g/star-g, except for that at the INVERCARGO11 load bus, which is made star-g/delta, convergence is much more
difficult. To obtain convergence, it is necessary to start the voltage at
I <8
'
;.
!
:
:f
:
: :
:.
'.
99.
jf.tM9.M9.
t9.
..
Equation number
:_
Magnitude
INVERCARG220
+ve
-ve
zero
INVERCARG011
+ve
ve
zero
ROXBURGH-220
+ve
ve
zero
MANAPOURI220
+ve
-ve
zero
MANAPOURI014
+ve
ve
zero
TIWAI-220
+ve
-ve
zero
NTH-MAKAR220
+ve
-ve
zero
INV-MAN--TMP
+ve
-ve
zero
ROXBURGH-011
+ve
-ve
zero
Phase ()
PQ
0.96285
0.01264
0.00043
-7.06927
-166.75234
38.40967
PQ
-7.98316
-165.72855
29.61255
0.95719
0.01287
0.00491
PQ
-1.78554
131.37500
-64.09227
0.99201
0.00284
0.00240
PQ
-1.63163
-156.61677
17.09903
1.02863
0.00362
0.00077
PV
1.02177
-156.61677
17.09903
1.05000
0.00181
0.00038
PQ
-8.10666
-166.20626
72.62988
0.95139
0.01441
0.00512
PQ
-6.65871
-166.16024
47.6068
0.96865
0.01212
0.00409
PQ
-6.65427
-166.78061
43.49631
0.96813
0.01197
0.00413
slack
1.00000
0.00142
0.00120
0.00000
131.37500
-64.09227
158
v.i* = s0
VCIC* = Sc
(5.2.35)
(v 0 +v + >a* + n ) = sa
(V0+a2V+)(a2I*
+al*) = Sb
(Vo + aV+)(aI* + a 2 l* ) = Sc
(5.2.36)
159
(5.2.37)
where
s 0 = s a + sb + sc
$+=Sa+a$b+a%
S_ = Sfl + a2Sb + aSc
(5.2.38)
This result confirms that there should be two solutions for the zero-sequence
voltage, if there is no path for zero-sequence current. To investigate this
effect further the test system of Figure 5.6 was solved using the three-phase
power flow.
In this system, the leakage reactance of the star-g/star transformer was set
to 0.01 p.u., and the leakage reactance of the star-g/star-g transformer was
varied from 10 p.u. to 0.01 p.u. in 100 steps. The load was held constant at
Sa = 49.5 -h/19.5 MVA, Sb = 50.5 + j"20.0 MVA and Sc = 50.5 + J20.5 MVA.
With the zero-sequence leakage reactance of the star-g/star-g transformer set to 10 p.u., two solutions to the power flow were obtained,
corresponding to approximately antiphase zero-sequence voltages at the
load bus. Starting from the larger of the zero-sequence voltage solutions, the
leakage reactance of the star-g/star g transformer was reduced in 100
steps down to 0.01 p.u., a more normal situation. The solution of each
intermediate power flow was used as the starting point for the next. The
outcome of this process is plotted as the top curve in Figure 5.7, which is a
plot of zero-sequence voltage magnitude as a function of leakage reactance.
The other two curves in Figure 5.7 show the effect of sliding the leakage
reactance down from 10 p.u. at the lower of the two zero-sequence voltagemagnitude solutions, and up from 0.01 p.u., starting with the zero-sequence
voltage set to zero. In the middle is a region of divergence. It is clear from
this diagram that there are always two solutions at load buses but that in
160
P+jQ
Figure 5.6 A test system with a variable path for zero-sequence current
normal circumstances one of the solutions is too far away from a purely
positive-sequence starting point to be converged towards. As the zerosequence impedance looking into the AC system from a load bus is
increased, the two solutions move closer together until they are of approximately equal magnitude. If the zero-sequence impedance is greater than
about 1 p.u., either solution can be converged towards, depending upon the
starting point.
In any real AC system there will always be a path for zero-sequence
Figure 5.7
4
5
6
zero seq impedance (p.u.)
10
ac system
measured dc current
alpha order
162
The basic control types represented in the present model are constant
current, constant power and minimum gamma. In the case of constantcurrent control, the controller is assumed to be a simple PI type, with a
low-pass current-transducer characteristic. The twelve firing and end-ofcommutation instants are included as part of the Newton solution. Under
current control, the firing instants are, therefore, modulated by ripple in the
DC current. The end-of-commutation instants are modulated owing to
distortion in the terminal voltage and DC current.
The firing controller is assumed to be referenced to an ideal PLO
tracking the terminal voltage. This means that all firings will be perfectly
equispaced if there is no ripple in the DC current, despite any distortion in
the terminal voltage.
The AC-side terminal voltage is described in harmonic-sequence components. Sequence components are chosen as there is usually no interaction
between zero-sequence voltage harmonics and the converter, thereby improving sparsity in the Jacobian matrix.
(5.4.1)
converter
Newton
solution
yes
calculate and factorise
converter Jacobian
use Jacobian to update
converter variables
calculate converter
complex power
evaluate power flow mismatches
power
flow
Newton
solution
calculate and factorise
power flow Jacobian
use Jacobian to update
load flow variables
calculate power-flow
Thevenin-equivalent
internal voltage
Figure 5.9
163
164
Converter Representation
Sc-P+jQ
Linearising
Inductor
(a)
Load-flow Representation
Vth
Converter Equations
(b)
Converter Equations
Figure 5.10
If the converter model is now solved with this new Thevenin voltage, an
updated value of Vc and I c will be obtained. For the next solution of the
power flow (Figure 5.10a), the specified power at the converter bus must be
updated to be
SCy = VCvI*
(5.4.2)
Likewise, when the Thevenin voltage for the converter model is calculated,
*>-%-%,
(5 4 3)
- -
166
Fundamental
Phase Voltages
Figure 5.11
frequency is equal to the current flowing out of the converter and any shunts
connected there
Ft/c = Ii + I 2 + I 3
(5-4.4)
where I x , I2,13 are defined in Figure 5.10c. The power-flow and converter
modules can readily supply these currents. The Jacobian for the unified
system is illustrated in Figure 5.11 as a plot of the structure. All other
elements in the Jacobian are small and have been discarded to improve
sparsity. The terms associated with current mismatch at fundamental frequency appear as horizontal and vertical bands coincident with the converter bus TIWAI
220. These terms, being in-phase components, contribute
approximately 300 elements to the Jacobian matrix over and above those
associated with the decoupled power flow and converter Jacobians. The
overheads associated with combining the two models into a single solution
are therefore minimal, especially since only one set of sparse storage vectors
need be defined and manipulated instead of two.
test
test
test
and
The voltage base for this test system was 220 kV, the power base was
100 MVA, the system impedance at the converter bus was 0.03 Z_70, the
converter power was approximately (610 +j'75)MVA and it was running
under constant-current control as a rectifier. The rectifier DC system and
filters corresponded to those at the rectifier end of the CIGRE benchmark
system. The rectifier delay angle was 15.
168
Table 53
Test no.
Type
NJR
TNIT
Time
J size
unified
dcp A
dcp B
dcpC
unified
dcp A
dcp B
dcpC
unified
dcp A
dcp B
dcpC
1
7
8
15
1
7
9
15
1
8
9
*
4
7
8
15
4
7
9
15
5
8
9
*
3.7
24.1
19.8
36.1
3.6
22.0
19.9
35.0
3.9
30.1
20.4
6971
6540
6542
6540
6991
6542
6542
6544
7355
7008
7008
*
type B where only one iteration of Newton's method was applied to the
converter and power-flow systems at each decoupled iteration. With the
linearising shunt absent, the decoupled method failed to converge in test 3,
when no filters were present. For all cases recorded in Table 5.3, the starting
point was from a converged, fundamental-frequency only, solution.
Convergence of the decoupled method is entirely dependent on the size
of the partial derivatives of the quantities that are being iteratively solved
for: in this case, the Thevenin equivalent voltage, Vth, and the converter
equivalent power, Sc. In three phases, and real and imaginary parts, there
are twelve quantities involved. The decoupled method is a fixed-point
iteration for these quantities
[SC]N+1 = converter (\Vthf)
[\y
N+1
= powerflow ([SJ
N+1
(5.4.5)
)
(5.4.6)
(5.4.7)
The average size of these partial derivatives is 0.3779 and 0.4178, for the
10
Figure 5.12
Iteration
20
Converter Convergence
10
Iteration
20
decoupled methods with and without the linearising inductance, respectively. Convergence of the decoupled methods is, therefore, uncertain at
best.
Even when the decoupled method does converge, a consequence of the
large coupling derivatives is that at each iteration the gains made by each
Newton subsystem are reduced. This leads to an oscillatory convergence
characteristic, illustrated in Figure 5.12. Figure 5.12 is a plot of the sum of
the magnitude of the mismatches (the 1-norm) of the power-flow and
converter systems at every step of the decoupled and unified methods.
Again, the advantage in using a linearising shunt is evident. A consequence
of the oscillatory convergence is that Jacobian matrices for the converter
and power flow must be recalculated at each iteration of the decoupled
method. This is the most significant contributing factor in the long execution time of the decoupled method, since the Jacobian must also be
refactored.
5 A3 Powerflow/converterinteraction
If interaction between the three-phase power flow and converter model is
not solved simultaneously, some degree of error must exist in the solution
obtained. For example, a common approach is to estimate the converter
170
10
15
20
25
30
35
Be
Left - Exact
Right - Approximate
15
20
25
Harmonic Order
Figure 5.13
Exact
Approx.
PQ
1.01217
0.00876
0.00780
1.01211
0.00705
0.00183
PQ
1.04252
0.00150
0.00044
1.04250
0.00077
0.00040
PQ
1.04470
0.00149
0.00036
1.04470
0.00192
0.00048
PV
1.05000
0.00096
0.00024
1.05000
0.00075
0.00018
Converter
1.00618
0.00803
0.00182
1.00625
0.01000
0.00989
PQ
1.01541
0.00673
0.00165
1.01547
0.00840
0.00741
PQ
1.01499
0.00655
0.00150
1.01504
0.00817
0.00720
Slack
1.05000
0.00075
0.00022
1.05000
0.00038
0.00020
5.5 Summary
To make it compatible with the converter-harmonic model of Chapter 2 and
thus enable an integrated Newton solution of the power flow and harmonic
interaction, the three-phase power flow described in this Chapter uses
172
Cartesian, rather than polar co-ordinates for the mismatches. The overall
solution is formulated in terms of phase components, which suits the P, Q
bus specifications; however, at generator buses, the mismatches are in
sequence components.
Interaction between the three-phase power flow and a three-phase harmonic converter model has been solved using both decoupled and full
Newton methods. The decoupled method displays good convergence if a
linearising shunt is present but is slow because the Jacobian matrices need
to be recalculated and factorised at every iteration. The full Newton method
(with constant Jacobian), is faster and displays more robust convergence.
The decoupled method is compatible with any existing three-phase power
flow, including a polar fast-decoupled type. A special but simple three-phase
power flow is required for the unified method as the converter model must
be framed in real variables. Overall, code for the unified method is shorter
and less complicated, as there is only one set of sparse storage, mismatch
evaluation, convergence checking, etc.
The effect of not modelling the power-flow/distortion interaction is a
moderate underestimate of unbalance in the power-flow solution, and a
large underestimate of distortion at low order noncharacteristic harmonics
in the converter model.
The AC-DC iterative algorithm can easily be extended into a generalpurpose model with the capability of including several AC systems, DC
systems, possibly with multiple terminals in each system, and all integrated
with the power-flow equations. This is essentially a software engineering
task, as each half pole contributes a block to the main diagonal of the system
Jacobian matrix, with diagonal matrices coupling to other harmonic sources
in the same system. The main task is to code the program so as to be
versatile, easy to use and modular.
5.6 References
1 ARRILLAGA, J., and HARKER, B.J.: 'Fast decoupled three-phase load flow',
Proc. IEE, August 1978, 125, (8), pp. 734-740
2 WASLEY, R.G., and SLASH, M.A.: 'Newton-Raphson algorithm for three phase
load flow', Proc. IEE, July 1974, 121, (7), p. 630
3 DENSEM, T.J., BODGER, P.S., and ARRILLAGA, J.: 'Three phase transmission
systems for harmonic penetration studies', IEEE Trans. Power Appar. Syst, February 1984, 103, (2), pp. 310-317
4 YACAMINI, R., and de OLIVEIRA, J.C.: 'Harmonics in multiple convertor
systems: a generalised approach', IEE Proc. B, March 1980, 127, (2), pp. 96-106
5 SMITH, B.C., WATSON, N.R., WOOD, A.R., and ARRILLAGA, J.: 'A Newton
solution for the steady state interaction of ac/dc systems', IEE Proc, Gener. Transm.
Distrib., March 1996, 143, (2), pp. 200-210
6 SMITH, B.C., WATSON, N.R., WOOD, A.R., and ARRILLAGA, J.: 'A solution
for the steady-state interaction of the ac/dc converter with weak ac and dc
systems'. ICHPS Conference, Las Vegas, October 1996
Chapter 6
6.1 Introduction
The simulation of electromagnetic transients associated with HVDC
schemes is carried out to assess the nature and likely impact of overvoltages
and currents, and their propagation throughout both the AC and DC
systems. Transient simulation is also performed for the purpose of control
design and evaluation. Transients can arise from control action, fault
conditions and lightning surges. The response to these events is initially
dominated by electrical resonances in lumped RLC components, and
travelling waves in distributed-parameter components such as overhead
lines and cables. Electromechanical transients, usually studied over longer
timescales, are considered in Chapters 7 and 8.
The complexity of the DC system and the AC-DC system interactions has
necessitated the development of special tools capable of adequate representations. The principal tools are DC simulators (transient network analyser,
TNA), and digital computers.
The DC simulator has long been the most powerful tool for the design of
HVDC systems. It is a versatile tool for simulating dynamic performance
under various conditions, ranging from fundamental frequency to high
frequency; however, it is costly and less flexible than digital simulation. With
the recent advances in digital computer technology, computers are playing
an increasingly important role. At present, both DC simulators and digital
computers are used, often in a complementary manner, in the design and
analysis of HVDC systems.1
The DC simulator is characterised by the scaled physical modelling of all
the power circuits associated with an HVDC system, although actual control
circuits can be used. Every component in the system is modelled either by
its scaled-down physical model or by an electronic model. Being a real-time
model, it is capable of simulating a large number of cases in a short time.
The DC simulator does not need any approximations or simplifications in
the bridge model, and can cover a wide range of frequencies of interest in
real-time operation. However, the extent of AC-network representation is
limited, and the frequency response of the simulator is limited to about that
of the switching surges and below. It is possible to minimise the AC-network
176
177
178
(L,ineqn. 6 2.15)
Figure 6.1
179
180
x=f(t,x)
which for a linear time invariant system becomes
x = [A]x + [B]u
(6.2.1)
where x is the vector of state variables, [A] is the system matrix and u the
vector of sources. The system matrix can provide useful information about
the range of time constants and resonant frequencies in the system. A
summary of state-space analysis is provided in Appendix IV. Eqn. 6.2.1 is
solved by numerical integration or evaluation of an explicit closed-form
solution. Numerical-integration techniques are discussed in Appendix V.
(iii) The state variable derivative xn + 1 is then estimated from state equation
*n + l=/0n + M* w + l)
(iv)
(6.2.3)
(6.2.4)
181
j:i-*J+i)}
(6.2.5)
182
source approach can cause very large voltage spikes when trying to overcome inaccurate initial conditions. Partitioning the inductor fluxes into state
and dependent variables is complicated and time consuming. An inductor
can still be dependent even if it is not connected directly to a radial node
consisting of inductive branches but has an intervening resistor/capacitor
network.
The identification of state variables can be achieved by developing a
node-branch incidence matrix where the branches are ordered in a particular pattern (e.g. current sources, inductors, voltage sources, capacitors,
resistors) and Gaussian elimination performed. The staircase columns
represent state variables.12 However, this is computationally expensive if it
is required frequently, such as in the case of converter switchings.
the vectors Va and Ilf representing the capacitive node voltage and inductive
branch currents, respectively, become dependent variables related to the
state variables by algebraic relations.
The resulting set of equations can be written as
state-related variables
(6.2.7)
(6.2.8)
183
h0
dependent variables
Tf __
TO
rjz
(6.2.9)
K{.V
(6.2.10)
(6.2.11)
'Xr
(6.2.15)
(6.2.16)
E, V and I represent the branch e.m.f. source voltage, node voltage and
branch current vectors, respectively.
R, L and C are the resistance, inductance and capacitance matrices,
respectively. The p operator in TCS represents a derivative w.r.t. electrical
angle rather than time.
The rest of the symbols and suffixes in the above equations represent the
following
a = node with at least one capacitive branch connected
/? = node with at least one resistive but no capacitive branches
connected
184
(6.2.17)
where
(6.2.18)
185
(a)
Time
The dependent state variables are then calculated at tx from the state
variables, and written to the output file. The next integration step will then
begin at tx with step length h as shown in Figure 6.3fe. This linear
approximation is sufficiently accurate over periods which are generally less
than one degree, and is computationally inexpensive. The effect of this
interpolation process is clearly demonstrated in a case with an extended
1 ms time step in Figure 6.5.
Upon switching any of the valves, a change in the topology has to be
reflected back into the main system network. This is achieved by modifying
the connection matrices. When the time to next firing is less than the
integration step length, the integration time step is reduced to the next
closest firing instant. Since it is not possible to integrate through discontinuities, the integration time must coincide with their occurrence. These
discontinuities must be detected accurately since they cause abrupt changes
in bridge-node voltages, and any errors in the instant of the topological
changes will cause inexact solutions.
Immediately following the switching, after the system matrices have been
reformed for the new topology, all variables are again written to the output
file for time tx. The output file therefore contains two sets of values for tx
immediately preceding and after the switching instant. The double solution
at the switching time assists in forming accurate waveshapes. This is
especially the case for the DC-side voltage which almost contains vertical
jump discontinuities at switching instants.
186
N\
I
(a) SO /is time step
N\
(b) 1 ms time step
187
Time
Figure 6.5
The firing angle was fixed at 25. Figure 6.4 shows the valve voltages and
currents for 50 fi$ and 1 ms (i.e. 1.0 and 21 degrees) time steps, respectively.
The system has achieved steady state even with steps as large as 20 times.
The progressive time steps are illustrated by the dots on the curves in
Figure 6Ab9 where interpolation to the instant of a valve-current reversal is
made and from which a half time-step integration is carried out. The next
step reverts back to the standard trapezoidal integration until another
discontinuity is encountered.
A similar case with an ideal AC system terminated with a DC source was
simulated using TCS. A maximum time step of 1 ms was used also in this
case. Steady-state waveforms of valve voltage and current derived with a
1 ms time step, shown in Figure 6.5, illustrate the high accuracy of TCS
both in detection of switching discontinuities and reproduction of the
50 /*s results. The time-step tracing points are indicated by dots on the
waveforms.
Further TCS waveforms are shown in Figure 6.6 giving the DC voltage,
valve voltage and valve current at 50 jus and 1 ms (although the actual step
was modulated under the program's control).
In the NETOMAC case, extra interpolation steps are included for the 12
switchings per cycle in the six-pulse bridge. For the 60 Hz system simulated
with a 1 ms time step, a total of 24 steps per cycle can be seen from the
Figure 6.46, where a minimum of 16 steps are required. The TCS cases
shown in Figure 6.6 have been simulated with a 50 Hz system. The 50 jus
case of Figure 6.6a has an average of 573 steps per cycle with the minimum
requirement of 400 steps. On the other hand, the 1 ms time step needed
only an average of 25 steps per cycle. The necessary sharp changes in
188
- valve I cunem(pu)
Rectified dc voliage(pu)
valve 1 voltage(pu)
rJVW
1
l\
ri
'
1
\
0.460
n/"i/
0.468
Figure 6.6
\
\
0.476
0.484
VV"
0.492
0.500
0.468
0.476
0.484
Time(s)
Time(s)
0.492
0.500
waveshape are derived directly from the valve voltages upon topological
changes.
When the TCS frequency was increased to 60 Hz, the 50 jus case used
fewer steps per cycle, as would be expected, resulting in 418 steps compared
to a minimum required of 333 steps per cycle. For the 1 ms case, an average
of 24 steps were required, as for the NETOMAC case.
The same system was run with a constant-current control of 1.225 p.u.,
and at 0.5 s a DC short circuit was applied. The simulation results with 50 fis
and 1 ms step lengths are shown in Figure 6.7. This indicates the ability of
the TCS to track the solution and treat waveforms accurately during
transient operations (even with such an unusually large time step).
0.48
dc fault application
x Rectified dc voitage(pu)
A dc current(pu)
0.60
0.52
0.54
0.56
Time(s)
(a) 1 ms time step
Figure 6.7
0.58
0.60
0.48
dc fault application
x Rectified dc voltage(pu)
A dc current(pu)
0.50
0.52
0.54
0.56
Time(s)
(b) 50 jus time step
0.58
0.60
189
NPLO
"order
Converter
Firing
Control System
VALFIR
I
Valve States
Feedback
Control
System
orP?
off
' measured
c .(fc'
'POggO'^-
190
Zero crossover points are detected by the change of sign of the reference
voltages and multiple crossings are avoided by allowing a space before
making the next crossing. Distortion in line voltages can cause difficulties
in proper zero-crossing detection, and the voltages must be smoothed
before being passed to the zero-crossing detector.
The time between two consecutive zero crossings, of the positive to
negative (or negative to positive) going waveforms of the same phase, is
denned here as the half period time, T/2. The measured periods are
smoothed through a first-order real-pole lag function with a user-specified
time constant. From these half-period times the AC-system frequency is
estimated, every 60 (30) for a six (12) pulse bridge.
Normally, the ramp for firing of the particular valve (c(l).. .c(6)) starts
from the zero crossing points of the voltage waveform across the valve. After
T/6 time (T/12 for 12 pulse), the next ramp starts for the firing of the next
valve in sequence. So, a of 0 for the incoming valve corresponds to the end
of T/6 time from the previous ramp, and so on.
It is possible that during a fault or due to the presence of harmonics in
the voltage waveform, the firing does not start from the zero crossover
point, resulting in a synchronisation error, B2, as shown in Figure 6.9. This
error is used to update the phase-locked oscillator which, in turn, reduces
the synchronising error, approaching zero at steady-state conditions. The
synchronisation error is recalculated every 60 (every 30 for 12 pulse).
The firing-angle order (oeorder) is converted to a level to detect the firing
instant as a function of the measured AC frequency by
Figure 6.9
q order (rad)
/AC(P.U.)
(6.2.19)
191
As soon as the ramp c(n) reaches the set level specified by T o , as shown in
Figure 6.9, valve n is fired and the firing pulse is maintained for 120. Upon
having sufficient forward voltage with the firing pulse enabled, the valve is
switched on and the firing angle recorded as the time interval from the last
voltage zero crossing detected for this valve.
At the beginning of each time step, the valves are checked for possible
extinctions. Upon detecting a current reversal, a valve is extinguished and
its extinction angle counter is reset. Subsequently, from the corresponding
zero-crossing instant, its extinction angle is measured, e.g. at valve 1 zero
crossing, y2 is measured, and so on. (Usually, the lowest gamma measured
for the converter is fed back to the extinction angle controller.) If the
voltage-zero crossover points do not fall on the time step boundaries, a
linear interpolation is used to derive them.
As illustrated in Figure 6.8, the NPLO block co-ordinates the valve-firing
mechanism, and VALFIR receives the firing pulses from NPLO and checks
the conditions for firing the valves. If the conditions are met, VALFIR
switches on the next incoming valve and measures the firing angle, otherwise it calculates the earliest time for next firing to adjust the step length.
Valve currents are checked for extinction in EXTNCT and interpolation of
all state variables is carried out. The valve's turn on time (Ton) is used to
calculate the firing angle and off time (Toff) is used for the extinction
angle.
To verify the operation of the firing controller implemented in the TCS
program, two cases are presented in Figures 6.10-6.12 using the same test
system as in section 6.2.6. A constant aorder of 15 was used in the case of
Figures 6.10 and 6.11, where a step change was applied at 0.5 s for five
cycles by increasing the aorder to 75.
- d c Current
Rectified dc voltage(pu)
Firing angle(rad)
Extinction angle(nid)
Rectified dc voltage
A dc Current
3.0
. . .
.
1.5
1
i f. u j i u .u n . ' i 1 V
v-
V-J
iiiinii
0.0
i.i.iih
iilillii
-1.5
0.450 0.475 0.500 0.525 0.550 0.575 0.600 0.625 0.65
Time(s)
(a) TCS
niiiii
Time(s)
response
192
0.45
0.50
0.55
0.60
0.65
Rectified dc voltage
0.70
0.50
0.55
Time(s)
0.60
0.65
0.70
0.75
Time(s)
(a) TCS
response
o dc Current(pu)
~~ dc Current(pu)
2.0 T
0.45
0.50
0.55
0.60
'"
'
0.65
0.70
Time(s)
(a) TCS response
-0.5
0.75
0.45
0.50
0.55
0.60
0.65
0.70
0.75
Time(s)
(b) EMTDC response
193
vm
Figure 6.13 Resistive branch
1"
R
(6.3.2)
dikm(t)
dt
(6.3.3)
- AO = i f
- \m(x)]dx
(6.3.4)
194
*- At
as illustrated in Figure 6.14. The right-hand side of eqn. 6.3.4 is approximated by the shaded trapezoidal, yielding
2^,(0-v m (0]
2L'
(6.3.5)
The branch current, therefore, consists of two components, one related to
the nodal voltages at the current time step which are to be solved, and the
other independent of the current-time step voltages. The independent
component can be represented as a current source, and the dependent
component as a branch resistance, yielding the Norton equivalent, shown
Utm
2L
A/
195
(6.3.6)
where
-2L
2L
A*
(6.3.8)
dt
(6.3.9)
/ - At)]
(6.3.10)
2C
2C
-g \yma -At)-
yk(t -
- At)
(6.3.11)
where
Ifcm(* - At) = ^
(6.3.12)
196
AL
2C
Figure 6.16 Capacitive branch
-Af
2C
Vk
2C
di
"dx~ Jt
dx
dt
(6.3.14)
(6.3.15)
(6.3.16)
where z = ^/L/C and s is the wave speed. Consider a wave starting at node
m and moving towards node k, as shown in Figure 6.17.
Multiplying eqn. 6.3.16 by z and adding it to eqn. 6.3.17 eliminates the
negative travelling wave, i.e.
v(t) + zi(t) = 2zf(x - st)
(6.3.18)
The left-hand side of this equation is a voltage that depends only on x st.
197
km
x=0
x=d
k
(6.3.19)
(6.3.20)
Thus
v(t - t) 4- zi(t -
T)
= t/(0 + zt'(0
(6.3.21)
(6.3.22)
(6.3.23)
where
(6.3.24)
Similarly
(6.3.25)
198
Figure 6.18
Transmission-line equivalent
where
(6.3.26)
The transmission line has therefore been replaced by two Norton equivalents (Figure 6.18).
In nodal form, eqns. 6.3.23 and 6.3.24 become
(6.3.28)
The linear system is constructed by combining the nodal equations for every
system component with due regard to the system topology.
Current sources are modelled merely as an additive term to the left-hand
side of eqn. 6.3.28. With voltage sources included, part of the right-hand
side of eqn. 6.3.28 is known, as opposed to the left-hand side. A partitioning
of Y is therefore required with voltage source nodes ordered last
(6.3.29)
where the history current terms have been included in the left-hand side,
VT are the specified voltage sources, and IN are the known Norton currentsource injections. This system is solved by first finding V from
[IN - BVT] = [A]V
(6.3.30)
(6.3.31)
199
Provided that neither the system topology nor the time step change, the
system conductance matrix [Y] is constant throughout the simulation. [Y] is
also sparse, since each node is connected to only a few others. The presence
of transmission lines may also cause [Y] to be block diagonal, affording a
breakdown of the system to be solved into smaller subsystems (section 6.4).
The linear system, represented by eqn. 6.3.28, is solved by Gaussian
elimination and LU decomposition.
A simpler alternative implementation of voltage sources is to assume that
they have a source impedance |zsource| > R ^ , and to model them as Norton
equivalents. This completely avoids the requirement of having to partition
the conductance matrix as in eqn. 6.3.29.
(6.3.32)
rnr^
>
V
- r '^A
L 1
(6.3.33)
onrn
J_ <t
s
m
Figure 6.19 Transmission-line stub solution for an inductor yielded by trapezoidal
integration
200
so that z = 2h/At, where L' and C are the line inductance and capacitance
per unit length. A wave leaving k at t At travels to 5 by t At/2.
From eqn. 6.3.22, in slightly different form
Vkm(t - At) + zikm(t -At) = zis(t - ^j
(6.3.34)
(6.3.35)
(6.3.36)
(6.3.37)
*
tan(a)Z(L'C'))
coV(L'C')_
1
coAt/2
ijteiL'
tan(coA*/2) jcoL
Consequently, the ratio of the impedance of the stub line to the impedance
of the inductor is
(6.3.39)
which has been plotted in Figure 6.20 as a function of the Nyquist frequency
1/2A2 of the simulation.
An error analysis for the trapezoidal integration of a lumped capacitance
is the dual of that for an inductor. The trapezoidal integration yields the
exact solution for an open-circuit stub line of transit time x = A//2 with
C7 = C, z = AJ/2C so that L = A*2/4C. Figure 6.20 is, consequently, also a
plot of the ratio of the actual admittance of the capacitor to the admittance
of the open-circuited stub line calculated by the trapezoidal integration.
201
>
_^_
._.
,j
"1
0.9
0.8
0.6
or
0.4
: Y;
0.3-
0.2-
0.1
Figure 6.20
633
0.1
0.2
0.3
0.4
0.5
0.6
0.7
: i\
0.8
0.9
Switching discontinuities
202
(a)
(b)
Figure 6.21
Switch representation
Unblocked
Blocked
Number of refactorisations
Simulation time
10 fis
50 fis
10 ^s
50 us
2570
2480
1
1
4 mm 41 s
1 min21 s
4 min 24 s
1 min 9 s
203
At (A) the diode should be switched off, but the negative current is only
detected at (B) and the conductance matrix is reformed at (C). This
approach is inadequate for HVDC applications, since any error in the
thyristor switching instants yields noncharacteristic harmonics on the AC
and DC sides. This is due to the fact that the converter acts as a modulator
of the large DC-side current so that a small variation in switching instant
yields a substantial current on the AC side. This effect is important for SSR
studies, configurations involving GTOs and the calculation of uncharacteristic harmonics.
An effective solution of this problem is to apply a linear interpolation at
(B) to find all the nodal voltages at a point very close to (A). An exact
solution for the zero crossing is a nonlinear problem; however, given the
close proximity of points t and t + At, a linear interpolation between them
introduces no significant error.
In the above example, the diode model must include logic which detects
that a switching has occurred between t and t + At, and an estimate, t + xv
of when. For a diode, a linear interpolation on the forward current yields
T =
Mf(t)
if(t)-if(t + At)
(6.3.40)
^ v(t) + [v(t
At
At) - v(t)]
(6.3.41)
204
Figure 6.23 A double interpolation scheme to find the switching instant and
resynchronise time steps
T-2At
T-At
T+At
T+2At
Figure 6.24 A circuit which will display voltage chatter after the diode switches off
205
(a)
T-2A/
T-Af
T+Af
T+2Ar
(b)
v
Figure 6.25 a Correct solution for voltage across the inductor
b Chatter voltage at the diode anode
(6.3.42)
(6.3.43)
post-processing averaging;
half time step interpolation;
critical damping adjustment method;
damping resistors.
206
i>
R
Figure 6.26 Damped inductor method for reducing chatter
-Xv(t-At)
(6.3.45)
where
R - 2L/A/
~ R + 2L/A*
is the decay constant. If R = 2L/A/, X = 0 and there is no chatter. However,
the simulation would be inaccurate because of the relatively small values of
damping resistance. Larger values of damping resistance lead to more
accurate simulations but increasingly persistent chatter after discontinuities.
Setting R = 2L/A/ is called critical damping, and substituting into eqn.
6.3.44
i{t) = i(t -At) + ^v(t)
(6.3.46)
Figure 6.27
207
trapezoidal integration
interpolation back to discontinuity, between (1) and (2)
trapezoidal integration with new conductance matrix
interpolation between (3) and (4) of half time step to
remove chatter
(5) -> (6) trapezoidal integration to get past t + At
(6) (7) interpolation between (6) and (5) to resynchronise time
steps
(7) - (8) normal integration resumes
208
1.0 Q
VL
0.1 H
Diode
0.00
0.01
0.02
0.03
Time(seconds)
0.04
209
10 -|
5 0
-5
10 -i
''.A
-5
0.01
0.02
0.03
Time(seconds)
0.04
6.4 Subsystems
Since the system-conductance matrix is refactorised every time switching
occurs, fast simulation can only be achieved if the factorisation is as efficient
as possible. This entails the use of techniques which utilise the sparsity of
the conductance matrix. The conductance matrix for a system of w-nodes is
of size n x n and, if factored in full, requires of the order of w3 floating-point
operations (0(n3) flops). Only a small fraction of the n2 elements of the
conductance matrix are nonzero, however.
Nonzero elements are present on the diagonal, and wherever two electrical nodes are connected by a nonzero conductance. Typically, each
electrical node is connected to only two others, and therefore contributes
three elements to the corresponding conductance matrix row. As a rough
210
10'
10'
10*
tioa
6
10"
102
10'
/cubic fit
50
100
150
200
250
No. of nodes
300
350
400
450
500
Figure 6.31 Even with sparsity, flops count for conductance-matrix factorisation
increases with the cube of the number of electrical nodes
211
compensation method;
partial-matrix reduction;
Sherman-Morrison formula;
subsystems.
(i) and (ii) are very similar, and have both been implemented in the EMTP
program, but were eventually discarded in favour of full-matrix reduction.
In (i), the system matrix is not modified after each switching. The correct
solution is obtained by injecting compensating currents at the switch
terminals obtained from a reduced Thevenin equivalent of the network.
The partial matrix reduction method confines the matrix reduction to
that part containing switches, by ordering nodes with switches last and
partitioning
ara
(6.4.1)
so that
= [D-CA ^
(6.4.2)
[I - BVS] = AV
(6.4.3)
(6.4.4)
0
0
212
(6.4.5)
(6.4.6)
(6.4.7)
and
and
(6.4.8)
213
(a)
Figure 632
214
subsystem 1
'
_L
subsystem 2
subsystem 3
subsystem 4
are obtained if the voltage and current at the point of connection are
stabilised, and if each component/model is represented in the other as a
linearised equivalent around the solution at the previous time step. In the
case of synchronous machines, a suitable linearising equivalent is the
subtransient reactance, which sould be connected in shunt with the machine
current injection.
Table 6.2 Factorisation overheads as a function of subsystem size for the CIGRE
benchmark
Time step
(A*)
Factorisations
Number of nodes
CPU time
(s)
40
40
40
40
40
40
9680
1
9680
1
9677
1
46
46
92
92
184
184
150
117
194
168
294
247
215
Next Page
216
Other switch types exist which are specifically designed for statistical
overvoltage studies. EMTP provides the means for multiple runs with some
parameter variations, for example, the time delay between switch openings
and closing can be systematically incremented or randomly varied with a
given distribution during these runs.
Switches are modelled in EMTP as short circuits when conducting, and as
resistances when not conducting. This means that the number of nodes in
the system varies with the switch positions during a simulation.
Chapter 7
Electromechanical stability
7.1 Introduction
Stable operation of the power-system network is dependent on the balance
of mechanical and electromagnetic forces keeping the generators in synchronism. A disturbance on the network may or may not result in the system
falling out of synchronism and electrically collapsing. The assessment of
such behaviour is the purpose of electromechanical-stability studies.
It is normally assumed that, prior to dynamic analysis, the system is
operating in the steady state and that a load-flow solution is available.
In line with the power-flow philosophy described in Chapter 2, the
network solution of the stability studies uses the nodal-matrix method.
Thus, for each network-loading component an injected current is calculated
by solving the relevant differential and algebraic equations. The nodal
voltages are then obtained from the current injections and matrix admittance. However, as the network voltages affect the loading components, an
iterative process is often required.
Two types of stability study are normally carried out. The subsequent
recovery from a sudden large disturbance is referred to as transient stability
and the solution is obtained in the time domain. The period under
investigation can vary from a fraction of a second, when first-swing stability
is being determined, to over ten seconds when multiple-swing stability must
be examined.
The term dynamic stability is used to describe the long-term response of
a system to small disturbances or badly-set automatic control. The problem
can be solved either in the time domain or in the frequency domain. In this
Chapter both transient stability and dynamic stability are considered together in the time domain, the latter as an extension of the former.
The main consideration of the electromechanical simulation is the rotorswing stability of synchronous machines. In this respect, as compared with
the rotor long-time constants, the AC and DC-transmission systems respond
rapidly to network and load changes. The time constants associated with the
network variables are extremely small and can be neglected without significant loss of accuracy. The synchronous machine stator time constants may
also be taken as zero. The relevant differential equations for these rapidly
changing variables are transformed into algebraic equations.
276
(7.2.1)
(7.2.2)
(7.2.3)
6 = co0 + 5
(7.2.4)
0 = (D0t + 3
(7.2.5)
where S is the initial position of the rotor with respect to the synchronously
rotating reference frame.
If 6 is the acceleration of the rotor, then from eqn. 7.2.1
= lmicb
where co is the shaft angular velocity.
(7.2.6)
P. = lJ&
(7-2.7)
(7.2.8)
If the rotational power losses of the machine due to such effects as windage
and friction are ignored, then the accelerating power equals the difference
between the mechanical power and the electrical power. From eqn. 7.2.2,
for a generator
P a = Pm - P e
(7.2.9)
M = P M - P e
(7.2.10)
^ ?
co0
(7.2.11)
where Gr is the three-phase MVA rating of the machine and coo is the
synchronous speed. Inserting eqn. 7.2.11 into eqn. 7.2.10 gives
***
O)0
Pm-r.
(7.2.12)
or
(7.2.14)
278
where N p is the number of pole pairs of the machine. Eqn. 7.2.13 then
becomes
-Pe]
(7.2.16)
(7.2.17)
Se = oj-co0
E = Ed+jEq
(7.2.19)
(7.2.20)
Electromechanical stability
>
279
> q axis
IdXd
[Ef-(Xd-X'd)ld-E'll]
rd0
[ _
where
Ey = field-winding voltage
X = synchronous reactance
TQ = transient open-circuit time constant.
(7.2.23)
(7.2.24)
280
(7.2.25)
(7.2.26)
rfO
-^r^
(7.2.27)
(7 2 28)
- -
where
E" = subtransient internal voltage
X" = subtransient reactance
T'o = subtransient open-circuit time constant.
Although the synchronous-machine model is commonly represented as a
voltage source behind a fixed impedance (either the transient or subtransient reactance and the armature resistance), it is actually modelled in
transient-stability programs as a current source in parallel with a fixed
admittance.
Advanced models can, to some extent, take into account saturation of the
machine. This is achieved in the TS program by modification of the current
injection used in representing the machine.5 The single-phase, fundamental-frequency nature of the stability programs will, however, always limit the
usefulness of this modelling.
Se-f(Ef)
Specifiec
voltage
1
Controlled busbar
voltage
0+TrP)
Ka(l+Tbp)
U+Tap)
fl.
+
Fornax
1
(Ke+Tep)
EF
max
^f max
Field
voltage
\
yv
aux
Other
Kx
U+Txp)
(1+TyP)
(1+Tzp)
1
U+Tdp)
Kfp
<1+Tfp)
Typel
signal
t00
282
R(J+T2P)
(1+TIP)
Machine T +
speed
I
tf=
I
O+T3p)
Boiler or
water power
G
\ Cmax
>(n) >Pgv
Power to
turbine
Ref. speed
Power delivered
through gate
(a)
+ Mechanical power
to generator
Power delivered
through valve
(b)
Figure 7.4
Mechanical power
to generator
PI, IP and LP
turbine power
a Hydro turbine
b Thermal turbine
Basic linear models of hydro and thermal turbines are shown in Figure
7.4. The hydro turbine model includes the penstock which gives the
characteristic lead-lag response of this type of turbine.
A rigorous formulation of the equations relating to the above controllers
can be found in Reference 5.
Electromechanical stability
k
283
Network
imaginary axis
direct axis
quadrature axis
cos d
sin 5
(7.2.29)
(7.2.30)
sin 3 cos (5
(7.2.31)
284
(a).
(b)
(7-2.32)
where
A suitable value for Y o is found by using the mean of direct and quadrature
admittances, i.e.
(7.2.33)
where
Thus, the unadjusted current components are
1
Ra
-x,
Ra
(7.2.34)
Electromechanical stability
285
i( X ;-x;')
2
(Ra + x;x;')
K - v.
(7.2.35)
-sin 2^ cos 23
E;-vr
sin 23
K-vm
cos 23
(7.2.36)
= kp(Wv(f)pf
(7.3.1)
(7.3.2)
286
pv
qv
pf
qf
Fluorescent lamp
Filament lamp
Heater
Induction motor (half load)
Induction motor (full load)
Reduction furnace
1.2
1.6
2.0
0.2
0.1
1.9
3.0
0
0
1.6
0.6
2.1
-1.0
0
0
1.5
2.8
-0.5
2.8
0
0
-0.3
1.8
0
A constant-impedance load is, therefore, totally included in the networkadmittance matrix and its injected current is zero. This representation is
extremely simple to implement, causes no computational problems and
improves the accuracy of the network solution by strengthening the diagonal elements in the admittance matrix.
Nonimpedance loads may be treated similarly. In this case, the steadystate values of voltage and complex power obtained from the load flow are
used to obtain a steady-state equivalent admittance (Yo) which is included
in the network admittance matrix [Y]. During the stability run, each load is
solved sequentially along with the generators, etc. to obtain a new admittance (Y)
Y
= jvji
<7-3-3)
The current injected into the network thus represents the deviation of the
load characteristic from an impedance characteristic
Iinj=(Y0-Y)V
(7.3.4)
Electromechanical stability
287
where [Iinj] is the vector of injected currents into the network due to
generators, converters and loads, and [V] is the vector of nodal voltages.
Any network (including loads) components represented by impedances may
be directly included in the network-admittance matrix, with the injected
currents (if any) due to these components set to zero.
288
(7.5.1)
*d**d + Moad
(7.5.2)
(A + Rd)
[(3V2/7r)flVtermcosamin - Vload]
[R, + ( 3 X ]
(7.5.3)
Delay angle a
load
Electromechanical stability
289
Natural voltage
characteristic
(slope = : f x c )
Constant current
control characteristic
(slope = -A)
where Ld represents the equivalent inductance in the load circuit. Substituting for Yd using eqn. 2.2.8 gives
Vload
(7.5.5)
290
where
K' = K + 73(1 - K 2 ) 1/2
R; = (3/rc)Xc + (3/2)Re
and Re is the transformer resistance.
The additional delay (al) before the valve conducts is
al = 6 0 - a - c o s " 1 (K)
(7.5.8)
K - cos(60-a)
'
2 cos(a)
- E
- C
30
- D
2//3 cos(a-30)
- F
cos(a-60)
1//3 O cos(Yn)
- F
C
E :
90-Y
< o < 90 -
<a
K - cos(a) - cos(150-* 0 )
K - oas(a) * cos(YQ)
120
150
180
Firing Angle
(7.5.9)
(7.5.10)
where
R; = (9/TI)XC + (3/2)Re
The total delay before the valve fires is 30 degrees.
The situation in modes D and F is similar to that of mode C, where there
are periods of AC and DC short circuits, the difference being that there is
no additional delay. Mode D represents abnormal rectifier operation, and
mode F represents abnormal inverter operation. The converter equations
292
nn
(7.5.11)
(7.5.12)
(7.5.13)
cos(2oO - cos(2(5')
(7 5 14)
where
5' = OL + ii + 30
a' = a - 30
A changing mode of operation can be considered as a discontinuity, and if
during one time step the mode of operation of a terminal changes, the
solution of eqn. 7.5.5 becomes invalid. The ideal method for implementing
mode changes is to find the exact instant when the mode of operation of
the terminal changes, and adjust the time step used so that for every
solution of eqn. 7.5.5, the terminal is in the same mode of operation over
the whole time step. This method is difficult to implement, since it requires
special logic to determine exactly when a mode change occurs, and then a
number of step length changes have to be made. Instead, when a mode
change is detected at the end of a DC time step, the DC currents are kept
constant, and the discontinuity equations presented in section 7.6.2 (i.e.
eqn. 7.6.6) are formed and solved. Thus, the DC-system variables are
continuous over the mode changing boundaries and any errors introduced
using this technique are small.
A mode change can also occur when solving the discontinuity equations:
in both cases, the value Vd changes. Inspection of eqns. 7.5.7, 7.5.10 and
7.5.12 shows that in the abnormal modes of operation, the value of Rc
changes; however, in eqn. 7.5.7 there is a nonlinear relationship between
DC voltage and current, not allowing a definite value to be determined for
Rc in mode B operation. The solution adopted is to leave Rc constant
Electromechanical stability
293
(7-5.15)
where
J
" iv t e r n
and
(7 5 17)
- -
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 approaches zero. However, when the
rectifier model is modified to account for the dynamic behaviour of the DC
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, a
sequential-solution technique can exhibit convergence problems,13 especially when the voltage is low. With small terminal voltages, the AC-current
magnitude of the rectifier load is related to the DC 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 currentphase angles.
Good convergence is achieved with a unified algorithm11 by reducing the
AC 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
294
LoJ
and
% = 1 + JO
(7.5.19)
The equivalent circuit shown in Figure 7.10 can now be applied to find the
rectifier current (Ip) by using the Newton-Raphson technique.
The effect of the rectifier on the rest of the system can be determined by
superposition
[V] =
(7.5.20)
where
fVff] = [Y]- 1 [ir nj ]
(7.5.21)
and [Ifnj] are the injected currents resulting from all other generation and
loads in the system.
If the network remains constant, vector [Z] is also constant and thus only
needs re-evaluation on the occurrence of a discontinuity.
The equivalent system of Figure 7.10 contains seven variables (Vterm, lp,
d, \j/9 a, Vd and Id). With these variables, four independent equations can be
formed. They are eqn. 2.2.8 and
(7.5.22)
V,. I , - V3aV t e r m . I p cos(0 - >/0 = 0
(7.5.23)
Eqn. 7.5.22 is complex and represents two equations. Substituting for V,,
and I p using eqns. 2.2.16 and 7.5.1 reduces the number of variables to five.
A fifth equation is necessary and with constant-current control, i.e. with the
Electromechanical stability
295
Dynamic
load element
"load
(7.5.25)
When the delay angle reaches a specified lower limit (oemin), the control
specification, given by eqn. 7.5.24 changes to
=0
(7.5.26)
Eqn. 2.2.8 is no longer valid. The DC current (Id) is now governed by the
differential eqn. 7.5.5. If the trapezoidal method is being used, this equation
can be transformed into an algebraic form and eqn. 2.2.8 is replaced by
ld = ka. E r . cos a kb = 0
(7.5.27)
(7.5.28)
S /2
2 k/i V
^ n a.ka.yteim(t)cosoi(t) + ' ' 'oad
(7.5.29)
296
where
kc = ( ^ )
+ 1/Tdc
(7.5.30)
and represents the time at the beginning of the integration step and h is
the step length. Commutation angle \i is not explicitly included in the
formulation, and since these equations are for normal operation, the value
of k in eqn. 2.2.16 is close to unity and may be considered constant at each
step without loss of accuracy. On convergence, fi may be calculated and a
new k evaluated suitable for the next step.
7.5.2 DC links
For as long as the DC link is operating normally, i.e. without commutation
failures, the converter steady-state equations can be used to simulate its
behaviour at each step of the stability study.
The initial steady-state operating condition of the DC link is obtained
from the power-flow solution, with which the control type and mode,
current and margin settings will have been established.
At each iteration of the stability study the control mode must be reconsidered. If the link operates on current control this can be done by first
assuming mode 1 (i.e. with the rectifier on c.c. control), and by combining
eqns. 2.2.19 and 3.3.6 a DC current can be determined as
dmodel
[1 + (Rd -
(mXJ/A}
Assuming this current to be valid, then DC voltages at each end of the link
can be calculated using eqns. 2.2.8 and 2.2.17. The DC link is operating in
mode 2 (i.e. with the inverter on c.c. control) if
V^ex-V^^O
(7.5.32)
[(1
For constant power control, under control mode 1, the DC current may be
determined from the quadratic equation
(7.5.34)
Electromechanical stability
297
outside
within
within/
greater
less
greater
less
less
less
where Pds is the setting at the electrical mid point of the DC system, i.e.
P * = (Ps, + P y 2
(7.5.35)
The correct value for Idmodei can then be found from Table 7.2. Control
mode 2 is determined using eqn. 7.5.33 and in this case the following
quadratic equation must be solved
r- *dmode2
^0**dmode2
* dmarg ~~ *ds ~
"
(7.5.36)
where
+
_V2
(7.5.37)
(7.5.38)
(7.5.39)
298
It is also possible that when the link is operating close to the changeover
between modes, convergence problems will occur in which the control mode
changes at each iteration. This can easily be overcome by retaining mode 2
operation whenever detected for the remaining iterations in that particular
time step.
The value of the DC current obtained from the above equations can
then be used to calculate the inverter and rectifier DC voltages and also the
AC-line currents. The remaining variables for the rectifier and inverter ends
can be found from a unified iterative solution of the set of equations
containing these variables, such as shown in eqn. 7.5.25 for the case of a
rectifier on constant-current control.
The active and reactive power at both converters can then be calculated
with
(7.5.40)
(7.5.41)
(7.5.42)
r
(7.5.43)
t a n (/>j
DC power modulation
It has been shown in the previous section that under the constant-power
control mode, the DC link is not responsive to AC-system terminal conditions, i.e. the DC-power transfer can be controlled disregarding the actual
AC-voltage angles. Since, generally, the stability limit of an AC line is lower
than its thermal limit, the former can be increased in systems involving DC
links by proper use of fast-converter controllability.
The DC power can be modulated in response to AC-system variables to
increase system damping. Optimum performance can be achieved by
controlling the DC system so as to maximise the responses of the AC system
and DC line simultaneously following the variation of terminal conditions.
1 1
AC system AL DC system
controller
controller
ai
DC network
DC control modes
AC/DC interface
0)
(iii)
Electromechanical stability
299
The AC-system controller uses AC and/or DC-system information to derive the current and voltage-modulation signals. A block
diagram of the controller and AC-DC signal conditioner is shown in
Figure 7.12.
(ii) The DC-system controller receives the modulation signals AI and AE
and the steady-state specifications for power P o , current I o and voltage
Eo. Figure 7.13a 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 DC link.
The rectifier current controller, Figure 7.13fc, includes signal limits
and rate limits, transducer-time constant, bandpass filtering and a
voltage-dependent current-order limit (VDCOL).
The inverter-current controller, Figure 7.13^, includes similar components plus a communications delay and the system-margin current
(ijFinally, the DC-voltage controller, including voltage restart dynamics, is illustrated in Figure 7.13^.
(iii) The DC current Id and voltage Ed derived in the DC-system controller
constitute the input signals for the AC-DC network model which
involves the steady-state solution of the DC system (neglecting the
DC-line dynamics which are included in the DC-system controller).
Here the actual AC and DC-system quantities are calculated, i.e.
control angles, DC current, voltage, active and reactive power. The
converter AC-system constraints are the open-circuit secondary voltages Ear and Eai.
AC/DC
signals
r3~^
^k i
AC system
controller
>
PTd
>
PTd
PTd
pTd
-A1,
}*
A/;
-AEr
300
dr
Current order
restart dynamics
\
=' ^1 ^ o ,
(a)
Current rate
limit
I
02
Current
limit
hi
Regulator
b4
dr
D(P)
1+pT
OmaxI
(b)
Current rate
limit
Communication
delay
Current
limit
k>5
(c)
E
0J
02
Regulator
l
Voltage rate
limit
04
(d)
N(P)
D(P)
O7
N(P)
D(P)
Electromechanical stability
301
302
start
J)
Electromechanical stability
main program
303
converter subroutine
solve machines
fort!]*
solve [V] = [ Y]-1 [ l]m
without converter
yes
obtain new
Thevenin impedance
vector
extract Thevenin
voltages from [V]
_L
I calculate mismatches
yes
convergence
calculate Jacobian for
given mode of
operation
>0
yes
block converter
I update variables
test mode of control
and mode of operation
calculate converter
variables
iter >max
yes
error message
"
(7.6.1)
304
y(t)
(1+pT)
iterative solution. However, the solution can be made direct by incorporating the differential equations into eqn. 7.6.1.
For example, consider the trivial transfer function shown in Figure 7.16.
The differential equation for this system is given by
py(t) = (G.z(0 ~y(t))/T
(7.6.2)
with the input variable being denoted by z to indicate that it may be either
integrable or nonintegrable.
The algebraic form of eqn. 7.6.2 has a solution at the end of the (n + l)th
step of
yn + 1 =cn+1
+mn + 1.zn+1
(7.6.3)
where
Cn+i
= (1 - 2bn+1)yn + bn+1.G.zn
m*+i=K+i-G
(7.6.4)
(7.6.5)
and
hH+l)
(7.6.6)
n+1
(7.6.9)
(7.6.10)
306
and fed as part op the input data of the TS program. However, for all normal
controller configurations tested, the algorithm correctly determined the
required current/power settings.
After a discontinuity, such as a fault application or switching on the AC
system, good starting conditions have to be determined to ensure that
the DC equations converge to the correct solution at the next time step.
The DC network contains only passive components, and the network can
experience an instantaneous drop in DC voltage. However, the DC-network
inductance prevents the network DC currents from changing instantly.
Immediately following a discontinuity, new DC voltages have to be determined for use as starting values in the DC system solution at the next time
step.
By summing the rate of change of DC current at every node in the DC
network, and setting the total change to zero, the following equation is
derived for the DC network at a discontinuity
[V] = [H]- 1 tC]
(7.6.11)
where
V
dsi
- 1/KC
K )
,*, g , gx
' " '
Electromechanical stability
307
and the terminal voltage is less than Vdsr, the mode choice is incorrect. If
the rectifier is assumed to be in the constant-current control, CCC, mode
and the DC voltage is greater than Vdsr the mode choice is incorrect. For an
inverter, if the DC voltage is greater than Vdsi the mode choice is incorrect.
For any valid set of AC voltages, the following algorithm is used to
determine the DC voltages of the HVDC link:
(i)
308
(7.6.14)
for inverter
h = - l
d s
(7.6.15)
1zVd
r ^7"
x;
sx:
for inverter
d
^7
^/.O.l I)
(7.6.19)
Qdc = Pdctan((t>)
(7.6.20)
Bs = -Qoc/VLn
(7-6.21)
(7.6.22)
Electromechanical stability
309
erm-Bs
(7.6.23)
(7.6.24)
I , = -V r e G + VimB
(7.6.25)
I6=-VimG-VreB
(7.6.26)
310
AV = 1 -
Id Xcpu + (ldFL/ld)[cos(y0
+ </>)- cosy]
or
*cpu ^
AV 1 - ^
x
h cvu + cos(j0 + </>)- cosy
(7.6.27)
where
AV = inverter sudden commutation-voltage reduction required to produce the theoretical onset of commutation failure; it refers to the
phase voltage for single-phase faults, specifically on the faulted
phase
for symmetrical 30 reductions, the per unit AV is the same for
phase and line quantities
Id = DC current at the beginning of the commutation
Id = DC current at the end of commutation
Xcpw = commutation reactance (in per unit), i.e.
x
=X =
b
V2W.
FL
Electromechanical stability
311
The two dependent factors AV and <j> can be derived from an iterative
solution of eqns. 7.6.27 and 7.6.28.
Prior to control action by the converters, any DC current increase caused
by the inverter AC-voltage reduction (i.e. the factor iyi^) is critical for
commutation failures. The current increase depends on the values of the
smoothing reactor, commutation reactance, line parameters and DC-side
filters.
For the short period between fault initiation and the first commutation
failure, the DC-current rise is relatively small and very nearly linear.
Extensive EMTDC simulations have been carried out on two schemes with
very different DC-side and commutation parameters to assess the range of
variation of iyi d to be expected.
In the New Zealand scheme, that ratio varied from 1.1723 to 1.1457 for
7=18 (the nominal value in current use) to 1.238 for y = 24, the critical
commutation margin being y0 = 8.
Substitution of these values in eqn. 7.6.27 results in voltage reductions of
between 12 and 14% for the case of a symmetrical fault.
Due to the zero-crossing phase shift caused by the voltage unbalance, the
probability of commutation failure is higher in the case of a single phase-toground as compared with the symmetrical faults.
However, the difference between the 3$ and single phase is, to a large
extent, counteracted by the larger DC-current increase caused by the
symmetrical fault. The combination of voltage magnitude and phase shift
following a single-phase fault reduces the voltage area available to the
commutation; this effect is more pronounced in the case of the starconnected transformer as compared to the delta connection. It is, therefore,
expected that the converter group with the star connection may be more
prone to commutation failures. As the standard HVDC converter contains
star and delta-connected groups, there is no need to consider the delta
connection when assessing the onset of commutation failure.
Thus, at every step of the stability study, the inverter voltages are checked
against the critical values to decide whether the quasi steady state model is
still applicable. If these voltages are below the critical levels, indicating the
onset of commutation failure, the prediction of DC-link performance needs
transient simulation, as described in the next Chapter.
312
AC-DC
0.0365
0.5968
2.5
2.5
345kV2ll.42kV
1196 MVA
n HMH
0.0365
MW
2H.42kV23OkV
1196 MVA
26.0
To
Rectifier
ac
system
0.5968
0.7406
warn
15.04
167.2
|3.23
0.0606
0.935
Electromechanical stability
TXvizcl
Clyde
Benmore
-HH*
Livingston
Aviemore
Table 73
Bus
ivu
Aviemore
Benmore
Clyde
Roxburgh
Livingston
Twizel
1.015
1.000
1.043
1.042
1.024
1.015
1.338
0.000
6.479
6.338
2.415
1.564
VDC
"DC
a
PDC
QDC
Rectifier
Inverter
500 kV
2kA
17.65
1000MW
583MVAr
490 kV
2kA
15
980 MW
545 MVAr
313
314
Electromechanical-stability programs are somewhat limited in their simulation of abnormal operating conditions. As described in section 7.5,
an abnormal operating condition exists when the commutation angle
exceeds 60 degrees. If this occurs at the inverter terminal, the inverter is
simply bypassed. The rectifier terminal is kept under operation and
maintains the current to at least 30% of its nominal setting. After a set
period, if the problem causing abnormal operation has been removed, the
inverter is restarted and the power ramped back up to normal using current
control.
The test system in the steady state has the rectifier AC-system bus-voltage
magnitudes and angles as shown in Table 7.3. The corresponding DC
steady-state conditions are shown in Table 7.4.
Roxburgh
-Aviemore
Benmore
Clyde
1.10
1.05&
1.00
"*...-'
^-
0.90
(a)
0.0
1.0
2.0
3.0
4.0
5.0
Time (s)
(b)
Figure 7.19
6.0
7.0
8.0
9.0
10.0
Electromechanical stability
315
2.1 -
Field >/olta
2.5 -
l.o -
Benmore
Clyde
1.7-
4/
0.90.5-
(a)
-Roxburgh
Aviemore
Benmore
o.o
1.0
2.0
3.0
4.0
5.0
6.0
Time (s)
(b)
Figure 7.20 Rectifier AC-system generator parameters
7.0
8.0
9.0
10.0
316
Roxburgh
Aviemore
(a)
- Roxburgh
Aviemore
Benmore
55
35
I 15
PQ
0.0
1.0
2.0
3.0
4.0
5.0
Time (s)
(b)
Figure 7.21 Rectifier AC voltage
6.0
7.0
8.0
9.0
10.0
Electromechanical stability
317
7.00
6.00
jg 5.00
i l 4.00
>
3.00-
2.00
1.00-
0.00
(a)
Rectifier alpha
(b)
- Rectifier P
- Rectifier Q
1200
0.0
1.0
2.0
3.0
4.0
5.0
Time (s)
(c)
Figure 7.22
Rectifier DC variables
6.0
7.0
8.0
9.0
10.0
318
7.00
6.00
8
40
>
3.00
2.00
5.00
I
1.00
0.00
(a)
Inverter gamma
100
60
40
20
0
(b)
- Inverter P
Inverter Q
1200
10.0
Electromechanical stability
-Roxburgh
Benmore
- Aviemore
Clyde
319
-1.0
(a)
1.0
2.0
3.0
4.0
5.0
6.0
7.0
8.0
9.0
10.0
Time (s)
(b)
and angles (angles relative to the reference bus Clyde). The effect at the
fault terminal is particularly evident as the voltage magnitude drops to zero
and the voltage angle rises nearly 50 degrees. The return to a stable voltage
angle can also be seen, with an initial oscillation due to the generator
controllers eventually decaying away.
Figures 7.22 and 7.23 display the effect of the fault on the DC variables
at each terminal of the link. At the instant of the fault, the rectifier terminal
is blocked and the inverter terminal bypassed. Once the fault is removed,
the rectifier terminal is unblocked and the DC current through the rectifier
ramped to 30% of its nominal value. The inverter terminal is then unbypassed and the DC real power ramped back up to its nominal value.
The machine response to the fault is shown in Figure 7.24. The field
voltages and machine angles vary and oscillate through the response of the
AVR but eventually settle back to steady-state values. The Aviemore
320
field-voltage response is somewhat different than that for the other generators owing to its larger forward regulator gain and lower upper-field
voltage limit.
7.8 Summary
This Chapter contains a concise description of the synchronous-generator
transient mechanical and electrical response to disturbances in the power
system. The modelling of other system components is also discussed with
emphasis on the AC-DC converters.
The structure of a basic transient-stability program is then described with
representation of the DC-link behaviour on the assumption that the link
maintains continuous controllability during the disturbance; however, a
check for the onset bf commutation failure has been included as part of the
solution. Prediction of such an event would indicate the limit of applicability
of the algorithm and the need to adopt the more detailed stability simulation discussed in Chapter 8.
7.9 References
1 KIMBARK, E.W.: 'Power system stabilityvol. Ill' (Synchronous Machines,
Dover Publications Inc., 1968)
2 CLARKE, E.: 'Circuit analysis of a-c power systemsvol. IF (John Wiley 8c Sons
Inc., New York, 1950)
3 PARK, R.H.: 'Two reaction theory of synchronous machinespart I: generalized
method of analysis', IEEE Trans., July 1929, 48, pp. 716-730
4 PARK, R.H.: 'Two reaction theory of synchronous machinespart IF, IEEE
Trans., June 1933, 52, pp. 352-355
5 ARRILLAGA, J., ARNOLD, C.P., and HARKER, B.J.: 'Computer modelling of
electrical power systems' (John Wiley 8c Sons, Chichester, England, 1983)
6 IEEE committee report: 'Computer representation of exciter systems', IEEE
Trans., June 1968, PAS-87, (6), pp. 1460-1464
7 IEEE committee report: 'Dynamic models for steam and hydro turbines in
power-system studies', IEEE Trans., 1973, PAS-92, (6), pp. 1904-1915
8 ARNOLD, C.P.: 'Solution of the multimachine power system stability problem'.
PhD thesis, UMIST, Manchester, UK, 1976
9 DANDENO, P.L., and KUNDUR, P.: 'A noniterative transient stability program
including the effects of variable load-voltage characteristics', IEEE Trans., 1973,
PAS-92, (5), pp. 1478-1484
10 BERG, G.L.: 'Power system load representation', Proc. IEE, March 1973, 120, (3),
pp. 344-348
11 ARNOLD, C.P., TURNER, K.S., and ARRILAGA, J.: 'Modelling rectifier loads
for a multi-machine transient stability programme', IEEE Trans., 1980, PAS-99,
pp. 78-85
12 GIESNER, D.B., and ARRILLAGA, J.: 'Operating modes of the three-phase
bridge converter', Int.]. Elect. Eng. Educ, 1970, 3, pp. 373-388
13 DOMMEL, H.W., and SATO, N.: 'Fast transient stability solutions', IEEE Trans.,
1972, PAS-91, (4), pp. 1643-1650
Electromechanical stability
321
Chapter 8
8.1 Introduction
The previous two chapters have described independently the electromagnetic and electromechanical behaviour of AC-DC power systems. In practice, however, the two are closely interrelated and the results obtained for
either of them independently will be in error.
It is, of course, possible to include the equations of motion of the
generators in the electromagnetic-transient program to provide more realistic solutions. However, considering the different time constants influencing the electromechanical and electromagnetic behaviour, such an approach
would be extremely inefficient. The electromagnetic-transient simulations
use steps of (typically) 50 //s, whereas the stability programmes use steps at
least 200 times larger.
To reduce the computational requirements the NETOMAC package1 has
two separate modes. An instantaneous mode is used to model components
in three-phase detail with small time steps in a similar way to the EMTP/
EMTDC programs.2 The alternative is a stability mode and uses r.m.s.
quantities at fundamental frequency only, with increased time-step lengths.
The program can switch between the two modes as required while running.
The HVDC converter is either modelled elementally by resistive, inductive
and capacitive components, or by quasisteady-state equations, depending on
the simulation mode. In either mode, however, the entire system must be
modelled in the same way. When it is necessary to run in instantaneous
mode, a system of any substantial size would still be very computationally
intensive.
A more efficient algorithm that combines the two programs described in
Chapters 6 and 7 is described here. 3 ' 4 It is a hybrid algorithm which takes
advantage of the computationally inexpensive dynamic representation of
the AC system in a stability program, and the accurate dynamic modelling
of the converter nonlinearities.
The slow dynamics of the AC system are sufficiently represented by the
stability program and, at the same time, the fast dynamic response of the
324
Detailed
analysis of dc
system (and possibly
some of ac system)
Figure 8.1
325
Stability Program
(a)
EMTDC Program
Stability Program
Stability
Program
(b)
With reference to Figure 8.2a, initially, the TSE hybrid reads in the data
files and runs the entire network in the stability program, until electromechanical steady-state equilibrium is reached. The quasisteady-state
representation of the converter is perfectly adequate as no fault or disturbance has yet been applied. At a selectable point in time prior to a network
disturbance occurring, the TS network is split up into the two independent
and isolated systems, system 1 and system 2.
For the sake of clarity, system 1 is classified as the AC part of the system
modelled by the stability program TS, and system 2 is the part of the system
modelled in detail by EMTDC.
The snapshot data file is now used to initialise the EMTDC program used,
instead of the TS representation of system 2. The two programs are then
326
8.2.2 Dataflow
Data for the detailed EMTDC model is entered into the program database
via the PSCAD graphics. Equivalent circuits are used at each interface point
to represent the rest of the system not included in the detailed model. This
system is then run until steady state is reached and a snapshot taken. The
snapshot holds all the relevant data for the components at that point in time
and can be used as the starting point when interfacing the detailed model
with the stability program.
The stability program is initialised conventionally through power-flow
results via a data file. An interface data file is also read by the TSE hybrid
and contains information such as the number and location of interface
buses, analysis options and timing information. The data flow diagram is
shown in Figure 8.4.
J_
yes
pass information
to detailed model
J_
solve EMTDC
extract information
from detailed solution
and include in stability
programme
T = T + step length
i
output results
327
328
positive
sequence
data for
whole
system
Figure 8.4
output
Dataflow
Figure 8.5
Hybrid interface
329
Interface
System 1
System 2
shown in Figure 8.6, where E t and Zt can represent the equivalent circuit
of system 1 and Ic and Z2 the equivalent circuit of system 2.
330
Ei
(a)
(b)
Figure 8.7
(8.3.2)
IF-
(8.3.3)
331
833
(8.3.4)
where
Yps = positive sequence voltage
Vfl, Vfc, Vc = phase voltage
a = 120 degree forward rotation vector (i.e. a = 1 L120)
Positive-sequence data from the stability program is converted to three
332
phase through simple multiplication of the rotation vector. For the voltage
Va = Vps
(8.3.5)
% = a2Vps
(8.3.6)
Vc = aVps
(8.3.7)
* ps '
2
\. =a V
o
ns '
zs
yo.o.oj
+ aV + V
ps
'
ns
V, = aVnv 4- d2Y
c
ps
'
zs
+V
ns '
zs
(8 3 Q\
\**J,KJ*%JJ
f8 S 101
\^ *-' y*j
where
Vns = negative-sequence voltage
Vzs = zero-sequence voltage
Similarly, any unbalance in system 2 can be accommodated in the transientstability program.
(8.3.11)
V = I2Z2
(8.3.12)
h = h +!c
(8.3.13)
333
(8.3.14)
(8.3.15)
and
0n = ey-(j)
(8.3.16)
where 0 is the displacement angle between the voltage and the current.
Thus, eqn. 8.3.15 can be written as
x
(8.3.17)
where
= 0Zi ~
If
Ei = E l r
then equating the real terms
E l r = ( 1 ^ ! cosj? + V)cos0v -f ( - 1 ^ ! sin jS)sin0v
(8.3.18)
(8.3.19)
B= -I^sinjS
(8.3.20)
334
(8.3.21)
where
the voltage angle 6y in the TS phase reference frame can be calculated, i.e.
(8.3.23)
(8.3.24)
- -
(8.3.26)
v
D = cos 9Z2 - lx sin <j)
(8.3.27)
(8.3.28)
(8.3.29)
335
Knowing the EMTDC voltage angle, 6y, allows calculation of the EMTDC
current angle, 6n, from eqn. 8.3.28. The magnitude value of Ex can then be
derived from eqn. 8.3.11.
The r.m.s. power can be extracted directly from the waveforms and,
therefore, Fourier transform or curve fitting methods are not necessary, i.e.
<8A2)
This method greatly reduces the computing time necessary for the data
extraction from EMTDC.
The choice of r.m.s. power was made on the basis that harmonic power
flow will result only from in-phase components of harmonic voltage and
current. The assumption was made that if a system contains only a lowresistive component, then the harmonic power flow is not significant.
This, however, is not valid for every situation and, particularly, at the
inverter end of an HVDC link, the effect of the resistive component of the
network is not insignificant. Certain harmonic frequencies in a network may
also be parallel resonant or close to parallel resonance and exhibit more
resistance than reactance. The presence of a transient may excite the
resonant frequency and greatly affect the results. It is important, therefore,
to model.the fundamental power flow accurately.
Another factor to consider is that the direction of fundamental power may
not necessarily be the same as the direction of harmonic power. With an
HVDC link, although fundamental power is drawn into a rectifier, much of
the harmonic power flow will be in the reverse direction. The r.m.s. real
336
800 ::
600 -
f\, ; \?V V? V
'
0.20
0.28
400
200 -
ti ; ' y
0 -200 -400-600
-800 10000.00
0.04
0.08
0.12
0.16
0.24
0.32
0.36
0.40
Time (s)
Figure 8.8
power measured is the difference between these two powers and so not
entirely representative of the fundamental power load of the HVDC link.
Conversely, at the inverter end, the r.m.s. real power will include harmonic
power flow into the AC system. This may exaggerate the amount of true
fundamental power from the inverter.
Figure 8.8 shows an EMTDC-simulated result of a single-phase fault at the
inverter end of the CIGRE HVDC benchmark model (see Appendix III).
The instantaneous power is shown, along with the total r.m.s. power over
discrete one-cycle intervals. For hybrid interfacing purposes, data is transferred at discrete intervals equal to the stability program time step, and
these intervals are significantly larger than the interval between the
discrete output points constituting the instantaneous power. The r.m.s.
power is then taken over a fundamental period to represent its use in a
hybrid situation.
Figure 8.9 shows an analysis of the single-phase fault comparing fundamental-frequency power with the total r.m.s. power. The fundamentalfrequency power was derived using the curve-fitting method described in
the next section to extract both fundamental voltage and current. The
comparison shows that, particularly during the fault time, there exists a
significant amount of Fh or harmonic r.m.s. power in the total r.m.s. power.
When the three-phase network is unbalanced the fundamental-frequency
337
900 800 -
700 -
'
'
0.28
0.32
600500
400
300
200
100
0-100
-2000.00
0.04
0.08
0.12
0.16
0.20
0.24
0.36
0.40
Time (s)
Figure 8.9 Fundamental power versus total r.m.s. power
(8.4.3)
and the negative and zero-sequence powers cause additional power loss in
the network. Figure 8.10 shows the sequence components of the fundamental r.m.s. power shown in Figure 8.9.
Fundamental-frequency negative-sequence currents, in the presence of
damper windings, can produce a braking torque which will retard the
rotor.8 Damper windings, however, also serve to lower the negative-sequence impedance of a machine which, in turn, reduces the negativesequence voltage.9 Which of these two opposing effects is dominant depends on the resistance of the damper windings. High-resistance windings
cause the braking torque to be the significant effect. The braking power of
the negative-sequence current is
Pb = \ I^R,
I2 (r - r)
(8.4.4)
(8.4.5)
338
000 -I
900
800 -
700
'
600
500
400
300 200 100
0.04
0.08
P-'
-1.
0.12
0.16
--,
0.20
0.24
0.28
0.32
0.36
0.40
Time(s)
Figure 8.10
where
lrn = negative-sequence rotor current
R,. = rotor resistance
Isn = negative-sequence stator (or armature) current
rn = negative-sequence resistance of the machine
rp = positive-sequence resistance of the machine
The negative-sequence resistance can be approximated from the rotor and
the armature resistance, i.e.
r ^(R s + |R r )
(8.4.6)
^ I (r __ r )
(8.4.7)
(8.4.8)
339
Fundamental power
1000
-200
0.00
0.04
0.08
0.12
0.16
0.20
0.24
0.28
0.32
0.36
0.40
Time (s)
Figure 8.11
where
I s = the effective very low-frequency value of the armature current
iDC = instantaneous DC component of armature current
The total r.m.s. power, then, is not always equivalent to either the
fundamental-frequency power or the fundamental-frequency positive-sequence power. A comparison of these three powers is shown in Figure 8.1 L
The difference between total r.m.s. power and positive-sequence power can
be seen to be highly significant during the fault.
The most appropriate power to transfer from EMTDC to TS is, then, the
fundamental-frequency positive-sequence power. This, however, requires
knowledge of both fundamental-frequency positive-sequence voltage and
fundamental-frequency positive-sequence current. These two variables contain all the relevant information and, hence, the use of any other power
variable to transfer information becomes unnecessary.
340
341
2.00
-5.00
0.00
0.02
0.04
0.06
0.08
0.10
0.12
0.14
0.16
0.18
0.20
Time (s)
Figure 8.13
342
ii
_____
(a)
- Staggered CFA
1.40
1.20
1.00
0.80
0.60
0.40
0.20
0.00
0.00
0.02
0.04
0.06
0.08
0.10
0.12
0.14
0.16
0.18
0.20
Time (s)
(b)
Figure 8.14
343
Staggered CFA
Weighted CFA (continuous)
30
20
10
i
i
i
-10
!* "
i
i
i
-20
-30
-40
-60
0.00
0.02
0.04
0.06
0.08
0.10
0.12
0.14
0.16
0.18
0.20
Time(s)
Figure 8.15
344
EMTDC
EMTDC
EMTDC
Time
Figure 8.16
Interfacing methods
345
Stability
step length
EMTDC
Dynamic program
step length
Figure 8.17
346
Stability
step length
EMTDC
Dynamic program
step length
Figure 8.18
is now also applied to the TS program which is then solved for a period until
it has again reached EMTDC's position in time (arrow 4). The normal
interaction protocol is then followed until any other discontinuity is
reached.
A full-period analysis after the fault has been applied is necessary to
extract accurately the fundamental-frequency component of the interface
variables. The mechanically-controlled nature of the AC system implies a
dynamically-slow response to any disturbance and so, for this reason, it is
considered acceptable to run EMTDC for a full period without updating the
system 1 equivalent circuit during this time.
347
Although the above concept has some advantages, it also suffers from
many disadvantages. The concept is proposed, in particular, for weak AC
systems. A weak AC system, however, is likely to have any major generation
capability far removed from the converter-terminal bus as local generation
serves to enhance system strength. If the generation is, indeed, far removed
out into the AC system, then the distance required for an interface location
to achieve considerably less phase imbalance and waveform distortion is also
likely to be significant.
The primary advantage of a hybrid solution is in accurately providing the
DC dynamic response to a transient-stability program, and in efficiently
representing the dynamic response of a considerably-sized AC system to the
DC solution. Extending the interface some distance into the AC system,
where the effects of a system disturbance are almost negligible, diminishes
the hybrid advantage. If a sizeable portion of the AC system requires
modelling in detail before interfacing to a transient-stability program can
occur, then one might question the use of a hybrid solution at all and
instead use the more conventional approach of a detailed solution with AC
equivalent circuits at the system cut-off points.
Another significant disadvantage in an extended interface is that AC
systems may well be heavily interconnected. The further into the system that
an interface is moved, the greater the number of interface locations
required. The hybrid interfacing complexity is thus increased and the
computational efficiency of the hybrid solution decreased. The requirement
for a detailed representation of a significant portion of the AC system serves
to decrease this efficiency, as does the increased amount of processing
required for variable extraction at each interface location.
The advantages of using the converter bus are:
The major drawback of the detailed solution is in not seeing a true picture
of the AC system, since the equivalent circuit is fundamental-frequency
Figure 8.19
Test system
348
PhaseB
PhaseC
(a) EMTDC
5
o
-1.4
0.98
1.00
1.02
1.04
1.06
1.08
1.10
1.12
1.14
Time (s)
(c) TSE Hybrid - Frequency dependent equivalent
Figure 8.20
Inverter AC voltage
1.16
1.18
349
350
-Phase A
PhaseC
(a) EMTDC
-50.0
0.98
1.00
1.02
1.04
1.06
1.08
1.10
1.12
1.14
Time (s)
(c) TSE Hybrid - Frequency dependent equivalent
Figure 8.21
Inverter AC current
1.18
1.18
351
352
compare variables of
EMTDC with TS
Ts=Ts+stepl
yes
isolate TS systems
land 2
set up equivalent
circuits
yes
state=l
Ts=0
Ts=Ts+stepl
call EMTDC for a
period T
convert EMTDC o/p
to positive sequence
state=-l
Ts=0
Ts=Ts+stepl
Figure 8.22
353
B
passTS system 1
variables to EMTDC
passTS system 1
variables to EMTDC
passTS system 1
variables to TS sys.2,
passTS system 1
variables to TS sys.2
pass EMTDC
variables to TS sys.l
Ts=0
Ts=Ts+stepl
state=3
Ts=0
Ts=Ts+stepl
yes
extract TS system
2 variables
compare variables of
TS sys.2 and EMTDC
terminate EMTDC
reconnect TS
systems 1 and 2
354
pass TS system 1
variables to EMTDC
pass TS system 1
variables to TS sys.2
call EMTDC for a
period of T/2
convert EMTDC o/p
to positive sequence
pass EMTDC variables
toTSsys.l
state =1
Ts=0
Ts=Ts+stepl
Figure 8.22
Table 8.1
State
-1
0
1
2
3
(Continued)
Interfacing states
Event
isolate TS network ready to interface with EMTDC
initial time through interfacing routine
normal interaction between TS and EMTDC
system network change signalledrun EMTDC for a full period past
the fault and return information to TS
TS catching up with EMTDC after a system-network change
355
Equivalent circuit
Actual system
25002000 -
1000-
1
i
1500-
5000-
100 200 300 400 500 600 700 800 900 1000 1100 1200
Frequency (Hz)
Figure 8.24 Frequency spectrum of the equivalent circuit in Figure 8.23
356
357
TSonly
2.50
-2.60
(a)
- EMTDC only
TS only
(b)
-TSonly
EMTDC only
1.20
1.00
3
0.80
0.60
0.20
0.00
1.50
1.65
1.80
1.95
2.10
2.25
2.40
2.55
2.70
2.85
Time (s)
(c)
Figure 8.25
3.00
358
TSonly
1.40
1.20 1.00
0.80
0.60
0.40
0.20
0.00
(a)
1 pu reference
- EMTDC only
-2.50
(b)
- T S E (EMTDC variable)
1 pu reference
2.50
1.50
1.50
Figure 8.26
1.65
1.80
1.95
2.10
2.25
2.40
2.55
2.70
2.85
3.00
359
TSonty
TSE (EMTDC variable)
DCci
1 00
0.800.60 0.400.200.001.50
1.65
1.80
1.95
2.10
2.25
2.40
2.55
2.70
2.85
3.00
Time (s)
-EMTDConly
TSonly
10.00
(a)
TSonly
- T S E (EMTDC variable)
10.00
-2.50
1.50
1.65
1.80
1.95
2.10
2.25
2.40
Time (s)
(b)
2.55
2.70
2.85
3.00
360
-TSE
1200
1.50
Figure 8.29
1.65
1.80
1.95
2.10
2.25
2.40
2.55
2.70
2.85
3.00
The rectifier DC currents, displayed for the three solutions in Figure 8.27,
show a very similar variation for the TSE and EMTDC solutions, except for
the region between t = 2.03 s and t = 2.14 s, but the difference with the
TS-only solution is very large. The corresponding DC-voltage responses are
shown in Figure 8.28.
Figure 8.29 compares the fundamental positive-sequence real and reactive powers across the converter interface for the TS and TSE solutions.
The main differences in real power occur during the link-power ramp.
The difference is almost a direct relation to the DC-current difference
between TS and TSE shown in Figure 8.26. The oscillation in DC voltage
and current as the rectifier terminal is deblocked is also evident.
As for the reactive power Q, prior to the fault, a small amount is flowing
into the system owing to a surplus MVAr at the converter terminal. The fault
reduces this power flow to zero. When the fault is removed and the AC
voltage overshoots in TSE, the reactive MVAr also overshoots in TSE and,
361
Avicmore
Benmore
0.0
1.0
2.0
3.0
4.0
5.0
6.0
7.0
8.0
9.0
10.0
8.9 Summary
The transient-stability (TS) algorithm described in Chapter 7 has been
modified to include a detailed simulation of the HVDC converter behaviour;
the latter is obtained with the EMTDC program. Much of Chapter 8
describes the interface of the TS and EMTDC programs, which has been
made sufficiently general for use with alternative versions of the TS and
EMTP solutions.
To demonstrate the advantages of the hybrid solution, two different
disturbances have been applied to a simple test system and the results
compared with those of the conventional electromechanical-stability and
electromagnetic-transient approaches.
First, a minor disturbance has been used to validate the hybrid technique
by showing that the response is practically identical to that of the electromechanical-stability program TS.
The major disturbance, a three-phase fault at the rectifier-converter
terminal, exemplifies the differences between the two types of solution. The
electromechanical solution shows the slower dynamic response of the AC
system, and the electromagnetic solution displays the fast dynamics of the
rapidly switched converter. In the hybrid method these responses are
362
combined, thus giving a more realistic overall picture of the entire system
performance. Electromechanical modelling of the AC system has been
improved by inclusion of the actual response of the fast-acting HVDC link,
and the link behaviour has been made more realistic through the inclusion
of the generator-controller response varying the AC-system voltage at the
converter terminal. Neither program on its own can provide a realistic
response of both AC and DC systems.
8.10 References
1 KULICKE, B.: 'Netomac digital program for simulating electromechanical and
electromagnetic transient phenomena in ac systems'. Siemens Aktienngesellschaft, E15/1722-101
2 WOODFORD, D.A.: Validation of digital simulation of dc links', IEEE Trans.,
September 1985, PAS-104, (9), pp. 2588-95
3 ANDERSON, G.W.J., ARNOLD, C.P., WATSON, N.R., and ARRILLAGA, J.: 'A
new hybrid ac-dc transient stability program'. Int. conf. on Power systems transients
(IPST), September 1995, pp. 535-540
4 ANDERSON, G.W.J.: 'Hybrid simulation of AC-DC power systems'. PhD thesis,
University of Canterbury, New Zealand, 1995
5 EGUILUZ, L.I., and ARRILLAGA, J.: 'Comparison of power definitions in the
presence of waveform distortion', Int.J. ofElectr. Eng. Educ, April 1995, 32, (2),
pp. 141-153
6 FRYZE, S.: 'Wirk-, blind- und scheinleistung in elektrischen stromkreisen mit
nichtsinusformigem verlauf von strom und spannung'. Elektrotech. Z., June 1932,
pp. 596-599
7 HEFFERNAN, M.D., TURNER, K.S., ARRILLAGA, J., and ARNOLD, C.P.:
'Computation of AC-DC system disturbancesparts I, II and III', IEEE Trans.,
November 1981, PAS-100, (11), pp. 4341-63
8 CLARKE, E.: 'Circuit analysis of ac power systemsvol. II' (John Wiley & Sons
Inc., New York, 1950)
9 KIMBARK, E.W.: 'Power system stabilityvol. Ill, synchronous machines'
(Dover Publications Inc., 1968)
10 REEVE, J., and ADAPA, R.: 'A new approach to dynamic analysis of AC networks
incorporating detailed modelling of DC systemsparts I and II', IEEE Trans.,
October 1988, PD-3, (4), pp. 2005-19
11 WATSON, N.R.: 'Frequency-dependent ac system equivalents for harmonic
studies and transient converter simulation'. PhD thesis, University of Canterbury,
New Zealand, 1987
Appendix I
for = 1 ^ N
m= l - N
'
At each iteration of the N - R method, the nonlinear problem is approximated by the linear-matrix equation. The linearising approximation can
best be visualised in the case of a single-variable problem.
In Figure I.I xp is an approximation to the solution, with error Axp at
iteration p. Then
f(xp + Axp) = 0
(1.2)
+ .
(1.3)
If the initial estimate 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(xp) + Axpf'(xp) = 0
(1.4)
or
a.5>
The new value of the variable is then obtained from
xp+1 =xp + Axp
(1.6)
364
^solution
Figure LI
f(xp) = -JAx*
(1.7)
3
dx
(1.8)
Newton-Raphson method
365
Sparsity programming
In conventional matrix programming, double-subscript arrays are used for
the location of elements. With sparsity programming,1 only the nonzero
elements are stored, in one or more single vectors, plus integer vectors for
identification.
For the admittance matrix of order n the conventional storage requirements are n2 words, but by sparsity programming 66 + 3n words are
required, where b is the number of branches in the system. Typically,
b = 1.5n 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 factorisation
To solve the Jacobian-matrix eqn. 3.3.7, represented here as
[AS] = U][AE]
(1.9)
for increments in voltage, the direct method is to find the inverse of (J] and
solve for [AE], i.e.
[AE] = (J]"1[AS]
(1.10)
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
366
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 |N 3 , compared to N 3 to find the inverse. With sparsity programming the
number of operations varies as a factor of N. If rows are normalised, N
further operations are saved.
Optimal ordering
In power-system power 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 minimise the accumulation of nonzero
terms, and hence conserve sparsity, rather than minimising round-off error.
The diagonals are used as pivots.
Optimal ordering of row eliminations to conserve sparsity is a practical
impossibility due to the complexity of programming and time involved.
Instead, two types of semi-optimal scheme are commonly used, i.e. preordering2 and dynamic ordering.3
1.3 References
1 OGBUOBIRI, E.C., TINNEY, W.F., and WALKER, J.W.: 'Sparsity-directed decomposition for Gaussian elimination on matrices', IEEE Trans., 1970, PAS-89,
(1), pp. 141-50
2 STOTT, B., and HOBSON, E.: 'Solution of large power-system networks by
ordered elimination: a comparison of ordering schemes', Proc. IEE, 1971, 118, (1),
pp. 125-34
3 TINNEY, W.F., and WALKER, J.W.: 'Direct solutions of sparse network equations
by optimally ordered triangular factorization', Proc. IEEE, 1967, 55, (11), pp.
1801-1809
INDEX
Index Terms
Links
A
Abnormal modes
AC load flow
291
61
control equations
77
converter variables
72
mismatch equations
per unit system
73
80
75
146152
74
sequential method
162
unbalanced
146
154
170
73
78
293
21
62
152
287
293
104
Additional delay
290
Admittance matrices
Analytic Jacobian
121
Asymmetry
170
280
312
27
34
B
Back-to-back
269
Backward Euler method
206
Bergeron
224
40
Index Terms
Links
C
Cable
218
228
241
23
34
135
170
Chatter
204
220
CIGRE
223
benchmark model
372
109
167
202
248
263
312
336
340
348
13
47
75
89
309
12
22
29
41
130
134
275
303
309
15
46
102
19
90
47
310
15
117
130
162
222
288
294
307
371
Commutating voltage
Commutation
angle
290
failures
margin (see also Extinction angle)
307
reactance
11
310
13
233
Index Terms
Links
78
191
222
15
46
223
17
63
68
86
117
121
128
133
145
154
167
180
293
301
297
Convergence
Convergence factor
139
Convergence tolerance
74
133
180
Converter control
78
189
215
252
258
223
289
23
38
45
77
162
222
288
307
23
46
77
96
162
297
12
46
119
characteristics
constant current
constant power
firing instant
124
power modulation
298
Converter transformer
13
17
24
36
72
94
230
connection
leakage reactance
magnetising characteristic
on-load tap changer
17
155
134
230
236
233
252
77
134
235
252
139
139
saturation
36
Index Terms
Links
Converter variables
18
71
75
Convolution
51
112
227
36
252
Cross-modulation
35
Current mismatch
116
166
Curve fitting
228
241
335
387
39
106
252
298
331
D
Damping
206
DC faults
264
DC ripple
10
123
69
84
163
80
91
162
Delay angle
15
19
46
73
110
120
130
288
Direct axis
124
278
Dynamic stability
275
dq axes
278
E
Eigenvalue
Electromagnetic Transients Program
124
224
385
310
324
Chapter 6
Chapter 8
Chapter 6
EMTDC Program 4
Chapter 5
Index Terms
Links
Equivalent PI model
286
ESCR
367
Euler
51
Exciter
Extinction angle
206
385
265
47
78
125
177
258
324
61
79
68
88
F
FACTS
201
145
162
Fast Fourier transform
339
341
36
87
94
32
41
46
77
87
110
120
190
288
12
46
120
124
190
100
Firing angle
168
13
292
15
35
341
Index Terms
Links
Fourier series
23
43
135
Fourier transform
41
177
335
276
278
283
328
333
31
177
225
324
348
355
27
109
275
61
86
182
207
259
Frames of reference
Frequency dependent
Frequency domain simulation
G
Gaussian elimination
199
Gibbs phenomena
136
GTO
203
H
Harmonic
currents
3136
domain
41
flow
impedances
135140
145
37
40
254
355
interaction
162
169
phasors
41
45
54
sources
136
transfers
30
History term
195
206
217
220
228
Hydro turbine
282
Index Terms
Links
I
Inertia constant
277
Instability
36
189
252
Integration
176
207
301
181
184
191
292
295
304
345
180
207
234
303
185
203
11
16
78
A-stable
385
206
implicit
301
numerical
385
Runge-Kutta methods
383
-stable
303
step length
trapezoidal method
Interpolation
Inversion
194
98
deionisation angle, see extinction
angle
extinction angle
Iterative methods
47
78
41
47
62
84
89
109
117
153
180
229
294
296
301
344
Index Terms
Links
J
Jacobian
303
numerical
117
analytical
121
Jacobian matrix
119
6568
79
82
91
117
153
166
K
Kron reduction
165
L
Lead-lag circuit model
282
169
Load characteristics
285
293
39
252
M
Magnetising reactance
MAP
104
MATLAB
210
Matrix sparsity
287
Mismatch equations
4347
219
241
78
112
150
Mismatch functions
65
113
Modulation theory
26
252
84
Multiterminal
307
301
Index Terms
Links
N
NETOMAC
186
216
61
66
79
91
145
229
293
295
363
43
46
65
67
114
116
145
153
34
93
Nodal formulation
196
198
Newton-Raphson
equations
67
flow diagram
66
Jacobian elements
79
Newtons method
312
631
194
212
233
238
282
285
307
351
Numerical Jacobian
117
Numerical stability
133
210
241
O
Ordering of sparse matrices
210
287
P
Parks transformation
278
74
217
Index Terms
Phase locked oscillator
Links
44
189
222
277
factor
15
292
331
flow
115
326
modulation
298
Q
QESCR
369
Quadrature axis
124
284
15
61
6671
85
95
101
298
308
335
16
78
R
Reactive power
360
Rectification
11
characteristics
77
controlled
16
delay angle
15
19
46
36
119
136
253
335
385
3
176
268
213
Index Terms
Links
S
Saliency
90
283
Sampling functions
41
51
39
252
38
124
150
338
36
102106
259
186
221
35
94
146
367
Six pulse
9
259
Smoothing reactor
34
116
Sparsity
ordering
programming
Speed governor
287
368
61
365
280
312
133
210
Stability,
numerical
241
representation of plant in
network
301
381
263
265
Step length
181
184
199
217
292
296
323
343
Stub line
Subharmonic
199
34
36
Index Terms
Subtransient reactance
Links
90
Swing curves
278
Switching system
130
Switching terms
119
214
133
22
276
angular momentum
277
controllers
280
damper windings
337
inertia constant
277
models
214
338
338
Parks transformation
278
saliency
saturation
subtransient reactance
90
337
265
283
280
90
transient reactance
279
338
Synchronous reactance
279
System damping
106
214
253
T
TACS (transient analysis of control
systems)
215
259
177
298
Index Terms
Links
Three phase
faults
251
316
338
356
power flow
145
representation
141
152
Chapter 6
275
292
323
37
252
230
269
278
reactance
279
stability
Transient stability program structure
216
275
301
324
361
302
327
334
352
Transmission line model
196
223
241
Transmission lines
196
211
286
224
286
Trapezoidal integration
301
302
Travelling wave
226
Triangular factorisation
365
equivalent PI
282
thermal
282
2736
189
259
110
Index Terms
Links
U
UMEC (unified magnetic equivalent circuit)
Unbalance
230
3
17
19
34
43
111
124
157
160
332
338
203
extinction
184
191
switching
184
Valve group
189
220
Voltage mismatch
114
V
Valve