Electrical Power and Energy Systems 72 (2015) 24–32

Contents lists available at ScienceDirect

Electrical Power and Energy Systems

journal homepage: www.elsevier.com/locate/ijepes

A discussion about optimum time step size and maximum simulation

time in EMTP-based programs
José Carlos G. de Siqueira a, Benedito D. Bonatto a,⇑, José R. Martí b, Jorge A. Hollman c,
Hermann W. Dommel b
a
UNIFEI – Federal University of Itajuba, Itajuba, MG, Brazil
b
UBC – The University of British Columbia, Vancouver, BC, Canada
c
BC Hydro, BC, Canada

a r t i c l e i n f o a b s t r a c t

Article history: This paper presents a discussion about the optimum time step size and maximum simulation time in
Received 1 February 2015 EMTP-based programs. This is of particular importance for new users of EMTP-based programs, since
Accepted 16 February 2015 the user is responsible for setting up these parameters before running a simulation case. The selection
Available online 26 March 2015
of the time step size affects the precision of the simulation. The time step size depends on the maximum
frequency expected in the phenomena, which is normally unknown, a priori. The calculation of eigenval-
Keywords: ues from the state space matrix obtained from the EMTP conductance matrix and history terms algorithm
EMTP
formulation allows the exact ideal values for the optimal time step size and a maximum simulation time.
Time step size
Time constant
In comparison, a general and robust algorithm is also presented here based on all the input data given for
Eigenvalue the circuit under simulation. The proposed calculation process is based on single- or multi-phase uncou-
Computer modeling pled or coupled circuits, with lumped or distributed parameters. Simulations are given demonstrating the
Transients effectiveness of the proposed rules.
Ó 2015 Elsevier Ltd. All rights reserved.

Introduction The methods described in this paper, follow the traditional

EMTP approach of using a constant time step calculated according
The existing versions of the EMTP software [1–9] do not provide to the maximum frequencies that need to be simulated accurately.
the user with an optimal estimate of the time step size ðDtÞ to be The trapezoidal integration rule is used in EMTP-based software,
used for a given simulation case. It is up to the user to choose ‘‘because it is simple, numerically stable, and accurate enough for
and adjust it properly. The same situation happens for the practical purposes’’ [1]. The effective suppression of numerical
maximum simulation time ðtmax Þ. oscillations in the EMTP, which are not solved by just reducing
Even though there is a close relationship between the time step the time step size, is well discussed in [13,14], with the Critical
size, the time constants, the propagation time, and the oscillations Damping Adjustment (CDA) method. The evaluation of accuracy
periods of voltages and currents in a simple circuit, this relation- and stability of different integration rules (integration solvers) as
ship is difficult to verify for large systems, because usually there functions of the time step size (or of the sampling frequency, or
is no prior knowledge of the analytical solution. of the per unit of the Nyquist frequency) is clearly discussed in
Programs like SPICE [10] change the time step continuously [13]. For an overall accurate and stable simulation, CDA takes
during the simulation, according the slope of change of the circuit advantage of the accuracy of the trapezoidal integration rule and
variables. This strategy, however, can result in increased solution the stability of the Backward Euler integration rule, for a given
times and in some cases in non-convergence [11], or numerical and fixed time step size.
instability conditions. Variable step size ODE solvers are discussed Unless the circuit equations are solved analytically (or the
in [12] and are beyond the scope of this paper. eigenvalues are somehow calculated for linear systems or for
linearized parts of nonlinear systems [15,16]), the time simulation
parameters are not known by inspection. It is possible, however, to
⇑ Corresponding author. establish ‘‘ranges’’ in which the time step size and the maximum
E-mail addresses: profjcgoulart@gmail.com (J.C.G. de Siqueira), bonatto@unifei. simulation time must be contained, as presented in [17,18]. This
edu.br (B.D. Bonatto), jrms@ece.ubc.ca (J.R. Martí), jorge.hollman@bchydro.com approach may help new and experienced EMTP-based program
(J.A. Hollman), hermannd@ece.ubc.ca (H.W. Dommel).

http://dx.doi.org/10.1016/j.ijepes.2015.02.007
0142-0615/Ó 2015 Elsevier Ltd. All rights reserved.
 J.C.G. de Siqueira et al. / Electrical Power and Energy Systems 72 (2015) 24–32 25

