388 IEEE TRANSACTIONS ON POWER APPARATUS AND SYSTEMS, VOL. PAS-88, NO.

4, APRIL 1969

annual availability function after a single, new machine is added, The data requirements are then primarily for the annual peak load
assuming that nothing else is altered. Let forecasts. It might be observed that the load model of the paper does
not require that the occurrence N of a load level of Li MW be an
C capacity of new machine, MW integer value. This might be useful in specifying the peak load varia-
A (M) annual availability of the basic system as a function of tion curve from historical data since, for instance, a model for a
the margin M 30-day period might include 35 load levels. This would permit using
M reserve margin more data to define the highest load levels in the peak variation
a availability of the new machine. curve. The load model also requires the specification of the mean
Thus, the new availability function becomes duration e of the peak periods. This value is a matter of judgment,
to be arrived at after a suitable study of load cycle data. Values of
A'(M) = aA(M - C) + (1 - a)A(M). (33) from 1/4 to over 1/2 of a day would seem to be appropriate for
various different systems.
This relationship may be used to consider the effect of adding dif- Concerning Mr. Watchorn's question about representing hydro
ferent machine capacities on the annual availability (i.e., the loss-of- units and plants (energy limits and head effects), we have not
load probability). A similar relationship may be worked out for the included these provisions in the analysis. However, we would like
frequency calculations. These are approximate since they assume no to note that the frequency and duration method will handle the same
change in maintenance outages for the old machines and that the problems that may be treated by loss-of-load probability.
new machine is not maintained during the year. We agree with Mr. Watchorn that the "economic criterion" would
The technique for constructing the load model suggested by Mr. be "by far the most meaningful aid to judgment." Even separate
Adler is similar to that used by the authors and their associates. economic criteria for bulk power supply and for distribution systems
That is, historical load data are examined statistically to establish would be of much value for system planning.
periodic (i.e., monthly or seasonal), per unit peak load variation Again we wish to thank the various discussers for their contribu-
curves, and monthly or seasonal peaks in per unit of the annual peak. tions. It is gratifying to observe the continuing interest in this area.

Digital Computer Solution of Electromagnetic

Transients in Single- and Multiphase
Networks
HERMANN W. DOMMEL, MEMBER, IEEE

Abstract-Electromagnetic transients in arbitrary single- or Among the useful features of this program are the inielusion of
multiphase networks are solved by a nodal admittance matrix nonlinearities, any number of switchings during the transient in
method. The formulation is based on the method of characteristics accordance with specified switching criteria, start from any
for distributed parameters and the trapezoidal rule of integration nonzero initial condition, and great flexibility in specifying
for lumped parameters. Optimally ordered triangular factorization voltage and current excitations of various waveforms.
with sparsity techniques is used in the solution. Examples and
programming details illustrate the practicality of the method. The digital computer cannot give a continuous history of the
transient phenomena, but rather a sequence of snapshot pictures
I. INTRODUCTION at discrete intervals At. Such discretization causes trulncation
errors which can lead to numerical instabilitY [3]. For this
THIS PAPER describes a general solution method for finding reason the trapezoidal rule was chosen for integrating the ordi-
the time responses of electromagnetic transients in arbitrary nary differential equations of lumped inductances and capaci-
single- or multiphase networks with lumped and distributed tances; it is simple, numerically stable, and accurate enough for
parameters. A computer program based on this method has practical purposes.
been used at the Bonneville Power Administration (BPA) and Branches with distributed parameters are assumed to be
the Munich Institute of Technology, Germany, for analyzing lossless; they will be called lossless lines hereafter. By neglecting
transients in power systems and electronic circuits [1], [2]. the losses (which can be approximated very accurately in other
ways, as will be shown) an exact solution can be obtained with
the method of characteristics. This method has primarily been
Paper 68 TP 657-PWR, recommended and approved by the used in Europe, where it is known as Bergeron's method; it was
Power System Engineering Committee of the IEEE Power Group
for presentation at the IEEE Summer Power Meeting, Chicago, Ill., first applied to hydraulic problems in 1928 and later to electrical
June 23-28, 1968. Manuscript submitted February 12, 1968; made problems (for historic notes see [5]). It is well suited for digital
available for printing April 10, 1968. The early stages of this work
were sponsored by the German Research Association (Deutsche computers [6]-[8]. In contrast to the alternative lattice method
Forschungsgemeinschaft) while the author was with the Munich for traveling wave phenomena [9] it offers important advantages;
Institute of Technology. for example, n-o reflection coefficients are necessary when this
The author is with Bonneville Power Administration, Port-
land, Ore. method is used.

Authorized licensed use limited to: ULAKBIM UASL - Yuzuncu Yil Universitesi. Downloaded on March 19,2021 at 08:04:28 UTC from IEEE Xplore. Restrictions apply.
DOMMEL: COMPUTER SOLUTION OF ELECTROMAGNETIC TRANSIENTS 389

The method of characteristics and the trapezoidal rule can TERMINAL k TERMINAL m
easily be combined into a generalized algorithm capable of
solving transients in any network with distributed as well as (a)
lumped parameters. Numerically this leads to the solution of a ik m(t) im k(t)
system of linear (nodal) equations in each time step. It will be
shown that lossless lines contribute only to the diagonal elements ek(j t k tr-)r)9 jemj(t)
of the associated matrix; off-diagonal elements result only from
lumped parameters. Thus a very fast and simple algorithm can (b)
be written when lumped parameters are excluded. However, no
such restrictions are imposed. Instead, the recent impressive Fig. 1. (a) Lossless line. (b) Equivalent impedance network.
advances in solving linear equations by sparsity techniques and
optimally ordered elimination [10] have been incorporated into Note that in (4) the expression (e + Zi) is constant when
an algorithm which automatically encompasses the fast solution (x - vt) is constant and in (5) (e - Zi) is constant when
of the restricted case and yet retains full generality. (x + vt) is constant. The expressions (x - vt) = constant and
(x + vt) = constant are called the characteristics of the differ-
II. SOLUTION FOR SINGLE-PHASE NETWORKS ential equations.
A digital computer solution for transients is necessarily a The significance of (4) may be visualized in the following
step-by-step procedure that proceeds along the time axis with a way: let a fictitious observer travel along the line in a forward
variable or fixed step width At. The latter is assumed here. direction at velocity v. Then (x - vt) and consequently (e + Zi)
Starting from initial conditions at t = 0, the state of the system along the line will be constant for him. If the travel time to
is -found at t = At, 2At, 3At, ... , until the maximum time tmn get from one end of the line to the other is
for the particular case has been reached. While solving for the d/v = dVL'C'
T= (6)
state at t, the previous states at t - At, t - 2At, . .. , are known.
A limited portion of this "past history" is needed in the method (d is the length of line), then the expression (e + Zi) encountered
of characteristics, which is used for lines, and in the trapezoidal by the observer when he leaves node m at time t - r must still
rule of integration, which is used for lumped parameters. In the be the same when he arrives at node k at time t, that is
first case it must date back over a time span equal to the travel em (t r) + Zim,k (t - r) = ek(t) + ZQ(-ik,m (t) )
-
time of the line; in the latter case, only to the previous step.
With a record of this past history, the equations of both methods (currents as in Fig. 1). From this equation follows the simple
can be represented by simple equivalent impedance networks. two-port equation for ik,m
A nodal formulation of the problem is then derived from these
networks. ik,m(t) = (I/Z)ek(t) + Ik(t - )
and analogous (7a)
Lossless Line
?im,k(t) = (1/Z)em(t) + Im(t - T)
Although the method of characteristics is applicable to lossy
lines, the ordinary differential equations which it produces are with equivalent current sources Ik and i,m7 which are known at
not directly integrable [8]. Therefore, losses will be neglected at state t from the past history at time t -r,
this stage. Consider a lossless line with inductance L' and
capacitance C' per unit length. Then at a point x along the line Ik (t -r) = - (1/Z)em (t - r) -i?,k (t - r)
voltage and current are related by (7b)
-de/8x= L'I(i/t) (la) Im(t - r) = - ik,m(t - r).
(1/Z)ek(t - r) -

