z-plane synthesis of response functions and

interpolators in z-transformelectromagnetic
transient analysis in power systems
Prof. W.D. Humpage, B.Sc, Ph.D., K.P. Wong, M.Sc, Ph.D., and T.T. Nguyen, B.E.
Indexing terms: Power systems and plant, Transmission and distribution plant, Mathematical techniques,
Matrix algebra, Power transmission, Electromagnetic transients, z-transforms
Abstract: In the original development of z-transform electromagnetic transient analysis in power systems,
transmission-line impulse responses and surge-impedance functions formed initially in the frequency domain
were mapped into the z-plane by bilinear transformation, although it is known that this transformation is
one that gives rise to a form of distortion error. In now providing means by which this error may be avoided,
methods are investigated in the paper by which transmission-line response functions can be synthesised
directly in the z-plane. Given the high accuracy in function definition which this can provide, an examination
is made of the relationships that are to be satisfied between the step length in the time-domain and trans-
mission-line aerial-and earth-mode transit times. In particular, it is shown that residual errors arising from the
frequency dependence of the earth-mode transit time are amenable to analysis and then, based on this, to
control by interpolation. Drawing on trial analyses of a 400 kV single-circuit untransposed line, the paper
then evaluates the further advances in z-transform electromagnetic transient analysis that these further
developments can achieve.
List of principal symbols
F(cS),F(z) = matrices of transmission-line forward-
impulse responses in the frequency domain
and the z-plane, respectively
Z
0
( CJ ) , Z
0
( Z) = surge-impedance matrices
X(co) = matrix of propagation coefficients
a(co), p(co) = real and imaginary parts of X(co), respect-
ively
v(to), V(Z) = vectors of transmission-line voltages
i'(co), i(z) = vectors of transmission-line currents
F
s
(z),F
r
(z) = forward characteristics at sending and
receiving ends, respectively
B
8
(z),B
r
(z) = backward characteristics at sending and
receiving ends, respectively
H(z) = general z-plane response function
H
s
(
z
i)
=
specified value of H(z) at the rth sampling
point
Q(6) = weighted error function to be minimised
9 = vector of function coefficients to be found
W(co
t
) = positive real weighting function
N = number of product terms in H(z)
M = number of frequency-domain and z-plane
sampling points
Pik,P2k = pole pairs of H(z);k= 1, 2, . . . ,N
r\k,
r
7k = zero pairs of H(z)\ k = 1, 2, . . . , N
T = wave-propagation transit time
At = step interval in time-domain recursive
sequences
/ = transmission-line length
m = integer value of ratio T/At
e = T/At m
co = angular frequency
z = z-transform operator
Subscripts s and r identify variables at the switching and
remote ends of the transmission line, respectively.
1 Introduction
While several different methods of electromagnetic transient
analysis in power systems [14] have now been developed
Paper 1868C (P9) received 15th June 1981
The authors are with the Department of Electrical & Electronic Engin-
eering, University of Western Australia, Crawley, Western Australia
6009
104
0143-7046/82/030104 + 07 $01.50/0
and widely applied, recent work [5,6] has investigated an
analysis formulation which is based on the z-transform. In
the new formulation, transmission-line equations in the
frequency domain are successively transformed first into the
z-plane and then into the time domain. Thereafter, the formu-
lation is wholly in the time domain, and the sequences in
solution, to which the steps of transformation through the
z-plane lead, are of recursive form. The earlier work has
sought to develop this analysis method and to evaluate its
principal properties in practical application.
Corresponding to the separate stages of transformation
from which the method derives, transmission-line forward-
impulse-response and surge-impedance functions, formed
initially in the frequency domain, were mapped to the z-plane
in the original work [5,6] by bilinear transformation. How-
ever, it is generally acknowledged that this transformation is
one which introduces a form of distortion error. High accu-
racy in response-function definition is achieved over an initial
range of frequency; beyond which progressively increasing
errors arise. When nominating the step length in the time
domain, it is necessary to bear this in mind, and to choose
the step length so that the distortion errors are negligibly
small over the range of frequencies relevant to the electro-
magnetic transient modes of system operation for which
analysis is to be carried out. The first purpose of this paper
is that of developing methods by which transmission-line
forward-impulse-response and surge-impedance functions may
be synthesised directly in the z-plane; thereby avoiding
altogether the errors arising when they are formed by bilinear
transformation from their counterparts in the frequency
domain.
It then remains to satisfy the requirement that the step
length in the time domain should be an integer submultiple
of wave transit times. The different wave transit times for
aerial and earth modes of propagation are usually accommo-
dated without difficulty in meeting this requirement; but the
frequency dependence of the earth-mode transit time rep-
resents a further source of errors. The second purpose of the
paper is that of showing that residual errors due to the non-
linear frequency dependence of earth-mode wave transit
times are amenable first to systematic analysis and then to
a substantial measure of control by interpolation.
In distinct, but closely related ways, direct z-plane response
synthesis and earth-mode interpolation offer further advances
IEE PROC, Vol. 129, Pt. C, No. 3, MA Y 1982
in the recently derived analysis fomulation [5, 6] . No previous
paper has referred to the developments in z-transform electro-
magnetic transient-analysis methods in power systems to
which this paper is devoted.
2 General form of z-plane functions
Transmission-line equations in the modal axes rearrange to
the scattering form [5, 6]
v
s
(to) Z
o
(oo)i
s
(cS) = F(oS) [v
r
(cS) + Z
o
(co)/
r
(co)] (1)
v
r
(co) Z
o
(oj)/
r
(co) = F(cS) [v
s
(CJ) + Z
o
(w)/
s
(CJ)] (2)
in which
F(cS) = exp{-A(co)/} (3)
In eqns. 1 and 2, Z
0
(co) is the matrix of surge-impedance
functions and F(co) is the matrix of forward impulse
responses. In eqn. 3, A(co) is the matrix propagation co-
efficients and / is the transmission-line length. Eqns. 1 and
2 transform directly into the z-plane to give:
v
s
(z)-Z
0
(z)i
s
(z) = F(z)[v
r
(z) + Z
0
(z)i
r
(z)]
v
r
(z)-Z
0
(z)i
r
(z) = F{z) [v
s
(z) + Z
0
(z)/
s
(z>]
where
F(z) =
Z
b
(z) = Zt{Z
0
(u>)}
(4)
(5)
(6)
(7)
The first purpose of this paper is to form F(z} and Z
0
(z)
directly as z-plane functions, rather than by deriving these
from rational-fraction forms of F(oS) and Zo(co) and the
z-transform, as eqns. 6 and 7 imply. As the steps involved
in achieving this are closely similar for F(z) and for Z
0
(z),
it is useful to adopt a general function H(z) in the detailed
working, and in terms of which both F(z) and Z
0
(z) can
subsequently be interpreted. On this basis, the factored
second-order form widely adopted in digital filter design is
used in which:
(8)
c
k
z~ +d
k
z
- 2
For this form, the'requirement is therefore one of finding
the a
h
, bk, Ck, d
k
and A coefficients of eqn. 8. In this work,
the quasi-Newton minimisation method, in conjunction with
a weighted mean-square error function, is used for this pur-
pose. The general form of eqn. 8 ensures that coefficient
values found by this process, and by which transmission-line
impulse response and surge-impedance functions are defined
in the z-plane, are always real.
then formed from:
M
Q{Q) =
f = i
(10)
in which p is positive and even. IV(co,) is a positive real
weighting function, and H(z
t
) is the value of the//(z) function
of eqn. 8 at the /th sampling point. 9 is the vector of unknown
coefficient values to be found so that:
i ' = a
Xt
bi
i
c
ll
d
lt
...
t
a
ki
b
h
,c
k
,d
h
,A
01)
In defining the H(z) function, it is required to find a set of
coefficient values 0* such that:
Q(6*) = minmax
M
- \H,(z
t
)\]\
P
(12)
In using the quasi-Newton minimisation algorithm for this
purpose, the case is taken for which p = 2, so that Q(d) is
of weighted mean-square error form. It is then convenient
to separate the^l coefficient in eqn. 8 by using:
H(z
t
) = AG(
Zi
)
so that:
T-2
(13)
(14)
Separating A correspondingly in the vector of coefficients
gives:
where
0' = ai,bi,ci,d
lt
. . . ,a
k
,b
k
,c
k
,d
k
05)
(16)
The weighted mean-square form of the error-function then
becomes:
Q(d>,A)=
M
07)
The differentiations dQ(<p,A)/dA and dQ((f>,A)l^(l> now
follow, and unknown coefficient values in the H(z) function
are found from the following iterative sequence
(a) At iteration step n + 1, calculate G(zf) using:
c(z
,
(18)
3 Quasi-Newton minimisation
3.1 Error function and coefficient-finding algorithm
At any frequency value co
f
, specified values for the H(z)
function of eqn. 8 are available from the known forward-
impulse-response and surge-impedance functions in the fre-
quency domain together with the z-transform operator:
= exp{/cJ,-Af}
(9)
In eqn. 9, At is the sampling interval, or step-interval, in the
time domain. Specified values formed in this way, at the
ith sampling point in the z-plane, are denoted by H
s
(z
t
). A
frequency-weigh ted error function [7] to be minimised is
in which a, b%, c and d are coefficient values from the
immediately previous step. The iterative sequence is started
when n = 0 from initially estimated values of the coefficients.
(b) From the known transmission-line response functions
in the frequency domain, calculateH
s
(z
t
) for / = 1, 2, . . .,M.
(c) Calculate the optimum value for the A coefficient
from:
M
in+l
\A\
M+l _
M
(19)
W(wi)IG(zi)l
n+l
IEEPROC, Vol. 129, Pt. C, No. 3, MA Y1982 105
(d) Calculate the differentials:
*\G(z
t
)\ *\G(z
t
)\ d\G(z
t
)\
da
n
h
for
k = 1 , 2 , . . . , 7 V
i = 1 , 2 , . . . , M
From these individual differentials form:
*\G(z,)\
for
A: = 1, 2, ...,N
(e) Calculate the gradient from:
A
n + 1
}
M
= 2
1 = 1
(20)
" I " . (*,)!]
30*
(/ ) Calculate coefficient values from [8]
(21)
in which //(0, ;4
n + 1
) approximates to the inverse of the
Hessian matrix, and is updated as in Reference 8. In eqn.
21, 7 is the value of i// which minimises Q(<Pk + \ps,A
n+1
)
where s is the direction of search as formed from:
s = -
(22)
(g) Repeat the steps (a)-(f) until a convergence criterion
is satisfied.
3.2 Guaran teeing stable system functions
Coefficient values formed from the procedure of Section
3.1 may not necessarily lead to F(z) and Z
0
(z) functions
which are themselves inherently stable. For stability, the
pole and zero locations of the functions should lie within
the unit circle in the .7-plane. .From these conditions, checks
and controls may be derived by which "tife iterative procedure
of Section 3.1 can be constrained so that the coefficient
values found from it always lead to stable system functions.
From the form of H(z) in eqn. 8, there are pairs of poles
and of zeros corresponding to each of the terms in the product
series. If p
l f e
and p
7
k are pole locations deriving from the
fcth term, then:
Pik
and
k = 1, 2, . . . , 7V
k = 1, 2, . . . , / V
(23)
(24)
There are similar relationships between zero locations and
the ah and b
h
coefficients.
At each iteration in the algorithm of Section 3.1, all pole
and zero locations are evaluated from coefficient values
calculated in step (/ ). When these satisfy the criteria for
stable-system functions, the sequence proceeds 1o the next
iteration without further intervention. If they are not satisfied,
constraints are introduced before continuing to the next
iteration. If, foi example, the pole denoted by p
ik
in eqns.
23 and 24 is found to lie outside of the unit circle, then the
coefficient values determining it are recalculated using:
+ P
Pik
d
k
= Pik/Pik
(25)
(26)
Corresponding recalculations are made where a zero location
is found to violate the conditions for stability. With the
recalculated coefficients, the iterative process is continued
until convergence is achieved. At convergence, all pole and
zero locations fall within the z-plane unit circle; and the final
set of coefficient values are those for which the F(z) and
Z
0
(z) functions defined in terms of them are inherently
stable-system functions. On taking inverse z-transforms, the
recursive sequences for solution in the time domain are
inherently stable sequences.
4 z-transform electromagnetic transient analysis with
direct synthesis of z-plane functions
In Fig. 1 is shown the earth-mode forward impulse response
function for a 160 km, 400 kV transmission circuit, for which
the principal data is that of the earlier work [5]. In the
original z-transform development [5], this is represented
in the s-plane by:
(1 +/ i s+/
2
s
2
) exp { smAt}
(l+f
3
s+f,s
2
)(l+
gl
s+g
2
s
2
)
;- for s =
(27)
If T is the wave-propagation transit time? then m in eqn.
27 is given by:
m = T/At
(28)
Coefficient values/! /
4
and g
x
and g
2
are formed by least-
squares minimisation, and the F(s) function is transformed
to the z-plane using:
z - \
(29)
In the z-plane:
F{Z) =
b
2
z~
2
+ b
3
z
- 3
b
A
z
- 4
(30)
l!
c c
! &
c <u
1.0
0.8
C.6
C.t,
0 2
0
2 3 A 5
frequency , rad/s x 10
Fig. t Distortion error in earth-mode forward impulse response
160 km, 400 kV single-circuit transmission line
from transmission-line data
calculated from 4th-order z-plane function derived by bilinear
transformation with step length SO jus
106
IEEPROC, Vol. 129, Pt. C, No. 3, MAY 1982
in which the a and b coefficients are found directly from the
f&n&g coefficients of eqn. 27.
Interpreted in the co plane using z = exp {jcoAt}, the
F(z) function of eqn. 30 is also shown in Fig. 1. The differ-
ence between this and the actual function calculated from
400 kV transmission-line data [5] is a measure of the distor-
tion error to which the bilinear transformation of eqn. 29
gives rise. Similar, though generally smaller, distortion errors
arise in the surge-impedance function when formed in the
z-plane by bilinear transformation from the frequency domain.
1 2 3 A 5 6
frequency, rad/s xlO
A
Fig. 2 Accuracy of forward impulse response representation from
direct z-plane synthesis
Earth mode of 160 km, 400 kV single-circuit transmission line
from transmission-line data
calculated from frequency-weighted 2nd-order function syn-
thesised directly in the z-plane
Turning now to where z-plane functions are formed directly
using the development of Section 3, Fig. 2 indicates that the
distortion error of Fig. 1 is now avoided. Moreover, the high
accuracy with which the transmission-line impulse response
is represented in Fig. 2 is achieved with a second-order form
for F(z), whereas fourth-order forms were used throughout
the earlier work. Taken further, first-order forms can also
be derived. From numerous trial studies, it is the authors'
preference to use second-order forms both for F(z) and for
Z
0
(z) in their application to ground modes, and to use first-
order forms for both in the case of aerial modes. A typical
accuracy of representation for the earth-mode surge
impedance is shown in Fig. 3.
Distortion errors in forming F(z) and Z
0
(z) by transform-
ation from corresponding frequency-domain functions by no
means transcribe proportionally into electromagnetic transient
solutions. As the earlier work confirmed, exceptionally high
solution accuracy Can" be achieved when fourtlvorder represen-
tations of response functions are used in conjunction with
>tep-length settings of about one-tenth of wave transit times.
Avoiding distortion errors allows a similar accuracy to be
ichieved with longer step lengths, thereby further reducing
computing time requirements. As the step length is progress-
ively increased, however, a time-shift form of error appears
which is unrelated to the accuracy by which the F(z) and
0 1 2 3 4
frequency . rad/s xiO*
1
Fig. 3 Surge-impedance function
Earth mode of 160 km, 400 kV single-circuit transmission line
from trnasmission-line data
calculated from frequency-weighted 2nd-order function syn-
thesised directly in the z-plane
Z
0
(z) functions are defined. This error is shown in Fig. 4
when the step length is increased to 110 jus. Here, the z-
transform solution is compared with a frequency-domain
solution for which the frequency separation is 200 rad/s,
the upper frequency limit is 10
s
rad/s, and the step length in
the numerical evaluation of inverse Fourier transforms is
20jus. For the single line solution of Fig. 4, the z-transform
computing time is 120 ms. This is less than 0.1% of the com-
puting time required for the frequency-domain solution with
which it is compared. The source of the time-shift error in
Fig. 4 is in the relationship between the step length and the
wave transit times. It remains to examine this source of
700
580
',60
340
o 220
100
-20
-140
-260
-380
-500
3
C >
O - *
n o
3 -
z o
a >
ti me. ms
Fig. 4 Switching overvoltage solution with extended step length
160 km, 400 kV single-circuit transmission line: phase a energisation
frequency domain
z-transform with 110 /is step length
errors and then to counter it. The origin and nature of the
errors suggest that some form of interpolation might be
considered in reducing them.
5 Interpolation
5.1 Step-interval and wave transit-time relationships
On separating out the propagation coefficient matrix A (co)
into its real and imaginary parts:
A(w) = a(co)+/p(w) (31)
Wave transit times for the individual modes are then given by:
T
t
= ftM//w i = 0, 1, 2, modes (32)
and where:
= diag(ft) i = 0, 1,2, modes (33)
A basic aspect of the derivation of z-transform electromagnetic
transient analysis is that of the step-length At being an integer
sub-multiple of the wave-transit time. Ideally, therefore:
m
{
At = T
t
i = 0, 1, 2, modes (34)
where m is a positive integer.
It is often possible satisfactorily to fulfil the requirements
of eqn. 34 by choice of At. For example, if the transit time
common to the two aerial modes for a particular single-circuit
transmission line is 450/zs, and that for the earth-mode is
500jus, then step-length settings of 50/JLS, 25/US, or 12.5 jus
allow eqn. 34 to be satisfied. In ^-transform electromagnetic
transient analysis, a setting of 50jus may safely be used.
However, this leaves the question of the frequency dependence
:>f the earth-mode transit time of which that for the trans-
mission line to which the response functions of Figs. 1-3
relate is shown in Fig. 5. This frequency dependence leads
:o a residual error in solution, the magnitude of which is
related to the contribution that the earth mode makes to the
complete solution. Always, however, where step-interval
values are less than about one-tenth of wave transit times,
the error is a small one.
IEEPROC, Vol. 129, Pt. C, No. 3, MAY 1982 107
Increasing the step length as a means of reducing computing
time is clearly limited by the requirement of satisfying the
integer relationships of eqn. 34. Where the choice of an
extended step length precludes eqn. 34 being satisfied strictly
for all modes, we use:
(ra
f
+ e,) Af = T
t
i = 0,1,2, modes
for
0 < e
t
< 1
(35)
(36)
700 r,
650
600
550,
1 3 U
frequency,
Fig. 5 Variation of earth-mode transit time with frequency
160 km, 400 kV single-circuit transmission line
5.2 The impulse response
In the frequency domain, the impulse-response F(co) is given
by:
F( CJ) = (37)
Drawing on eqns. 32 and 35
for a transmission line of length /.
It has previously been shown [5] that the imaginary part
of X(co) can be decomposed into components which have
separately linear and nonlinear dependences on frequency.
On this basis:
(38)
(39)
(40)
(41)
(42)
(43)
The impulse response matrix F(cS) is then given by:
where
F
m
(co) = exp [-{a(oj)+/p
n
(co)}/] exp {-
and
= exp { jetjAt)
for
e = diag {e,} / = 0, 1, 2, modes
In these terms, F
m
(CJ) is the forward impulse-response matrix
when the wave transit times and the step interval have integer
ratios, and eqn 34 is satisfied. When this is not the case, the
forward impulse response has the product form of eqn. 40.
5.3 Forward an d back ward charac ten's tics
From the scattering form of transmission-line equations in
the z-plane, forward characteristics at the sending and receiv-
ing ends of the line F
s
(z) and F
r
(z) are defined by:
F
r
(z) = v
r
(z) + Z
0
(z)i
r
(z) (45)
The backward characteristicsB
s
(z) and B
r
(z) are given by:
B
s
(z) = v
s
(z)-Z
0
(z)i
s
(z) (46)
B
r
(z) = v
r
(z)-Z
0
(z)i
r
(z) (47)
Forward and backward characteristics are then interrelated
by:
B
s
(z) = F(z)F
r
(z) (48)
B
r
(z) = F{z)F
s
{z) (49)
On using the product form for F(z), as in eqn. 40:
B,(z) = F
m
(z)F
e
(z)F
r
(z) (50)
B
r
(z) = F
m
(z)F
e
(z)F
s
(z) (51)
When F
e
(z) = 1, and F
m
(z) is known from the developments
of Section 3, eqns. 50 and 51 transform directly into the
time domain by taking inverse z-transforms. In achieving
transformation when F
e
(z)^ 1, it is useful to take the scalar
form of F
e
(z) in which case:
F
e
(z) = z (52)
More generally, eqn. 52 applies for each mode. We now
approximate z~
e
by the polynominal form:

1
+. . . +. 4
n
( e ) z - (53)
To assist in achieving a clear presentation of the steps in
analytical development now required, a lst-order truncation
of eqn. 53 is used, in which case:
B.(z) = F
m
(z)[A
0
(e) + A
l
(e)z-
1
]F
r
(z)
B
r
(z) = F
m
(z)[A
0
(e)
(54)
(55)
Following the observations of Section 4, we now use F
m
(z)
in eqns. 54 and 55 in the form:
+az~
l
+bz
+cz~
l
+dz ~
2
(56)
The derivation of coefficients, using the developments of
Section 4, is, therefore, that for a second-order function where
the A coefficient of eqn. 8 in Section 2 is 1.0. Substituting
for F
m
(z) from eqn. 56 into eqn. 54 and crossmultiplying:
[1
1
+dz~
2
]B
s
(z) =
(57)
Expanding the RHS of eqn. 57 and rearranging:
B
s
(z) = A
0
(e)z-
m
F
r
(z)+{aA
0
(
(58)
108
IEEPROC, Vol. 129, Pt. C, No. 3, MAY 1982
On taking inverse z-transforms of both sides of eqn. 58:
*. ( ' ) =A
0
(e)F
r
(t
n
-mAt)
+ {bA
0
(e) + A,(e)}F
r
{t
n
-(rn + 2)At}
+ bA
l
(e)F
r
{t
n
-(m + 3)At)
-cB
s
(t
n
-At)-dB
s
(t
n
-7At)
Similarly, using eqn. 55:
B
r
(t
n
) = A
0
(e)F
s
(t
n
-mAt)
+ {aA
0
(e) + A
1
(e)}F
s
{t
n
-(m + l)At}
+ {bA
0
(e) + A, (e)}F
s
{t
n
- (m + 2)At]
+ bA
1
(e)F
s
{t
n
-(m
(59)
(60)
- cB
r
(t
n
- At) - dB
r
(t
n
- 2At)
These sequences are identical in form with those of the earlier
work [5]; only coefficient values are modified by the presence
of A
0
(e) and A^e). Given these coefficient values, the sol-
ution proceeds precisely as in the original development [5].
But now a linear interpolator is embedded in the formulation.
By this, the solution errors arising when the time-domain
step length is not a precise submultiple of modal wave transit
times are controlled. Although the steps of including a process
of interpolation in solution have been developed for the
particular case when the interpolator is linear, the develop-
ment follows along the same lines when an interpolator of
higher order is to be included In particular, it appears from
numerous experiments that a quadratic interpolator is in
many cases the most satisfactory choice in the present appli-
cation.
5.4 Lagrange coefficients
To complete the interpolation scheme, it remains only to find
values for the coefficients A
0
(e) and A
x
{e) in eqns. 59 and
60. More generally, it is required to find values for the coef-
ficients of the polynomial by which z~
e
is approximated.
Here we use the Lagrange coefficients [9] for this purpose. If
there are q terms in the polynomial and q is even, then:
(61)
k = 0, 1 , . ..,(q- 1)
When q is odd:
{(q-\)/2 + k}\{(q-\)l2-k}\(e-k).
k = 0, 1, . . . , ( </ - ! ) (62)
Knowing coefficient values from these forms, the recursive
sequences for solution, as in eqns. 59 and 60 for the particular
case of a linear interpolator, are defined. For clarity, we have
given here the developments for the solution of the forward
and backward characteristics. In turn, these are solved from
the relationships of eqns. 44-47 involving transmission-line
current and voltage variables and the suige-impedance func-
tion. The detailed steps here are those of the original work
[5].
IEEPROC, Vol. 129, Pt. C, No. 3, MA Y 1982
6 Solutions with interpolation
While the interpolation scheme of Section 5 leading to modi-
fied coefficients in the time-domain recursive sequences may
be applied in principle to all modes, trial applications indicate
that it has its greatest effectiveness when applied to earth
modes. The discontinuities in aerial-mode-solution values
can in some cases reduce the error control which interpolation
can provide. By comparison, the more rounded and damped
responses for earth modes allows interpolation safely to be
applied in all applications. Moreover, as the aerial-mode
transit times are almost constant, the step interval can be
selected to satisfy the integer relationships of eqn. 34 for
these modes. Thereafter, the earth-mode interpolator can
counter very effectively the time-shift form of errors when
the step length is not an integer submultiple of the transit
time in this mode. Included in this deviation from an integer
relationship is the dependence of the earth-mode transit
time on frequency.
Fig. 6 Solution with interpolation
Switching transient for 160 km, 400 kV single-circuit transmission
line
frequency domain
z-transform with earth-mode interpolator and 110 MS step
length
In Fig. 6, solutions are shown for the same conditions
and step length as those of the solution of Fig. 4; but now
with a quadratic interpolator embedded in the earth-mode
sequences. The effectiveness of the interpolator in countering
the solution errors of Fig. 4 is clear from Fig. 6. The step
length is 110/is, leading to a computing-time requirement of
120ms for a 10 ms solution period. The accuracy in solution
is characteristic of z-plane electromagnetic transient analysis
and is achieved for all variables of which only one is shown in
Fig. 6.
7 Conclusions
The paper reports further developments in z-transform electro-
magnetic transient-analysis methods in power systems. A
direct derivation of transmission-line forward impulse
responses and surge-impedance functions in the z-plane
coincides with that of most function synthesis in digital
filter analysis and design, and represents the preferred ap-
proach. In the present developments it has led to z-plane
functions of lower order than those of the initial work, and to
longer step-length settings. Both measures have a direct bearing
on computing-time requirements in analysis. By them, further
and substantial reductions in computing time have been
achieved, while at the same time preserving the high accuracy
and stability inherent in z-transform analysis. The detailed
formulation, the recursive solution sequences that follow
from it and the principal properties of the analysis method
are those of the original work [5]. Direct z-plane synthesis
procedures relate specifically to the derivation of coefficients
which first define transmission-line functions in the z-plane,
and then the time-domain recursive sequences for solution to
109
which these lead. As such, the procedures are those of the
preparatory phase of analysis which follows transmission-line
parameter calculations and precedes electromagnetic transient
analysis itself.
While progressively increasing step-interval settings in the
time domain is not usually consistent with the requirement
for values to be integer submultiples of modal wave transit
times, interpolation can counter the time-shift form of error
which appears in solution when integer ratios are not achieved.
The paper has shown how an interpolator may be embedded
in the z-transform analysis formulation for this purpose.
It modifies certain coefficients in the recursive sequences
in solution, but it does not change their essential structure.
It also avoids superimposing an additional stage of numerical
processing on the basic z-transform algorithm, which proposals
for using interpolation might initially suggest. Once developed
for a single transmission-line model, as in this paper, interp-
olation is of considerable benefit in multinode analysis [10]
where a very wide range of transit times can be encountered.
8 Acknowledgments
The authors are grateful to the Australian Research Grants
Committee for financial support, and to the West Australian
Regional Computing Centre for running their programs. The
generous support of power systems research in the Department
of Electrical and Electronic Engineering at the University
of Western Australian by Prof. A.R. Billings, together with
his professional co-operation at all times, is gratefully
acknowledged. T.T. Nguyen gratefully acknowledges the
award of a University Research Scholarship.
9 References
1 BICKFORD, J.P., and DOEPEL, P.S.: 'Calculation of switching
transients with particular reference to line energisation', Proc. IEE,
1967,114, (4), pp. 465-477
2 BATTISSON, M.J., DAY, S.J., MULLINEUX, N., PARTON, K.C.,
and REED, J.R.: 'Calculation of switching phenomena in power
systems', ibid., 1967,114, (4), pp. 478-486
3 DOMMEL, H.W.: 'Digital computer solution of electromagnetic
transients in single- and multi-phase networks', IEEE Trans., 1969,
PAS-88, pp. 388-396
4 SEMLYEN, A.,and DABULEANU, A.: 'Fast and accurate switching-
transient calculations on transmission lines with ground return
using recursive convolutions', ibid., 1975, PAS-94, pp. 561-569
5 HUMPAGE, W.D., WONG, K.P., NGUYEN, T.T., and SUTANTO,
D.: 'Z-transform electromagnetic transient analysis in power sys-
tems', IEE Proc. C, Gen., Trans. & Distrib., 1980, 127, (6), pp.
370-378
6 HUMPAGE, W.D., WONG, K.P., and NGUYEN, T.T.: 'Time convol-
ution and z-transform methods of electromagnetic transient analysis
in power systems', ibid., 1980, 127, (6), pp. 386-394
7 DECZKY, A.G.: 'Synthesis of recursive digital filters using the
minimum p-error criterion', IEEE Trans., 1972, AU-20, pp. 257-
263
8 FLETCHER, R.: 'A new approach to the variable metric algorithm',
Computer J., 1970, 13, pp. 317-322
9 SC1IAFER, R.W., and RABINER, L.R.: 'A digital signal processing
approach to interpolation', Proc. IEEE, 1973,61, pp. 692-702
10 HUMPAGE, W.D., WONG, K.P., and NGUYEN, T.T.: 'Development
of z-transform electromagnetic transient analysis methods for
multinode power networks', IEE Proc. C, Gen., Trans. & Distrib.,
1980,127, (6), pp. 379-385
Abstracts of papers published in other Parts of the IEE PROCEEDINGS
The following papers of interest to readers of IEE Proceedings Part C: Generation, Transmission & Distribution have appeared in
other Parts of the IEE Proceedings:
Eddy-current distribution in cylindrical structures caused
by rotating magnetic fields
P.P. YANNOPOULOS-LASCARATOS and PROF. J.A.
TEGOPOULOS
IEE Proc. B, Electr. Power Appl., 1982, 129, (2), pp. 64-74
Analytical solutions are given for the eddy currents (magni-
tude and phase) induced in a cylindrical shell of infinite
length and. thickness by a. 3-phase winding of finite length
placed coaxially inside the shell. The flux densities and
issociated loss are also determined. Solutions are obtained
in integral form by solving Maxwell's equations and applying
appropriate boundary conditions. Use is made of the magnetic
/ector potential. The loss is calculated by applying Poynting's
theorem on the internal surface of the shell. Numerical results
are obtained on a digital computer for different radii of the
shell by using the Gauss-Laguerre method. The limitations
of the method are discussed. The conducting material is
nonmagnetic and the solutions are valid for all frequencies
smaller than those of the optical spectrum. Tests were made
on an aluminium shell and measurements of field values in
the air were compared with the corresponding theoretical
results.
Dynamic response of a power-transformer winding under
axial short-circuit conditions
D.A. QU1NNEY
IEE Proc. B, Electr. Power Appl., 1982, 129, (2), p. 87
In a previous paper, a mathematical model was developed to
consider the dynamic response of a power-transformer winding
after it had been subjected to an axial short-circuit. In this
paper, an analytic solution, which is valid for small times t,
is produced which justifies the basic assumptions made in the
original paper.
110
IEE PROC, Vol. 129, Pt. C, No. 3, MAY 1982