users to minimize the simulation time without compromising the Ranges of time constants
accuracy of the results. In this case, the main focus of the study
can remain on the analysis of the transient system response, rather This paper develops an algorithm to estimate an adequate Dt
than on a trial and error approach to adjust the time step size and automatically from the circuit. In a first simple approach, we will
the maximum simulation time by running multiple simulations. assume that we have no knowledge of the way the circuit
Moreover, when multiple simulations cases are run through components are connected.
script-based files, and performed considering the variance of net- Here we ‘‘estimate’’ a minimum and maximum value for all the
work or control parameters, this technique seems to be promising R’s combined into a single parallel (minimum) or series branch
and particularly useful to capture the overall phenomena. (maximum); similarly for all the L’s and all the C’s, according to
In principle, it is well known for EMTP-based simulation that (7)–(12). For n R’s, L’s, or C’s we have:
the maximum frequency expected to be present in the simulation
1
(2) should be at least five times less than the Nyquist frequency (1). RMIN ¼ 1
ð7Þ
This provides a reasonable accuracy (error in the order of 3.31%) in R1
þ R12 þ    þ R1n
the modeling of the discretized lumped electrical components 1
LMIN ¼ 1
ð8Þ
(inductor and capacitor) using the trapezoidal rule of integration L1
þ L12 þ    þ L1n
[1,2]. This academic and well known rule (included here for
1
completeness) results in (3) for the calculation of a reasonable time C MIN ¼ 1
ð9Þ
C1
þ C12 þ    þ C1n
step size Dt.
RMAX ¼ R1 þ R2 þ    Rn ð10Þ
1
f Ny ¼ ð1Þ LMAX ¼ L1 þ L2 þ    Ln ð11Þ
2 Dt
1 1 1 1 C MAX ¼ C 1 þ C 2 þ    C n ð12Þ
f max ¼  f Ny ¼  ¼ ð2Þ
5 5 2Dt 10Dt Eqs. (7), (8) and (12) assume that each type of circuit compo-
1 T min
Dt ¼ ¼ ð3Þ nent is all connected in parallel whereas (9)–(11) assume that each
10f max 10 type of circuit component is all connected in series.
In terms of frequency content, Fig. 1 [3,19] illustrates the wide Correspondingly, we can estimate the maximum and minimum
speed of phenomena in electrical power systems. values for inductive or capacitive time constants (i.e., ðT L ÞMAX ;
A common practice for EMTP-based simulation is to always ðT L ÞMIN ; ðT C ÞMAX ; ðT C ÞMIN Þ and for the system natural modes of
adopt a time step size of oscillations (i.e., ðT N ÞMAX ; ðT N ÞMIN Þ.
For a circuit composed by a resistance and an inductance, its
Dt ¼ 50 ls ð4Þ
time constant can be calculated using the expression (13),
which then implies that:
L
TL ¼ ð13Þ
1 1MHz R
f Ny ¼ ¼ ¼ 10 kHz ð5Þ
2  50 ls 100 For a circuit composed of a resistance and a capacitance, the
1 1 1 1 time constant can be calculated using the expression (14),
f max ¼  f Ny ¼  ¼ ¼ 2 kHz ¼ 2  103 Hz ð6Þ
5 5 2Dt 10  50 ls
T C ¼ R:C ð14Þ
If one chooses, for example a Dt of 1 microsecond, even though
very small, it might not be adequate for example, when simulating For a circuit composed of an inductance and a capacitance, its
very high frequency transients in SF6 gas insulated substations. For natural oscillation time period is given by:
the very common Dt of 50 microseconds, only transients up to pffiffiffiffiffiffi
2 kHz will be properly simulated. T N ¼ 2p LC ð15Þ

Fig. 1. Typical frequency ranges for electromagnetic and electromechanical transient phenomena [3,19].
26 J.C.G. de Siqueira et al. / Electrical Power and Energy Systems 72 (2015) 24–32

Thus, it is possible to determine the minimum and maximum 1

where, in general, k1 ¼ 100 , and k2 ¼ 5.
limits for the max and min ‘‘time constants’’, which can be either Observe that, by using the factor k1 ¼ 1=100 of the minimum
inductive or capacitive, (T) in (23), instead of 1/10 as in (3), the transient response at
LMIN t ¼ Dt is closest to the ideal analytical value for t ¼ 0þ , specially
ðT L ÞMIN ¼ ð16Þ
RMAX for simple RLC cases, used to validate the initial version of this
LMAX empirical algorithm.
ðT L ÞMAX ¼ ð17Þ Remember that, the time calculated values in the vector (T),
RMIN
ðT C ÞMIN ¼ RMIN :C MIN ð18Þ except for very simple first order cases, are not the real time
constants. In general, the minimum time calculated value in the
ðT C ÞMAX ¼ RMAX :C MAX ð19Þ
vector (T) represents the lower limit, whereas the maximum time
From (15), it is possible to calculate the minimum and calculated value in the vector (T) represents the upper limit of the
maximum time values related to the system natural oscillations boundaries of all possible time constants of a circuit for all possible
modes, as presented in (20) and (21), combinations, i.e. the vector (T) gives the range in which the real
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
ðT N ÞMIN ¼ 2p LMIN C MIN ð20Þ time constants will inevitably be within.
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi In cases with distributed parameters transmission lines using
ðT N ÞMAX ¼ 2p LMAX C MAX ð21Þ the Bergeron method, then a factor k1 ¼ 1=10 instead of 1/100 is
For transmission lines, the proper characteristic impedance (Z c Þ used in (23), which seems to be good enough. Moreover, in such
and the propagation time (sÞ for each transmission line are the cases with distributed parameters transmission lines using the
parameters that can be considered. However, the characteristic Bergeron method, if the maximum time calculated value in the
impedance (Z c Þ is not used, because it seems that identifying the vector (T) is the sMAX , a factor k2 ¼ 10 instead of 5 might be used
sMIN and sMAX among all transmission lines in a case, and include in (24), which guarantees that, after all wave reflections, steady
them in the vector (T) is good enough in this initial version of state conditions are reached within the simulation window.
the empirical algorithm. Next, test cases with simulation results illustrate the applica-
Independent current and voltage sources, representing external tion of the method. The procedure can be automated in EMTP-
forcing functions of the simulated system, namely generators with based programs. Further developments are possible by integrating
arbitrary f(t), are accounted in this method by their corresponding intelligent logics for power systems and control components,
time period, based on their respective frequencies, using (22): based on their specific time requirements and the necessary nearly
2p 1 optimal Dt and t max .
T1 ¼ ¼ ð22Þ
x f
In this case, trapezoidal integration and the Nyquist theorem
require a value of Dt of at least T 1 =10 [1,2]. When the circuit is
Test cases simulations
excited by voltage or current step-functions (DC sources), from
(22) it results T 1 ! 1, which then are not considered.
Simple RLC circuit
For the case of non-linear circuits that, in EMTP simulations, are
quite common (e.g. non-linear inductances of machines
Fig. 2 presents a RLC single-phase circuit where the minimum
(transformers, generators, CTs), surge arresters, etc.), it is possible
and maximum of the series and parallel combination of the circuit
to extend the empirical method, considering that, in general,
parameters were calculated to obtain the integration step size (dis-
EMTP-based programs represents the nonlinear characteristic
crete time Dt) and the maximum simulation time ðt max Þ.
point-by-point as piecewise linear, from which maximum and
The goal is to get good simulation results for the voltage and
minimum values can be calculated. Then, those values can be used,
current transient in the circuit. The circuit is excited by a voltage
accordingly, in the computation of (7)–(12). For example, in the
step-function (DC source) with zero initial conditions, where
case of a surge arrester, from its voltage versus current character-
R ¼ 100 X; L ¼ 5 H; C ¼ 250 lF, and eðtÞ ¼ 10 V.
istic, the maximum value for its resistance is directly calculated
The minimum and maximum values of the equivalent R, L and C
from the first slope, whereas the minimum value for its resistance
are calculated from Eqs. (7)–(12):
is calculated from the last pair of values (voltage divided by
current). The minimum value for the surge arrester resistance
can then be used in the computation of (7), and the maximum RMIN ¼ 100X; RMAX ¼ 100 X
value for the surge arrester resistance can then be used in the LMIN ¼ 5 H; LMAX ¼ 5 H
computation of (10).
C MIN ¼ 250 lF; C MAX ¼ 250 lF
The method can be extended to power system stability controls
[20] simulated with EMTP-based programs, and circuits with Thus, it is possible to determine the minimum and maximum
time-varying topology (e.g., power electronics applications [21]). limits for the max and min ‘‘time constants’’ from (16)–(19), and
In this case, closing and opening times of switches (i.e. parameters the minimum and maximum time values related to the system
tclose and topen Þ, and their differences have also to be considered, for natural oscillations modes, from (20) and (21). When the circuit
reasons such as being able to assure that the selected time step size
is compatible with minimum switching periods (i.e. maximum
switching frequencies in power electronics-based simulation), vR(t) vL(t)
and also assure that the maximum simulation time will include
the effects of all switching.
i(t) R L
With all the time calculated values sorted in ascending order in
a vector (T), the ‘‘optimum’’ time step size and the maximum
e(t) C vC(t)
simulation time can be estimated as in (23) and (24), respectively: -
Dt ¼ k1  minðTÞ ð23Þ
t max ¼ k2  maxðTÞ ð24Þ
Fig. 2. Test case: RLC single-phase circuit.
 J.C.G. de Siqueira et al. / Electrical Power and Energy Systems 72 (2015) 24–32 27

is only excited by a voltage step-function (DC source), as in this Branch Voltage v (2,1)
simple case, from (22) it results T 1 ! 1, which then is not 10

considered.
Once the time step limits have been calculated, it is now
possible to build a table (Table 1) with these limits organized in
ascending order, as in a vector (T).
Since the calculated time values in the previous analysis are 5
only the limits or boundaries, a practical approach would be to

voltage [V]
choose as the initial optimal time step size ðDtÞ one submultiple
of the smallest calculated time value in Table 1, and as the maxi-
mum simulation time ðt max Þ, a multiple of the biggest calculated
time value in Table 1.
Applying (23) and (24) to Table 1, the results (25) and (26) are 0

obtained.

1 1
Dt ¼  minðTÞ ¼  25 ms ¼ 250 ls ð25Þ
100 100
t max ¼ 5  maxðTÞ ¼ 5  222 ms ¼ 1:11 s ð26Þ -5
0 0.2 0.4 0.6 0.8 1
time [s]
The effectiveness of this criterion can now be tested by simulat-
ing the circuit in an EMTP-based program and analyzing, for Fig. 4. Inductor voltage.
example, the voltage and current in all the elements (please, see
Figs. 3–5). Observe that, by using 1/100 of the minimum TÞ in
(25) (instead of 1/10 as usually calculated with (3)), the transient
response at t ¼ Dt for the voltage across the inductor (Fig. 4) is Branch Voltage v (1,0)
15
closest to the ideal analytical value (vL ðtÞ ¼ 10½V, for t ¼ 0þ ).
This is easily verified in the plotting of Fig. 4, and can also be veri-
voltage [V]

10
fied by zooming over the transient period. Moreover, the initial
estimate for tmax is sufficient to contain the full transient period
within the simulation window. This is relevant for engineering 5

overvoltage and overcurrent analysis, when traveling wave phe-

nomena, resonances, or ferroresonances may show up after some 0
0 0.2 0.4 0.6 0.8 1
elapsed simulation time.
time [s]

Branch Current i (1,0)

Table 1 0.06
Calculated time values for the circuit of Fig. 2.
0.04
current [A]

ðT C ÞMIN 25 ms
ðT C ÞMAX 25 ms 0.02
ðT L ÞMIN 50 ms
ðT L ÞMAX 50 ms 0
ðT N ÞMIN 222 ms
ðT N ÞMAX 222 ms -0.02
0 0.2 0.4 0.6 0.8 1
time [s]

Fig. 5. Capacitor voltage and current.

Branch Voltage v (3,2)
6

4 Since with modern computers, getting an accurate transient

voltage [V]

response is normally more important than minimizing the

2 simulation time, the suggested approach provides a good initial
0
estimation for Dt and tmax .
In any case, users should always have the option of properly
-2 choosing these parameters according to their own experience
0 0.2 0.4 0.6 0.8 1
and criteria, e.g., the specific case under study, background
time [s]
technical experience, and the objectives of the simulation.
Branch Current i (3,2)
0.06
vS(t)
0.04
5 Ls 4 3 Zc=400[ Ω] 2 1
current [A]

vL(t)
0.02 τ =500[µs]

0 Rs i(t) L
e(t) t=0[s] vC(t)
C
-0.02 t=20[ms]
0 0.2 0.4 0.6 0.8 1
time [s]

Fig. 3. Resistor voltage and current. Fig. 6. Transmission line terminated in a LC circuit with switching operation.
28 J.C.G. de Siqueira et al. / Electrical Power and Energy Systems 72 (2015) 24–32

Transmission line terminated in a LC circuit with switching operation x 10

5
L4 Branch Voltage v(2,1)
2

Fig. 6 presents a single-phase transmission line terminated in a

1
LC circuit energized by an AC voltage source behind a parallel RL cir-

voltage [V]
cuit as an equivalent source impedance, where: eðtÞ ¼ 100; 000:
0
cosðxtÞ [V], f = 60 [Hz], Rs = 50 [X], Ls = 400 [mH], t close 1 ¼ 0 [s],
0 6 0 11
L ¼ 2:10 [H/m], C ¼ 1:25:10 [F/m], length = 100 [km], Zc = 400 -1
[X], s = 500 [ls], L = 100 [mH], C = 1 [lF], topen 2 ¼ 20 [ms].
From the circuit parameters, time values were calculated as -2
presented in Table 2. The integration step size (discrete time Dt) 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08

and the maximum simulation time ðt max Þ were proposed by the time [s]
empiric algorithm using (27) and (28).
The computer simulation results for this case are presented in L4 Branch Current i(2,1)
1000
Figs. 7–10, illustrating the electromagnetic transients towards
the steady state period.
500

current [A]
1 1
Dt ¼  minðTÞ ¼  50 ls ¼ 5 ls ð27Þ
10 10 0

t max ¼ 5  maxðTÞ ¼ 5  16:67 ms ¼ 83:33 ms ð28Þ

-500

Lightining, overvoltages, and insulation case study -1000

0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08
time [s]
Insulation coordination studies for practicing engineers can be
made with the EMTP, or with similar programs, for high-voltage Fig. 8. Terminal inductor voltage and current.
substations. The line diagram for a typical study of a 220 kV
substation is shown in Fig. 11 [22]. Lines A, B, C, and D represent
busbar sections, whereas line 5–3 represents the overhead line x 10
5
L4 Branch Voltage v(2,1)
struck by lightning. Line 1 in the diagram is an overhead line which
is so long that no wave returns from the remote end during the
1
duration of t max ¼ 10 ls of this study.
voltage [V]

Table 2
Calculated time values for the circuit of Fig. 6. -1

ðT C ÞMIN 50 ls 1 2 3 4 5 6 7 8 9 10
ðT C ÞMAX 50 ls time [s] -3
x 10
ðsÞMIN 500 ls
ðsÞMAX 500 ls
ðT L ÞMIN 1.60 ms L4 Branch Current i(2,1)
ðT N ÞMIN 1.7772 ms
ðT N ÞMAX 4.4429 ms 500
current [A]

ðT L ÞMAX 10 ms
ðT 1 ÞMIN 16.67 ms 0
ðT 1 ÞMAX 16.67 ms
-500

5
Nodal Voltage (2) 0 0.005 0.01 0.015 0.02 0.025
x 10
3
time [s]

Fig. 9. Terminal inductor voltage and current (zoom).

2

Line 2–4 represents a cable which enters the generator step-up

1
transformer of an underground power plant directly through the
voltage [V]

bushing. The question to be answered in this case is whether a

0 surge arrester located at the cable terminal 2 in the open air
substation can protect the transformer at the end of the cable.
The transformer is represented by its surge capacitance, in par-
-1
allel with a damping resistance approximately equal to the surge
impedance of the transformer winding. The capacitance in 1 repre-
-2
sents a capacitive voltage transformer. The nonlinear characteristic
of the surge arrester and its breakdown voltage are shown in [22].
If the lightning stroke is assumed to hit the phase conductor
-3 directly, 30 m away from tower in 5, and if assumed a piecewise
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08
linear current impulse with imax ¼ 8 kA at 1 ls, and i ¼ 4 kA at
time [s]
50 ls, it is obtained the overvoltages and currents shown in
Fig. 7. Nodal voltage at node 2. Figs. 12–14.
 J.C.G. de Siqueira et al. / Electrical Power and Energy Systems 72 (2015) 24–32 29

5
x 10 C3 Branch Voltage v(1,0)
4

2
voltage [V]

-2

-4
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08
time [s]

C3 Branch Current i(1,0)

1000

500
current [A]

-500

-1000 Fig. 12. Voltage in [kV] in location 1 and 2 [22].

0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08
time [s]

Fig. 10. Terminal capacitor voltage and current.

Tower with volt-me

Busbar secons A, B, C, D with
i characterisc of insulator string
Z = 370 Ω , c = 300 m/μs
A
l = 30m
t Z = 370 Ω , c = 300 m/μs
3
l = 30 m 5 l = 300 m B
Z = 30 Ω , c = 150 m/μs 4
2
l = 150m
C surge
Z = 370 Ω , c = 300 m/μs 1 arrester 5000 Ω
l=∞ 3Ω 1000 pF
4400 pF D

Fig. 11. Line diagram of 220 kV substation [22].

Fig. 13. Lightining stroke current in [kA] and discharge current in the surge arrester
at location 2 [22].
Fig. 12 shows the voltage in location 1 and 2, while Fig. 13
shows the lightning stroke current wave shape, and the discharge
current in the surge arrester.
Fig. 14 shows the voltage at the transformer terminal 4, and at
node 3. It can be seen that the overvoltage at the transformer
terminal (location 4) is somewhat higher ð722:4 kVÞ than at the
surge arrester (location 2) ð640:5 kVÞ, because of the cable connec-
tion in between.
The volt-time characteristic of the 220 kV insulator string at the
tower location 5 can also to be considered as explained in [22]. For
imax ¼ 8 kA there is no flashover across the insulator. For
imax ¼ 10 kA the overvoltage in 5 would have intersected the
volt-time characteristic, and a flashover to the tower would have
occurred.
All these results were simulated using a Dt ¼ 0:1 ls (the
minimum propagation time), and tmax ¼ 10 ls.
However, if for simplicity we ignore the surge arrester and the
volt-time characteristic that determines possible flashover to the
tower, and then apply the suggested algorithm for the calculation
of the optimal time step size and maximum simulation time, we
will find the range of calculated time values of Table 3, in
Fig. 14. Voltages in [kV] at the transformer terminal 4, and at node 3 [22].
ascending order.
30 J.C.G. de Siqueira et al. / Electrical Power and Energy Systems 72 (2015) 24–32

Table 3 Physical meaning arises when also plotting all the curves.
Calculated time values for the circuit of Fig. 11. Nevertheless, a good engineering approach is to pay attention before
ðsÞMIN 0.1 ls in the modeling of the system, and all the simplifications made in
ðsÞMAX 1 ls the equivalent circuit for the phenomena under investigation.
ðT C ÞMIN 5 ls
ðT C ÞMAX 22 ls
Calculating continuous-time eigenvalues from the digital EMTP
formulation

By applying (23) and (24) to Table 3, the algorithm automati- Step-by-step eigenvalue analysis with EMTP discrete-time
cally assumes for this case that k1 ¼ 1=10, and k2 ¼ 5, resulting solutions is presented in [15,16]. Essentially, a discrete-time
in (29) and (30). state-space matrix is first calculated using (31), from which
1 1 discrete-time eigenvalues ðzi Þ and their respective eigenvectors
Dt ¼  minðTÞ ¼  0:1 ls ¼ 0:01 ls ð29Þ ðv i Þ are calculated.
10 10
t max ¼ 5  maxðTÞ ¼ 5  22 ls ¼ 110 ls ð30Þ Based on (32)–(35) continuous-time eigenvalues ðki Þ are now
calculated using (36). As long as (31) does not change (e.g.,
A holistic and practical engineering comprehension of this case switching events cause changes in the conductance matrix ½GAA )
is possible by running the EMTP-based simulation considering the the continuous-time eigenvalues remain the same as in (36).
whole transient period, i.e. with t max ¼ 200 ls. Fig. 15 presents the Moreover, the continuous-time eigenvalues do not depend on
voltages at the stroke location, and at the transformer terminal the simulation time step size [6,7]. The equations relating the
(location 4), when not including the surge arrester, and running network topology and the EMTP conductance matrix [GAA] are
the simulation for this case with tmax ¼ 200 ls. given by:
When the maximum instantaneous values or peak values of h i
Ad ¼ ½K ^ þ ½g^½L½GAA 1 ½Lt ð31Þ
each nodal voltage are detected, as presented in Table 4, this
reveals that the maximum reached overvoltages (up to 1480 kV, h i  2   
  
2

  1
i.e. 1:48 MV in some nodes) occur as expected just at the beginning Ad ¼ ½I þ Ac  ½I  Ac ð32Þ
Dt Dt
of the lightning transient wave propagation. At the transformer h i
d
terminal (location 4), without surge arrester the overvoltage is A :v i ¼ zi :v i ð33Þ
 c
1,168.61 [kV], whereas for the case with surge arrester, the A :v i ¼ ki :v i ð34Þ
overvoltage would be limited to 722.404 [kV], based on the EMTP  2 
Dt
þ ki
modeling for this case. 
zi ¼ 2  ð35Þ
Dt
 ki
2
Dt
½zi  1
ki ¼ ð36Þ
½zi þ 1
where
h i
Ad = discrete-time state-space matrix;
h i
b = branch parameters according to branch nature, i.e.,
K
b kk ¼ þ1 for inductive branches, K
K b kk ¼ 1 for capacitive branches,
b b
K kk ¼ 0 for resistive branches, K jk ¼ K b kj ¼ 0;
½g^ = branch parameters for history current source calculation,
defined according to branch nature, i.e., g ^kk ¼ 2 D2Lt for inductive
2C
branches, g^kk ¼ 2 Dt for capacitive branches, g^kk ¼ 0 for resistive
branches, g ^jk ¼ g^kj ¼ 0;
½L = Graph connectivity or branch-node incidence matrix,
where for any branch connected between nodes k, and
m; v km ¼ þ1v k  1v m . Then the coefficients +1 and 1 are then
used to build ½L such that ½v km  ¼ ½L  ½v , where ½v km  is the vector
of branch voltages, and ½v  is the vector of nodal voltages;
Fig. 15. Voltages in [kV] at the stroke location, and at the transformer terminal 4,
½GAA  = conductance matrix after partition to exclude nodes with
when not including the surge arrester, during the whole transient period.
voltage sources connected to;
½Lt = transpose of ½L;
Table 4 Dt = time step size;
Maximum nodal voltage magnitudes at the circuit of Fig. 11. ½I = unit matrix, where all elements are zero, except the
Node Maximum voltage with Maximum voltage without surge diagonals, which are equal to 1;
 c
name surge arrester arrester (MV) A = continuous-time state-space matrix;
STROKE 1.48000 MV 1.48000 v i = eigenvector associated with the eigenvalue with index i;
NODE 5 1.48000 MV 1.48000 zi = discrete-time eigenvalue of index i;
NODE 3 638.333 kV 1.16509
ki = continuous-time eigenvalue of index i.
NODE 3 641.416 kV 1.16549
TOP The relevance of this formulation is that the dynamics of the
NODE 2 640.456 kV 1.16651 system can be determined for every event, and if desired, at every
NODE 4 722.404 kV 1.16861 instant of the simulation, as the topology changes by switches
NODE 1 648.637 kV 1.16329
operations, or as we change regions in non-linear elements.
NODE 1 651.179 kV 1.16370
LOW
Also, fast and slow subsystem areas can be identified for use
by multisystem programs, like MATE – Multi-Area Thévenin
 J.C.G. de Siqueira et al. / Electrical Power and Energy Systems 72 (2015) 24–32 31

Equivalent [23], Multi-level MATE [24], and Latency exploitation the network topology is also suggested. This algorithm does not
[25]. These techniques allow the use of different time-step sizes make the simplifications of the empirical method and has the
in each separate subsystem, with further resynchronization [26] potential of providing very accurate pre-determinations of the
at specific times, thus optimizing the overall solution in EMTP- optimum time step and solution times.
based off-line or real-time digital simulations [27]. User-independent determination of the simulation parameters is
The continuous-time eigenvalue calculation with (36) can also of particular importance for new users of EMTP-based programs.
be used to find the exact values for optimum time step size and Existing programs require the user to specify these parameters,
maximum simulation time. The eigenvalues calculated with this which can require extensive experience in performing transients
technique are exact. simulations. A future challenge will be the extension of these
Considering the branches associated with inductive and techniques to adapt the time step size dynamically in EMTP-based
capacitive components and with transmission lines, for each com- programs as the topology and parameters of the system change
plex continuous-time eigenvalue in (37), it is possible to calculate during the simulation.
the corresponding exact time constant as in (38), and the damped
oscillation period as in (39). Acknowledgement
ki ¼ ri  jxdi ð37Þ
The authors gratefully acknowledge the financial support from
1 CAPES, CNPq, FAPEMIG (Brazil) and NSERC (Canada). An earlier
si ¼  ð38Þ
ri version of this paper was published at PSCC 2014 – 18th Power
2p Systems Computation Conference, August 18–22, 2014, Wroclaw,
Tndi ¼ ð39Þ
xdi Poland, organized by Power Systems Computation Conference
With these time values organized in ascending order in a vector and Wroclaw University of Technology.
(T), we can now derive the exact optimum time step size (23) and
the maximum simulation time (24) as in the simplified procedure References
first presented.
By applying (31)–(39) to the case of Fig. 2 results: [1] Dommel HW. Digital computer solution of electromagnetic transients in
single- and multiphase networks. IEEE Trans Power Apparatus Syst 1969;PAS-
k1 ¼ 10 þ j26:4575 ð40Þ 88:388–99.
[2] Dommel HW. Electromagnetic transients program reference manual.
k2 ¼ 10  j26:4575 ð41Þ Vancouver: Department of Electrical Engineering, The University of British
Columbia, Canada; 1996.
1
s1 ¼ s2 ¼  ¼ 100 ms ð42Þ [3] Greenwood A. Electrical transients in power systems. John Wiley & Sons; 1991.
10 [4] Microtran web site. <http://www.microtran.com/>.
2p [5] Alternative Transients Program (ATP) web site. <http://www.emtp.org/>.
Tnd1 ¼ Tnd2 ¼ ¼ 237:48 ms ð43Þ [6] EMTP-RV web site. <http://emtp.com/>.
26:2475 [7] PSCAD web site. <https://hvdc.ca/pscad/>.
1 1 [8] RTDS web site. <http://www.rtds.com/index/index.html>.
Dt ¼  minðTÞ ¼  100 ms ¼ 1 ms ð44Þ
100 100 [9] OPAL-RT Technologies web site. <http://www.opal-rt.com/>.
[10] Cadence. Cadence PSpice A/D and advanced analysis. <http://www.
t max ¼ 5  maxðTÞ ¼ 5  237:48 ms ¼ 1:19 s ð45Þ cadence.com/products/orcad/pspice_simulation/pages/default.aspx>.
In this simple case, when comparing the result obtained by (25) [11] Intusoft. Solving SPICE convergence problems. <http://www.intusoft.com/
articles/converg.pdf>.
with the one obtained by (44) and also the result obtained by (26) [12] Barros FJ. Comparing synchronous and asynchronous variable step size explicit
with the one obtained by (45), one has to assess the additional ODE solvers: a simulation study. 21st International workshop on principles of
computational effort of performing the most accurate analysis with advanced and distributed simulation (PADS ’07); 2007. p. 32–7. http://dx.doi.
org/10.1109/PADS.2007.17.
that of the simplified analysis of Section ‘Ranges of time constants’.
[13] Martí JR, Lin J. Suppression of numerical oscillations in the EMTP. IEEE Trans
More complex cases have to be developed to further assess the Power Syst 1989;4(2):739–47.
validity of the proposed algorithms. [14] Lin J, Martí JR. Implementation of CDA procedure in the EMTP. IEEE Trans
Power Syst 1990;5(2):394–402.
The most important issue, however, is the emphasis on the
[15] Hollman JA. Step by step analysis with EMTP discrete time solutions. Ph.D.
learning process aiming to optimize the analysis of the transient thesis, The University of British Columbia, Vancouver; 2006.
simulation. If optimal time step size is not used, relevant phenom- [16] Hollman JA, Martí JR. Step-by-step eigenvalue analysis with EMTP discrete-
ena may not be represented in the simulation. time solutions. IEEE Trans Power Syst 2010;25(3):1220–31. ISSN: 0885-8950.
[17] de Siqueira JCG. Introduction to numerical analysis of linear electrical circuits
Moreover, if optimal maximum simulation time is not used, with time invariant lumped parameters. Federal School of Engineering of
overvoltages and or overcurrents which may arise along time Itajuba, Itajuba, MG; 1992. p. 99.
due to traveling wave and or resonance phenomena may not be [18] de Siqueira JCG, Bonatto BD, Martí JR, Hollman JA, Dommel HW. Optimum
time step size and maximum simulation time in EMTP-based programs. In:
properly considered in the engineering insulation study or PSCC 2014 – Power Systems Computation Conference, organized by Power
overcurrent protection coordination, eventually resulting in time, Systems Computation Conference and Wroclaw University of Technology,
equipment or system failures with technical, economic, environ- Wroclaw, Poland, August 18–22; 2014.
[19] Henschel S. Analysis of electromagnetic and electromechanical power system
mental or human losses. transients with dynamic phasors. Ph.D. thesis, University of British Columbia,
Vancouver, BC, Canada; 1999.
[20] Kundur P. Power system stability and control. New York: McGraw-Hill; 1994.
Conclusions [21] Bonatto BD, Dommel HW. EMTP modelling of control and power electronic
devices – electromagnetic transients programs helping the analysis of the
power interaction on either the load or the network side. LAP – Lambert
This paper has provided algorithms to automatically determine Academic Publishing AG & Co; 2010. 175p. ISBN: 978-3-8383-2790-7.
the time step size and the simulation length in EMTP-type of [22] Martí JR, Dommel HW. Line models for lightning studies. Transactions
programs. A simple empiric method was proposed to calculate an Engineering & Operation Division, Canadian Electrical Association, vol. 28;
1989.
estimate for the optimal time step size ðDtÞ and for the maximum
[23] Martí JR, Linares LR, Hollman JA, Moreira FA, OVNI: integrated software/
required simulation time ðt max Þ. This would help the user by e.g. hardware solution for real-time simulation of large power systems. In:
minimizing running time or providing a warning if an inappropri- Proceedings of the 14th power systems computation conference (PSCC02),
ate time-step has been set. Seville, Spain, June 24–28; 2002. <http://www.pscc02.org>.
[24] Armstrong ML, Marti JR, Linares LR, Kundur P. Multilevel MATE for efficient
A more exact algorithm, based on the determination of the cir- simultaneous solution of control systems and nonlinearities in the OVNI
cuit continuous time eigenvalues from the circuit parameters and simulator. IEEE Trans Power Syst 2006;21(3):1250–9.
32 J.C.G. de Siqueira et al. / Electrical Power and Energy Systems 72 (2015) 24–32

[25] Moreira FA, Martí JR. Latency techniques for time-domain power Computation Conference (PSCC02), Seville, Spain, June 24–28; 2002. <http://
system transients simulation. IEEE Trans Power Syst 2005;20 www.pscc02.org>.
(1):246–53. [27] Matar M, Iravani R. FPGA implementation of the power electronic converter
[26] Linares LR, Marti JR. A resynchronization algorithm for topological changes in model for real-time simulation of electromagnetic transients. IEEE Trans
real time fast transients simulation. In: Proceedings of the 14th Power Systems Power Delivery 2010;25(2):852–60.