- ai/d = C' (9e/ot). (lb) Fig. 1 shows the corresponding equivalent impedance network,
which fully describes the lossless line at its terminals. Topo-
The general solution, first given by d'Alembert, is logically the terminals are not connected; the conditions at the
i (x, t) = fl (x- vt) + f2 (X + Vt) (2a) other end are only seen indirectly and with a time delay r
through the equivalent current sources I.
e (x, t) = Z f (x - Vt) - Z-f2 (X + Vt) (2b)
with fi (x - vt) and f2 (x + vt) being arbitrary functions of the Inductance
variables (x - vt) and (x + vt). The physical interpretation of For the inductance L of a branch kc, m (Fig. 2) we have
fi (x - vt) is a wave traveling at velocity v in a forward direction
and of f2(x + vt) a wave traveling in a backward direction. ek- em = L (dik,m/dt) (8a)
Z in (2) is the surge impedance, v is the phase velocity
which must be integrated from the known state at t - At to
Z= V/L'/C' (3a) the unknown state at t:
v = 1/1VL'C'. (3b) ik,m(t) = ik,m(t - At) +
(ek - em) dt. (8b)
Multiplying (2a) by Z and adding it to or subtracting it from L AtZ
(2b) gives Using the trapezoidal rule of integration yields the branich
e(x, t) + Z i(x, t) = 2Z fi(x - vt) (4) equation
e(x, t) Z.i(x, t) = -2Z.f2(x + vt). (5) ik,m (t) = (At/2L) (eAk (t) - em (t)) + Ik,m (t - At) (9a)

Authorized licensed use limited to: ULAKBIM UASL - Yuzuncu Yil Universitesi. Downloaded on March 19,2021 at 08:04:28 UTC from IEEE Xplore. Restrictions apply.
390 IEEE TRANSACTIONS ON POWER APPARATUS AND SYSTEMS, APRIL 1969
NODE k - - NODE m
NODEk 4....° NODEm
(a)
Ik t At)
(a)
Ik m(t-At

ek(t R em(t
Ikot) R- jm DATUM DATUM
DATUM DATUM
(b) (b)

Fig. 2. (a) Inductance. (b) Equivalent impedance network. Fig. 3. (a) Capacitance. (b) Equivalent impedance network.

(2) (1)
INITIALLY:
ik m(t) l-V4V \ Y ] TRIANGULAR FACTORIZATION
l r0 - - -
+
Ik't R emRt) \ HIN EACH TIME STEP:
[j
(2)(2) FORWARD SOLUTION
BACK SUBSTITUTION
DATUM DATUM

Fig. 4. Resistance. Fig. 5. Repeat solutions of lillear equations.

where the equivalent current source Ik,m is again known from known [3] and will not be explained here. The result is a system
the past history: of linear algebraic equations that describes the state of the system
at time t:
Ik,,m(t- At) E Yle (t) I Ei (t) I EII (12)
= ik,m(t - At) + (At/2L) (ek(t - At) - em(t - At)). (9b) with
The discretization with the trapezoidal rule produces a trun- [Y] nodal conductance matrix
cation error of order (At)3; if At is sufficiently small and cut in [e (t)] column vector of node voltages at time t
half, then the error can be expected to decrease by the factor [i (t)] column vector of injected node currents at time t
1/8. Note that the trapezoidal rule for integrating (8b) is (specified current sources from datum to node)
identical with replacing the differential quotient in (8a) by a [I] known column vector, which is made up of known
central difference quotient at midpoint between (t - At) and t equivalent current sources I.
with linear interpolation assumed for e. The equivalent impe-
dance network corresponding to (9) is shown in Fig. 2. Note that the real symmetric conductance matrix [Y] remains
unchanged as long as At remains unchanged. It is, therefore,
preferable, though not mandatory, to work with fixed step width
Capacitance At. The formation of [Y] follows the rules for forming the nodal
For the capacitance C of a branch k, m (Fig. 3) the equation admittance matrix in steady-state analysis.
In (12) part of the voltages will be known (specified excit-
ekc (t) -em ()=J i,k,m (t ) dt + ek¢(t - At) - em (t - t ) ations) and the others will be unknown. Let the nodes be sub-
C t_^t divided into a subset A of nodes with unknown voltages and a
can again be integrated with the trapezoidal rule, which yields subset B of nodes with known voltages. Subdividing the matrices
and vectors accordingly, we get from (12)
ik,m (t) = (2C/At) (ek (t) - em (t) ) + 'k,m (t - At) (lOa)
LYAA][YAB] [eA (t)] iA (t)] IA]
with the equivalent current source Ik,m known from the past
history: [CYBA][YBB] [eB (t)] [iB (t0] _CIB]
Ik,m(t- At) = ik,m(t- At) from which the unknown vector EeA (t)] is found by solving
-
(2C/At) (ek(t At) - - em (t -
At)). (lOb) [YAA][eA (t)] = [Itota] - [YAB][eB (t)] (13)
An equivalent impedance network is shown in Fig. 3. Its with
form is identical with that for the inductance. The discretization EItotai]I= UiA (t )- EIA ].
error is also the same as that for the inductance.
This amounts to the solution of a system of linear equations in
Resistance each time step with a constant coefficient matrix [YAA], pro-
vided At is not changed. The right sides in (13) must be recalcu-
For completeness we add the branch equation for the resistance lated in each time step.
(Fig. 4):
ik,m(t) = (1/R) (ek(t) - em(t))* (11) Practical Computation
Equation (13) is best solved by triangular factorization of
Nodal Equations the augmented matrix [YAA], [YAB] once and for all before
With all network elements replaced by equivalent impedance entering the time step loop. The same process is then extended
networks as in Figs. 1-4, it is very simple to establish the nodal to the vector [Itotal] in each time step in the so-called forward
equations for any arbitrary system. The procedures are well solution, followed by back substitution to get [eA (t)], as indi-

Authorized licensed use limited to: ULAKBIM UASL - Yuzuncu Yil Universitesi. Downloaded on March 19,2021 at 08:04:28 UTC from IEEE Xplore. Restrictions apply.
DOMMEL: COMPUTER SOLUTION OF ELECTROMAGNETIC TRANSIENTS 391

cated in Fig. 5. Only a few elements in [YAA], [YAB] are non- READ DATA; SET INITIAL CONDITIONS;
zero; this sparsity is exploited by storing only the nonzero t1O
elements of the triangularized matrix. The savings in computer
time and storage requirements can be optimized with an ordered BUILD UPPER PART OF TRIANGULAR MATRIX
elimination scheme [10].
Should the nodes be connected exclusively via lossless lines, BUILD REDUCED MATRIX
with lumped parameters R, L, C only from nodes to datum, . I

then [YAA] becomes a diagonal matrix. In this case the equations |CHECK SWITCHES FOR CHANGES
could be solved separately node by node. Some programs are
based on this restricted topology. However, the sparsity tech-
T
HAVE ANY SWITCH POSITIONS CHANGED? -% -~~YES
nique lends itself automatically to this simplification without j NO
having to restrict the generality of the network topology. (IS THIS I S
The construction of the column vector ['total] is mainly a
bookkeeping problem. Excitations in the form of specified current
INOLALTER REDUCED MATRIX FOR SPECIFIC
sources [iA(t)] and the past history information in -[IA] are SWITCH POSITIONS; BUILD LOWER
PART OF TRIANGULAR MATRIX
entered into [Etotal] before going into the forward solution;
after [Itotal] has been built, using the still available voltages IF NONLINEAR ELEMENTS:
from the previous time step, specified voltage sources [eB (t)] FIND VECTOR I Z]
are entered into [e (t)]. Excitation values may be read from
cards step-by-step or calculated fromn standardized functions UPDATE PAST HISTORY, ENTER -[IA] INTO
(sinusoidal curve, rectangular wave', etc.). The excitations may Itt,nl AND OUTPUT DESIRED CURRENTS
be any combination of voltage and current sources, or there may
be no excitation at all (e.g., discharge of capacitor banks).
After having found EeA (t)], the past history records are updated
while constructing the vector [Etotal] for the next time step
(see flow chart in Fig. 6). Some practical hints about recording
the past history and about nonzero initial conditions may be
found in Appendixes I and II.
Approximation of Series Resistance of Lines
The simplicity of the method of characteristics rests on the
fact that losses are neglected. This simplicity also holds for the
distortionless line, where R'/L' = G'/C' (R' is the series re-
sistance and G' the shunt conductance per unit length); the
only difference is in computing Ik (and analogous Im):
Ik(t- T) = exp (- (R'/L')r) (- (1/Z)em(t- T) -im,k(t -T)). I LOWER PART OF TRIANGULAR MATRIX
*
Unfortunately, power lines are not distortionless, since G' is BACK S-UBSTITUTION; ej= ei IF SWITCH
usually negligible (or a very complicated function of voltage if j-i CLOSED AND i>j
corona is to be taken into account).
The distributed series resistance with G' = 0 can easily be iIFNONLINEAR ELEMENTS: CORRECT VOLTAGES
approximated by treating the line as lossless and adding lumped
resistances at both ends. Such lumped resistances can be inserted OUTPUT VOLTAGES
in many places along the line when the total length is divided
into many line sections. Interestingly, all cases tested so far Fig. 6. Flow chart for transienlts program.
showed no noticeable difference between lumped resistances
inserted in few or in many places. The voltage plot in Fig. 13
was practically identical for lumped resistances inserted in 3, with
65, and 300 places. In its present form, BPA's program auto- =(Z -4R)/ (Z + 4R)-
matically lumps R/4 at both ends and R/2 at the middle of the
line (R is the total series resistance); under these assumptions The real challenge for a better line representation is the fre-
the equivalent impedance network of Fig. 1 is still valid and quency dependence of R' and L', which results from skin effects
only the values cbange slightly (I,,, analogous to Ik): in the earth return [11] and in the conductors; BPA plans to
explore this further (see Section IV).
Z \= L'/ C' + 1 R
Ik(t - r) = ((1 + h)/2) {Ik from eq. (7b)} Switches
The network may include any number of switches, which may
+ ( (1 -h)/2) { Im from eq. (7b)} change their positions in accordance with defined criteria. They
are represented as ideal (R = 0 when closed and R = oo when
1 Since the trapezoidal rule is based on linear interpolation, a open); however, any branches may be connected in series or
rectangular wave of amplitude y will always be interpreted as having parallel to simulate physical properties (e.g., time-varying or
a finite rate of rise y/At in the first step in the presence of lumped
inductances and capacitances. current-dependent resistance).

Authorized licensed use limited to: ULAKBIM UASL - Yuzuncu Yil Universitesi. Downloaded on March 19,2021 at 08:04:28 UTC from IEEE Xplore. Restrictions apply.
392 IEEE TRANSACTIONS ON POWER APPARATUS AND SYSTEMS, APRIL 1969

SET A SET 8 EQ
EQ
WITHOUT

SWITCHES

LOWER PART O
TRIANGULA

Fig. 8. Solution for nonlinear parameter.

MATRIX SWITCHES MATRIX

(a) (b)
Fig. 7. Shaded areas show computation. (a) Initially. (b) After
each change.
I--

lr:
IYAA] [CZI] [JZ.] 4
Q
0
CD
et

0.
I z:r] [ Z "'
Fig. 9. Disconiiected subnetworks in [YAA1- Fig. 10. Influence of At.

With only one switch in the network, it is best to build the and the nonlinear equation in the form of the given characteristic,
matrix for the switch open and to simulate the closed position ek(t) em(t) f(ik,m(t)).
- =
(16a)
with superimposed node currents [2]. With more switches in the
network, it is preferable to build [YAA], [YAB] anew each time BPA's program represents the nonlinear characteristic point-by-
a change occurs. However, it is not necessary to repeat the point as piecewise linear (Fig. 8), but any mathematical function
entire triangular factorization with each change. Nodes with could be used instead.
switches connected are arranged at the bottom (Fig. 7). Then The nonlinear characteristic (16a) is that of a nonlinear,
the triangular factorization is carried out only for nodes without current-dependent resistance. -If it is to represent a lightning
switches (upper part of triangular matrix). This also yields a arrester, then ik,m = 0 until ek(lineaI) (t) - em(Iinear)(t)I reaches
reduced matrix for the nodes with switches (assumed to be the specified breakdown voltage of the arrester.
open). Whenever a switch position changes, this reduced matrix For a time-varying resistance, (16a) must be replaced by the
is first modified to reflect the actual switch positions (if closed: simpler equation
addition of respective rows and columns and retention of the em (t ) = R(lR ) * ik,m (t )
ek (t ) - (16b)
higher numbered node in place of two nodes), then the triangular
factorization is completed (lower part of triangular matrix). where R (tR) is given as a function of the time tR (e.g., in the
This scheme is included in the flow chart of Fig. 6. form of a table). The time count tR may be identical with the
time t of the transient study, or it may start later according to a
Nonlinear and Time-Varying Parameters defined criterion.
The characteristic of a nonlinear inductance is usually specified
With only one nonlinear parameter in the network, the so- as 41 = f (ik,m). The total flux is
lution can be kept essentially linear by confining the nonlinear
algorithm (usually an iterative procedure) to the branch with V (t) =
(ek (t) - em (t)) dt + (0).
the nonlinear parameter. To accomplish this the nonlinear pa-
rameter is not included in the matrix; its current ik,m is simulated With the trapezoidal rule of integration this becomes
with two additional node currents:
ek (t ) em (t) = (2/At)f (ik,m (t)) - c (t- At),
- (16c)
i. = ik,m and ik = ik,m
which simply replaces (16a). The value c must be updated with
Let [z] be the precalculated difference of the mth and kth
columns of [YAAJ-', which is readily obtained with a repeat c (0) = (2/At) <(0) + el(0) - em(0)
solution of (13) by setting [total] = {0, except +1.0 in mth and then recursively
and -1.0 in kth component} and [eB(t)] = 0. Ignoring the
nonlinear parameter at first, we get [eA(linear) (t)] from (13); c (t - At) = c (t - 2At) + 2 (ek (t - At) - e (t - At)).
the final solution follows from superimposing the two additional
currents ik = im = -ik,m: Generally, when a network contains more than one nonlinear
parameter, the entire problem becomes nonlinear and its iterative
[eA (t) ] = ] +EZ]. ikm (t).
EeA(linear) (t]) (14) solution quite lengthy. The algorithm remains simple, however,
if the network is topologically disconnected into subnetworks,
The value ikm in (14) is found by solving two simultaneous each containing only one nonlinear parameter. Disconnections
equations, the linear network equation (Thevenin equivalent) give [YAA] a diagonal structure with submatrices on the diagonal
) - emlinr) () + (Zk - Zm) ik,m(t)
(Fig. 9). Note that topological disconnections are quite likely
e ek(t)
- e ek
in networks containing lossless lines, since they, as well as
(15) lumped parameters from node to datum or to nodes with known

Authorized licensed use limited to: ULAKBIM UASL - Yuzuncu Yil Universitesi. Downloaded on March 19,2021 at 08:04:28 UTC from IEEE Xplore. Restrictions apply.
DOMMEL: COMPUTER SOLUTION OF ELECTROMAGNETIC TRANSIENTS 393

voltage, do not introduce off-diagonal elements into [YAA]. For

each subnetwork there is an independent equation (15) and each
vector [z] has zeros outside of that subnetwork. (Fig. 9 sym-
bolizes nonlinear parameters I-IV in four effectively disconnected
subnetworks.) Therefore each nonlinear parameter can be treated
separately and exactly as above. In the superposition (14),
each subnetwork will have its own [z] and ik,m. However, it is UM=
INk
L__ __ _ DATUM-
GROUND IN m
possible to compress these columns (Ezr]-[ziv] in Fig. 9) into
one single column, if an address column is added to indicate Fig. 11. Mutual couLpling.
the number of the subnetwork; the latter is necessary to insert
the right current ik,m into (14). BPA's program automatically [ik,m(t)] = [S]'([ek(t)]- [em(t)]) + [Ik,m(t - \t)] (17a)
checks for violations of the disconnection rule while computing
this single column together with the address column. with [Ik,m (t - At)] from the recursive formula
[Ikm (t -At)] = [H]([ek (t- At)]- [em (t -At)]
Accuracy
+ [S][lkm (t - 2A t)]) - [lk,r(t - 2At)]. (17b)
To arrive at (13), approximations have to be made only for
lumped inductances and capacitances. Lossless lines and re- All matrices in (17) are symmetric:
sistances are treated rigorously. [S] = [R]+ (2/At)[L]
In practice, a truncation error is introduced by a lossless line
whenever its travel time -r is not an integer multiple of At. [H]= 2([S]1-- [S]-E[R][S]-').
Then some kind of interpolation becomes necessary in computing The only difference compared with a single branch is, that in
Ik(t - r) and Im(t - -). One option in BPA's program uses building [YAA], [YAB] in (13), a matrix [S]-' is entered instead
linear interpolation, because in most practical cases the curves of a scalar value. Also in each time step a vector Elk,i,] enters
e(t) and i(t) are smooth rather than discontinuous. For cases into [Itotal] instead of a scalar value.
with expected discontinuities, another option rounds the travel If Fig. 11 is part of a multiphase 7r-equivalent representing
time r to the nearest integer multiple of At. Both options raise a line section, then each set of terminals will be capacitance
travel times r < At to At; otherwise the equivalent impedance connected. These capacitances are actually single branches; thus
network of Fig. 1 could not be used any more. no new formula is necessary. BPA's program treats them as a
The trapezoidal rule of integration, used for lumped in- matrix entity [C] to speed up the solution.
ductances and capacitances, is considered to be adequate for
practical purposes, especially if the network has only a few Lossless Multiphase Line
lumped parameters. Compared with the alternative of stubline
approximations [9], the results are more accurate. It is well Equation (1) is also valid for the multiphase line if the scalars
known that the trapezoidal rule is numerically stable and has are replaced by vectors [e], [E] and matrices EL'], EC']. By
almost ideal round-off properties [12, p. 119]. When the step- differentiating a second time, one of the vector variables can be
width At is chosen sufficiently small to give good curve plots eliminated, which gives
(points not spaced too widely), linear interpolation, on which
the trapezoidal rule is based, should be a good approximation. E[2e(x, t)/&x2] = [L'][C']I[2e(x, t)/9t2] (18a)
Both requirements go hand in hand. The choice of At is not [a2i(X, t)/aX2] = EC'][L'][02i(x, t)/0t2]. (18b)
critical as long as the oscillations of highest frequency are still
represented by an adequate number of points. Changing At The solution of (18) is complicated by the presence of off-
influences primarily the phase position of the high-frequency diagonal elements in the matrices, which occur because of mutual
oscillations; the amplitude remains practically unchanged (see couplings between the phases. This difficulty is overcome if the
Fig. 10 which resulted from an example similar to that of Fig. 12). phase variables are transformed into mode variables by similarity
Higher accuracy could be obtained with the Richardson ex- transformations that produce diagonal matrices in the modal
trapolation [12, p. 118]. Here, the integration from (t - At) equations [2], [13], [14]. This is the well-known eigenvalue
to t would be carried out twice, with At in one step and with problem. Each of the independent equations in the modal domain
At/2 in two steps, and both results extrapolated to At = 0. can then be solved with the algorithm for the single-phase line
The amount of work in each time step is thus tripled and the by using its modal travel time and its modal surge impedance.
work for the initial triangular factorization is doubled, since two The transformation matrices, which give the transition to the
matrices, built for At and At/2, are necessary. phase domain, will generally be different for voltages and cur-
rents, e.g.,
III. MUTUAL COUPLING AND MULTIPHASE NETWORKS [epha,s] = [Te][emode] (19a)
Lumped Parameters with Mutual Coupling [iphae] = [Ti]j[ode]. (19b)
To include mutual coupling with lumped parameters the The columns in [T.], [Ti] are always undetermined by a con-
scalar quantities of a single branch are simply replaced by matrix [2], [15] is:
stant factor, if not normalized. A helpful relation
quantities for the set of coupled branches. Consider the three
coupled branches in Fig. 11 with a resistance matrix [R] and [Tijunnormaiized = [C'][T]. (19c)
an inductance matrix EL]. They could represent the series
branches of a three-phase r-equivalent with earth return; in If all diagonal elements in [L'] are equal to L'^lf and all off-
this case [L] as well as [R] would have off-diagonal elements diagonal elements are equal to L'mutuai (analogous for [C']),
(mutual coupling). Applying the trapezoidal rule of integration then a simple transformation is possible, even if the inductances
[2] yields: are frequency dependent [15]:

Authorized licensed use limited to: ULAKBIM UASL - Yuzuncu Yil Universitesi. Downloaded on March 19,2021 at 08:04:28 UTC from IEEE Xplore. Restrictions apply.
394 IEEE TRANSACTIONS ON POWER APPARATUS AND SYSTEMS, APRIL 1969

I I

180-MILE LINE
0- -89
C~~~~~~~~~~~
J A~~~~~~~~
_
18 MULTI - TT-EQUIVALENTS
SWITCH CLOSING: I 1
A lOims
B 4ms 14 ms
C 4ms 14ms

Fig. 12. Sequential closing. Network and results at the receivin-g end. Line energizing: 180-mile line,
transposed at 60 and 120 miles. RLC for 60 Hz.

a single line. Work is in progress at BPA to incorporate the

frequency dependence approximately into the method of char-
[Te] = [Ti]= 'i I.. 1 (20) acteristics; then, instead of one value from the past history,
several weighted samples will go into the computation of Ik
and Im. The weights would have to be chosen to match the
frequency spectrum derived from Carson's formula [11] or from
where M is the number of phases. measurements on the line. In a similar approach [18], the earth
It can be shown [2] that the phase current vector [ik,m] return mode is passed through two RC filters before entering
entering the nodes at terminal k toward m can again be written the node, while the others are attenuated without distortion.
as a linear vector equation
V. EXAMPLE S
[ik,m (t) ] = [G][ek (t ) ] + [Ik]
Two simple cases are used to illustrate applications of the
and analogous for [i4,k]. Equation (21) is derived from a set program. Fig. 12 shows the results for sequential closing of a
of modal equations, subjected to the transformations (19). In three-phase, open-ended line. The curves were automatically
building EYAA], [YAB] in (13), a matrix EG] is entered instead plotted by a Calcomp plotter. For this study, the line was
of a scalar value 1/Z. The vector ElI], which enters [IA], is represented by 18 multi-r-equivalents with (coupled) lumped
calculated from the past history of the modal quantities. Since parameters. Fig. 13 shows the voltage at the receiving end of a
the span (t - ') for picking up the past is different for each single-phase line (320 miles long, R' = 0.0376 Q/mi, L' = 1.52
mode, a time argument was deliberately omitted in writing mH/mi, C'= 0.0143 ,F/mi), that is terminated by an induct-
Elk]. Even though the nodal equations are in phase quantities, ance of 0.1 H and excited with a step function e (t) = 10 V.
the past history must be recorded in modal quantities. The solid curve results from representing the line with 32
lumped-parameter equivalents, the dashed curve from a dis-
IV. FREQUENCY-DEPENDENT LINE PARAMETERS tributed-parameter representation.
Skin effects in the earth return and conductors make the line VI. CONCLUSIONS
parameters R' and L' frequency dependent [11], [14]. In
multiphase lines, this affects primarily the mode associated with A generalized digital computer method for solving transient
earth return. It is not easy to take the frequency dependence phenomena in single- or multiphase systems has been described.
into account and at the same time maintain the generality of the The method is very efficient and capable of handling very large
program. Methods using the Fourier transform [15], [16] or networks. Further work is necessary to find a satisfactory way
the Laplace transform [17] are usually restricted to the case of to represent frequency dependence of line parameters.

Authorized licensed use limited to: ULAKBIM UASL - Yuzuncu Yil Universitesi. Downloaded on March 19,2021 at 08:04:28 UTC from IEEE Xplore. Restrictions apply.
DOMMEL: COMPUTER SOLUTION OF ELECTROMAGNETIC TRANSIENTS 395

eliminating the inner nodes in the connection [2]. Likewise,

single- or multiphase iz-equivalents with series [R] and [L]
matrices and with identical shunt [C] matrix at both terminals
are treated as one element. If the system has identical network
elements (e.g., in a chain of wr-equivalents), then the data are
specified and stored only once.
APPENDIX II
INITIAL CONDITIONS
BPA's computer program has two options for setting nonzero
initial conditions. Voltages and currents at any point in a study
can be stored and used again as initial conditions in subsequent
studies that take off from that point (usually with a different
At). They can also be computed for any sinusoidal steady-state
condition with a subroutine "multiphase steady-state solution."
The first option must be used if the steady-state solution is
Fig. 13. Single-phase line with inductive termination. nonsinusoidal because of nonlinearities. In this case a transient
study is made once and for all over a long enough time span to
settle to the steady state. This gives initial conditions for all
APPENDIX I subsequent studies.
RECORDING THE PAST HISTORY
ACKNOWLEDGMENT
The equivalent current sources I in Figs. 1-3 constitute that
part of the past history, known from preceding time steps, that The author wants to thank his colleagues at the Bonneville
has to be recorded and constantly updated. They are needed Power Administration, notably Dr. A. Budner, J. W. Walker,
in building the vector ['total]. For each inductance and ca- and W. F. Tinney, for their help and for their encouragements.
pacitance a single value Ik,m (t - At) must be recorded, for each The idea of weighted samples to incorporate the frequency
lossless line a double list Ik, I. for the time steps t - At, t - 2At, dependence of line parameters is due to Dr. A. Budner, and the
... ,t-r. subroutine to get ac steady-state initial conditions was written
In updating Ik,m for inductances an-d capacitances, it is faster by J. W. Walker.
to use recursive formulas:
REFERENCES
Ik,m (t - At) = (Ikt,m (t - 2At) + 2x) [1] 11. Prinz and H. Dommel, "tberspanuiungsberechinung in
Hochspannungsnetzen," presented at the Sixth Meeting for
(+ for inductance, - for capacitance), with x = G (ek (t -At)- Industrial Plant Managers, sponsored by Allianz Insurance
em (t At)) and G = At/2L for inductance and G = 2C/At for
- Company, Munich, Germany, 1964.
capacitance. [2] H. Dommel, "A method for solving transient phenomena in
multiphase system," Proc. 2nd Power System Computation
These formulas are easily verified by expressing the currents Conference 1966 (Stockholm, Sweden), Rept. 5.8.
in (9b) and (lOb) by (9a) and (lOa), respectively. To assure [3] F. H. Branin, Jr., "Computer methods of network analysis,"
Proc. IEEE, vol. 55, pp. 1787-1801, November 1967.
correct initial values in the very first time step, Ik,m must be [4] L. Bergeron, Du Coup de Belier enHydraulique au Coup deFoudre
preset before entering the time step loop en Electricite. Paris: Dunod, 1949. Transl., Water Hammer in
Hydraulics and Wave Surges in Electricity (Translating Commit-
Ikcm (initial) = i%m (O) -G (e (O) - em (0)). tee sponsored by ASME). New York: Wiley, 1961.
[5] H. Prinz, W. Zaengl, and 0. Vdleker, "Das Bergeron-Verfahren
zur Loesung von Wanderwellen," Bull. SEV, vol. 16, pp. 725-
The initial conditions e (0) for voltages and i (0) for currents are 739, August 1962.
part of the input. [6] W. Frey and P. Althammer, "Die Berechnung elektromag-
For a lossless line the values Ik, 'm must be recorded for netischer Ausgleichsvorgaenge auf Leitungen mit Hilfe eines
Digitalrechners," Brown Boveri MItt., vol. 48, pp. 344-355, 1961.
t- At, t 2At, - ...back to t r; they are stored in one
, -
[71 P. L. Arlett and R. Murray-Shelley, "An improved method
double list, where the portion for each line has its length adjusted for the calculation of transients on transmission lines using a
to its specific travel time r. After [e (t)] has been found, the
digital computer," Proc. PICA Conf., pp. 195-211, 1965.
[81 F. H. Branin, Jr., "Transient analysis of lossless transmission
double list is first shifted back one time step (entries for t - At lines," Proc. IEEE, vol. 55, pp. 2012-2013, November 1967.
become entries for t 2t, etc.); then Ik/(t - r), Im(t - r) are
-
[91 L. 0. Barthold and G. K. Carter, "Digital traveling-wave solu-
tions," AIEE Trans. (Power Apparatus and Systems), vol. 80,
computed and entered into the list. Physically, the list is not pp. 812-820, December 1961.
shifted; instead, the starting address is raised by 1 modulo [10] W. F. Tinney and J. W. Walker, "Direct solutions of sparse
{length of double list} [8]. The initial values for 'k, Im must network equations by optimally ordered triangular factoriza-
ation," Proc. IEEE, vol. 55, pp. 1801-1809, November 1967.
be given for t = 0, -At, -2A1, -r. The necessity to know
... , [11] J. R. Carson, "Wave propagation in overhead wires with ground
them beyond t = 0 is a consequence of recording the terminal return," Bell Syst. Tech. J., vol. 5, pp. 539-554, 1926.
conditions only. If the conditions were also given along the line [12] A. Ralston, A First Course in Numerical Analysis. New York:
McGraw-Hill, 1965.
at travel time increments At, then the initial values at t = 0 [13] A. J. McElroy and H. M. Smith, "Propagation of switching-
would suffice. surge wavefronts on EHV transmission lines," AIEE Trans.
(Power Apparatus and Systems), vol. 81, pp. 983-998, 1962
BPA's computer program has features that help to speed up (February 1963 sec.).
the solution. Thus a series connection of resistance, inductance, [14] D. E. Hedman, "Propagation on overhead transmission lines
and capacitance is treated as a single branch. This reduces the I-theory of modal analysis," IEEE Trans. Power Apparatus
and Systems, vol. PAS-84, pp. 200-211, March 1965; discussion,
number of nodes; the respective formulas can be derived by pp. 489-492, June 1965.

Authorized licensed use limited to: ULAKBIM UASL - Yuzuncu Yil Universitesi. Downloaded on March 19,2021 at 08:04:28 UTC from IEEE Xplore. Restrictions apply.
396 IEEE TRANSACTIONS ON POWER APPARATUS AND SYSTEMS, APRIL 1969

[15] H. Karrenbauer, "Ausbreitung von Wanderwellen bei ver-

schiedenen Anordnungen von Freileitungen im Hinblick auf
die Form der Einschwingspannung bei Abstandskurzschlues-
sen," doctoral dissertation, Munich, Germany, 1967.
[16] M. J. Battisson, S. J. Day, N. Mullineux, K. C. Parton, and
J. R. Reed, "Calculation of switching phenomena in power
systems," Proc. IEE (London), vol. 114, pp. 478-486, April
1967; discussion, pp. 1457-1463, October 1967.
[17] R. Uram and R. W. Miller, "Mathematical analysis and solu- L1 1n 1 1,1 I' I I s
tion of transmission-line transients I-theory," IEEE Trans.
-

Power Apparatus and Systems, vol. 83, pp. 1116-1137, Novem- R,=Or) -SOns R,=1000 -SOns
ber 1964.
[18] A. I. Dolginov, A. I. Stupel', and S. L. Levina, "Algorithm
[-J
L... ["J.

ET T
and programme for a digital computer study of electromagnetic
transients occurring in power system" (in Russian), Elektri- _ I I I_' " I I EI:]=--
chestvo, no. 8, pp. 23-29, 1966; English transl. in Elec. Technol. R,=3000 -SOns Rd=500S , 50 ns

(USSR), vols. 2-3, pp. 376-393, 1966. .6000-

mercury relay pie 5 2,32 kO- d,ider

.n,rntnnn 4-sc.
u

Fig. 14. Measured step response of a low-impedance voltage divider.

Discussion

W. Zaengl and F. W. Heilbronner (Hochspannungsinstitut der

Technischen Hochschule Muinchen, Munich, Germany): Dr. Dommel
is to be congratulated for these lucid elaborations of the treatment
of electromagnetic transients. In order to demonstrate how effective
this method is, we wish to append two examples of a single-phase
application of the algorithm as described and the verification by
experiments: 1) evaluation of the step response of an impulse voltage
measuring circuit and 2) computation of the voltage breakdown in
sparkgaps.
1) In high-voltage measuring techniques voltage dividers are used
which cannot be constructed coaxially and are, because of voltages
up to some million volts, of big dimensions. Therefore the voltage to
be measured is led to the divider by metallic pipes, at the input end
of which, in general, a damping resistor is connected.
For this purpose the equivalent circuit of the total measuring Fig. 15. Calculated step response of the test setup
circuit is best represented by a lossless line (for the metallic pipe), according to Fig. 14.
on which traveling wave phenomena occur, and lumped parameters
(for damping resistor and voltage divider). An analytic general
solution to get the step response of this network is not possible.
In Fig. 14 the used measuring circuit is sketched with its dimen-
sions. The 2.32-kU divider consists of stacked resistors. The output
voltage, reduced by a factor of 100, is measured by an oscilloscope
(Tektronix 585). Four oscillograms of the output voltage are given,
resulting from various damping resistors Rd, if a voltage step gener-
ated by a mercury relay occurs at the input end of the measuring
circuit.
In Fig. 15 the equivalent circuit of the test setup with its data is
given and the results of the digital computation of the step response
G(t) with the program outlined in the paper. The surge impedance
Z = 272 ohms and the travel time r = 20 ns result from the geo-
metric dimensions of the pipe. The divider is represented by a multi-
section network of a total of five T quadripoles and an input shunt
capacity C, = 5 pF. In the calculation a step width At of 2-10-
seconds was used. The comparison shows a very good agreement with
the experimental results of Fig. 14.
2) Whereas the solution of the foregoing problem requires no
specific modification of the straightforward procedure as described
in the paper, in the case of voltage breakdown, nonlinearities have to
be taken into account [19]. One means of evaluating the voltage u
at a time t during breakdown of a gap was given by Toepler [20]:
u(t) = kla-i(t) /f i(t) df (22) Fig. 16. Test setup. Front left: screened measuring cabin; front
center: damied capacitive divider; front right: 80-cm rod-rod gap;
i.e., the resistance of the spark is inversely proportional to the center: 3-million-volt impulse generator (capacitive divider is used
amount of charge which has flowed into the gap (a = gap spacing as load capacitance and is standing in front of the generator).
in cm, k = constant in the range of 10-4 V-s/cm, i(t) = current in
amperes, t = time in seconds).

Manuscript received July 3, 1968.

Authorized licensed use limited to: ULAKBIM UASL - Yuzuncu Yil Universitesi. Downloaded on March 19,2021 at 08:04:28 UTC from IEEE Xplore. Restrictions apply.
 DOMMEL: COMPUTER SOLUTION OF ELECTROMAGNETIC TRANSIENTS 397

Equivalent Circuit of Load Capacitance Rod-Rod Gap Damped-Capacitive

3MY Impuls Generator Used as Voltage Spaced 80 cm Divider
Divider Lower Electrode
-3021) 1 m Above Ground.
(Uc.d
Shunt for Current (ud.c.d.=4000)
Measuremen t.
Fig. 17. Equivalent impulse circuit of Fig. 16.

O 2jus

(a)

(b)

1~~~~~~~~ . . .. .. .+.R

(c)

Fig. 18. Oscillograms from the voltage breakdown of a 80-cm rod-rod Fig. 19. Calculated voltages in different points and gap current at
gap (temperature: 20'C, 716 mm Hg); horizontal deflection 10-6 breakdown according to Fig. 17. (a)-voltage of capacitive divider;
seconds/division. (a) Capacitive divider: 138 kV per vertical (b)-gap current j(t); (c)-gap voltage u(t.
division. (b) Damped capacitive divider: 183 kV per division.
(c) Current shunt: 1060 amperes per vertical division.

Using the trapezoidal rule of integration, (22) can be rewritten as

Thus, in terms of the paper, the voltage across a sparkgap between
int( -) a((-At) )) At(23) the nodes k and m will be
R(tIAl) (ek (t) em (t) ) (linear)
where int(I-At) is the value of the integral in the denorninator of (22) (25)
R(t-t) + (Z,k Zk,m Zm,k + Zm,m)
At). This is the equivalent expression Off (ik,,m (t) )
- -

up to the time (t -

in (16a). Since the solution in connection with (15) would be of the In order to start the process, in (24) a certain initial value of int(1-1
quadratic type, it was found sufficient to linearize the problem and is needed. This means in physical terms, that by some predischarges
take the resistance of the previous time step (t - At): the gap must have been ionized and thus assumed some conductivity.
Experience has shown that for a start the value of R(t-,t) might be
chosen a thousand times higher than the biggest resistance in the
int(R-2( . (j att) + j(t-AO) .24 circuit.

Authorized licensed use limited to: ULAKBIM UASL - Yuzuncu Yil Universitesi. Downloaded on March 19,2021 at 08:04:28 UTC from IEEE Xplore. Restrictions apply.
398

As a demonstration, in Fig. 18 three oscillograms (Tektronix 507)

of the breakdown of a 80-cm rod-rod gap in the test circuit of Figs. 16
and 17 are given. The calculated values with k = 0.3. 10-4 Vs/cm
and At = 20@10-9 seconds, multiplied by the corresponding divider
ratios, are plotted in Fig. 19. They correspond fairly well to the
oscillograms. The voltage resulting from the damped capacitive
divider is within +41 percent of u and is therefore not plotted.
Two conclusions may be drawn from a comparison of Figs. 18
and 19 and are stated without further explanation: 1) A damped
capacitive divider [21] reflects the gap voltage much better than a
purely capacitive divider, and 2) the common equivalent circuit of
:
IEEE TRANSACTIONS ON POWER APPARATUS AND

-ROW K BUILT
SYSTEMS, APRIL 1969

PART OF MATRIX
ALREADY TRIANGULARIZED
AND STORED IN PACKED
FORM

FROM BRANCH TABLE

Fig. 20. Triangularization scheme.

a divider may be too rough in the cases where higher harmonics
occur. Then an equivalent circuit as in item 1) would be necessary.
The described application of the transient algorithm in high- H. W. Dommel: The author wants to thank Dr. XV. Zaengl, Mr. F.
voltage impulse circuits has led the discussers to various secondary W. Heilbronner, Mr. D. G. Taylor, and Mr. M. R. Payne for their
problems and suggestions, of which two can be sketched here in valuable discussions, which illustrate the usefulness of Bergeron's
general terms only. method in traveling wave studies and also raise some interesting
1) In problems with many nodes, computer storage might be too questions.
small for building up the matrix [Y]. Thus the method of diacoptics One of the main differences between the author's computer pro-
is of help, especially when two major parts of the circuit are con- gram and that of Mr. D. G. Taylor and Mr. M. R. Payne is the
nected by a single lead which can be represented by a lumped subdivision of the line into sections of travel time r = At in the latter.
parameter (inductance in Fig. 17). It appears that considerable savings in computer time (but not in
2) If sudden changes of network parameters occur, e.g., the break- storage requirements) are possible when such subdivisions are
down of a sparkgap on account of a certain overvoltage, where the avoided and multiple past histories are stored. It must be admitted,
resistance changes from the order of megohm to ohm in fifty to some however, that lumped series resistances can be included more easily
hundred nanoseconds, it might be desirable to make the time step in more places with the line being subdivided, even though this can
At smaller and increase it again when the rate of change is no longer always be done with the author's program in the definition of the
of importance. Thus it is necessary to adapt the stepwidth At to the model at the expense of more input data. Interestingly enough, test
rate of voltage change in the network. examples showed very little or no difference at all between the inser-
tion of lumped series resistances in few or many places (section
REFERENCES Approximation of Series Resistance of Lines). Therefore, the auto-
[19] F. Heilbronner and H. Kiirner, "Ein Verfahren zur digitalen matic insertion at three places (terminals and midpoint) was felt
Berechnung des Spannungszusammenbruchs von Funken- to be adequate. This observation might not be true for all cases.
strecken," ETZ-A, vol. 89, pp. 101-108, 1968. Also, not too much significance has been placed on the approximation
[20] M. Toepler, "Funkenkonstante, Ztindfunken und Wander- of distributed resistance by lumped, series resistances in developing
welle," Arch. f. Elektrotech., vol. 16, pp. 305-316, 1925. the program, since the final objective has been the approximation of
[21] W. Zaengl, "Das Messen hoher, rasch verhnderlicher Stoss- the frequency dependence in the zero-sequence mode. This has not
spannungen," doctoral dissertation, Munich, 1964. been included yet, but preliminary tests with a weighting function
representation look promising.
Mr. D. G. Taylor and Mr. M. R. Payne use a stub-line (short
line) representation for lumped (series and shunt) L and shunt C.
It can be shown that this stub-line representation for shunt L and
D. G. Taylor and M. R. Payne (Central Electricity Generatiiig shunt C is equivalent with the integration of (8b), and the respective
Board, London, England): We have also programmed the Bergeron equation for C, by the trapezoidal rule over two time steps from
method for single-phase switching problems and are currently en- t- 2At to t (no such simple equivalence was found for series L).
gaged in extending the treatment to multiconductor systems. Since the author's method for lumped L and C is based on the
Lumped L and shunt C have been represented as short lines and trapezoidal rule of integration over one time step only from t - At
special "hyper-nodes" have been introduced to deal with series R to t, it is more accurate than stub lines. The stub-line representation
and series C. Only one past history is stored which necessitates sub- is very helpful, however, in studies involving more than one non-
dividing lines into sections of equal traveling time. Processed system linear element. As described in the section Nonlinear and Time-
data together with past and present values of voltage and current Varying Parameters, more than one nonlinear element can be handled
are stored in a structured file (in core) which is passed, using list- in closed form only if they are separated by elements of finite travel
processing techniques, in order to advance the solution by one time time. A stub-line representation accomplishes just such a separation.
step. As an example, a case involving a lightning arrester connected to a
One advantage of subdividing lines over storing multiple past nonlinear inductance (transformer with saturation) can be solved
histories is that series resistance can be introduced between all by modeling the total inductance as a linear and nonlinear inductance
sections; we have found this to be desirable in cases where the in series, with the linear inductance placed on the side of the lightning
response is oscillatory and the degree of attenuation is important. arrester and treated as a stub line.
The author's comments on this point would be appreciated. Mr. D. G. Taylor and Mr. M. R. Payne raise the question of errors
A source of approximation which should be mentioned arises from introduced either by making all travel times an integral multiple of
the necessity for all traveling times to be integral multiples of the At or by using interpolation between past values. It is true that this
time increment At. This also applies to the method of multiple past question is even more critical in multiconductor systems with small
histories since any interpolation between values is invalid. The differences in mode propagation velocities. Interpolation is indeed
problem is made more severe in multiconductor systems by the questionable if sudden changes occur. However, the presence of
propagation velocities in the modes being different, in some cases inductances and capacitances often, though not always, smoothes
by small but significant amounts. How does the author take this out sudden changes; then interpolation is a good and valid approx-
into account in making his initial choice of time increment, in imation. Sudden changes may also be introduced through stub-line
particular for systems including asymmetrical multiconductor con- representations and not lie in the nature of the problem. In cases
figurations? where sudden changes do occur, the user has an option in which all
In conclusion the author is to be congratulated on his adaptation travel times are rounded to the nearest integral multiple of At. As
of the problem for use with ordered-elimination techniques which of now, the step width At must be chosen by the user.
have already made such an impact on steady-state analysis. We look In the first part of their discussion, Dr. W. Zaengl and Mr. F. W.
forward to the author's further developments in this field, particu- Heilbronner show how closely computed results can agree with test
larly with regard to the treatment of frequency dependence. results. This speaks at least as much for their good engineeririg

Manuscript received July 3, 1968. Manuscript received August 8, 1968.

Authorized licensed use limited to: ULAKBIM UASL - Yuzuncu Yil Universitesi. Downloaded on March 19,2021 at 08:04:28 UTC from IEEE Xplore. Restrictions apply.
 IEEE TRANSACTIONS ON POWER APPARATUS AND SYSTEMS, VOL. PAS-88, NO. 4, APRIL 1969 399

judgment in selecting an equivalent model as it does for the useful- the triangularized matrix, and finally the elements Y'k,k, Y'k k+k+, -
ness of the computer program. Their effort to include the dynamic of this transformed row are added in packed form to the triangular-
law of spark gaps into the program should be of interest to high- ized matrix. In a way, the method does have a built-in tearing
voltage engineers. feature similar to diacoptics in cases involving lines with distributed
As to the specific questions raised, it is felt that the sparsity parameters, which disconnect the network topologically. This dis-
technique used (optimally ordered elimination with packed storage connection is more than tearing in diacoptics, since it is a true dis-
of nonzero elements only) is more efficient than the method of connection where no reconnection effect has to be introduced at a
diacoptics. It was probably not made clear in the paper that the later stage of the algorithm. Thus, the use of a stub-line representa-
matrix [Y] is never built explicitly. Rather, a branch table is used tion for the inductance in Fig. 17 with surge impedance Z = L/At
to store the information for the matrix [Y]. As indicated in Fig. 20 and travel time r = At, might reduce the storage requirements
for the kth elimination step, the original row k is built from a search beyond those already achieved through sparsity. The possibility to
of the branch table (therefore, only one working row is necessary), change At during the computation would indeed be desirable. It is a
then the elements to the left of the diagonal are eliminated with the straightforward programming task, involving changes of [Y]. Due
information contained in the already available rows 1, *, k - 1 of to lack of time, it has not been incorporated so far.

Nonlinear Programming Solutions for Load-Flow,

Minimum-Loss, and Economic Dispatching
Pro bles
ALBERT M. SASSON, MEMBER, IEEE

Abstract-A unified approach to load-flow, minimum-loss, and The load-flow problem [1], [2] was first solved by a simplified
economic dispatching problems is presented. A load-flow solution Newton-Raphson approach which involved the power system
is shown to coincide with the minimum of a function of the power nodal admittance matrix. As the equations are not quadratic,
system equations. An unconstrained minimization method, developed the simplified approach together with an iterative process was
by Fletcher-Powell, is used to solve the load-flow probl m. The
method always finds a solution or indicates the nonexistence of a justified. Later [3] an iterative Gauss-Seidel approach was
solution. Its performance is highly independent of the reference- successfully used. Further improvements were based on using
slack bus position and requires no acceleration factors. Several con- the nodal impedance matrix [4] and the mesh impedance
strained minimization techniques that solve the minimum-loss and matrix [5]. More recently Newton's technique has been used
economic dispatching problems are investigated. These include the [6], [7] claiming extremely rapid convergence. Even if much
Fiacco-McCormick, Lootsma, and Zangwill methods. The technique progress has been made in load-flow analysis, there are situations
finally recommended is shown to be an extension of the method which cause difficulties in obtaining solutions with some of these
used to solve the load-flow problem. The approved IEEE test sys- methods. The position of the reference-slack bus, the choice of
tems, and other systems whose response to conventional methods acceleration factors, the existence of negative line reactances, a
was known, have been solved.
large ratio of long-to-short line reactance for lines terminatinig
in the same bus, and certain types of radial systems are the
INTRODUCTION cause of much instability in methods of solution of the load-flow
problem. When a divergent solution is obtained, it is not clear
whether the divergence has been due to instability in the method
UCH WORK has been done in the fields of load-flow used or to the fact that there may not be a solution at all. The
analysis and economic dispatching; some papers have pre- methods presented in this paper are quite insensitive to many
sented methods that obtain a minimum-loss solution. Each of of the factors which cause instability to existing methods and
these problems has been solved independently from the others. give a definite answer as to whether a solution exists or not.
The methods discussed in this paper present a unified approach The approach is based'on the construction of a function of the
which demonstrates that all three problems fall into a single power-system equations, whose minimum coincides with the
class of optimization problems. solution of the equations.
The minimum-loss problem has also been solved in various
Paper 68 TP 673-PWR, recommended and approved by the ways. The first approach to the problem was to solve the eco-
Power System Engineering Committee of the IEEE Power Group nomic dispatching problem minimizing losses at the same time.
for presentation at the IEEE Summer Power Meeting, Chicago,
Ill., June 23-28, 1968. Manuscript submitted February 7, 1968; This ignored the possibility of minimizing losses by an optimal
made available for printing May 14, 1968. use of the reactive capabilities of the system as a whole. It is
The author is with the Imperial College of Science and Tech- from the second point of view that this paper considers the
nology, London, England, and the Instituto Teenologico y de
Estudios Superiores de Monterrey, Monterrey, N. L., Mexico. minimum-loss problem. One of the first attempts [8] was to

Authorized licensed use limited to: ULAKBIM UASL - Yuzuncu Yil Universitesi. Downloaded on March 19,2021 at 08:04:28 UTC from IEEE Xplore. Restrictions apply.