XVth International Symposium on High Voltage Engineering

University of Ljubljana, Elektroinštitut Milan Vidmar, Ljubljana, Slovenia, August 27-31, 2007 T10-539.pdf

International comparison of software for calculation of lightning impulse parameters based

on a new processing algorithm
J. Hällström1* , A. Bergman2, Ding 3, F. Garnacho 4, R. Gobbo 5, T. Kato 6, Y. Li7, A. Nilsson2, G. Pesavento 5, S. Sato 8, A. X. Yu3
1
Helsinki University of Technology (TKK), Espoo, Finland
2
Technical Research Institute of Sweden (SP), Borås, Sweden
3
Beijing Hua Tian Mechanical-Electrical Institute (BHT), Beijing, China
4
Electrotechnical Official Central Laboratory (LCOE), Madrid, Spain
5
University of Padova, Padova, Italy
6
Doshisha University, Kyoto, Japan
7
National Measurement Institute (NMI), West Lindfield, Australia
8
Utsunomiya University, Tokyo, Japan
*Email: [email protected]

Abstract: A new algorithm has been proposed to Schering Institute (Germany), LCOE (Spain) and NGC
calculate the parameters of full lightning voltage (United Kingdom) in coordination with the CIGRE
impulses. The new algorithm enables the application of WG 33.03, in order to remove the present ambiguities in
the test voltage factor (also referred to as k-factor in the evaluation of the test voltage when oscillations
some literatures) for calculation of the equivalent test and/or an overshoot appear superimposed in a lightning
voltage of impulses with superimposed impulse.
oscillations/overshoots. The new algorithm at the same Within the project breakdown tests were performed
time provides a robust procedure for obtaining time on different materials (air, oil, XLPE, and SF 6) to study
parameters of the impulses from not only smooth the influence of oscillations and overshoot
waveforms but also waveforms with varying degrees of superimposed on a full (close to double exponential)
distortions in the front part of the impulses. These lightning impulse. The results showed that:
distortions include oscillations on the impulse front and • All investigated materials showed similar behaviour.
overshoots in the peak region. A critical part of the new • The results obtained match with earlier results [1].
algorithm is a 4-parameter fitting procedure to obtain • The 50% disruptive discharge voltage, U50, of the
the base curve, which is used for calculation of the test material under investigation for a full lightning
voltage curve. Another important part of the algorithm impulse without oscillations, Ut, lies between the
is applying a filtering procedure in the calculation of the extreme value of the lightning impulse with
test voltage curve. superimposed oscillations, Ue, and the peak value of
the mean curve of the impulse, Ump. How far Ut is
The new algorithm was tested in different laboratories placed from Ue and Ump depends on the oscillation
using different programming languages and different frequency. The higher the frequency of the
techniques for realising the fitting and filtering routines. oscillation or overshoot the smaller its influence on
The paper reports the results obtained from the the breakdown voltage. This leads to the empirical
participating laboratories using the proposed algorithm. equation with which the test voltage of a material
The results obtained by the participating laboratories subjected to a lightning impulse with oscillations or
using existing software based on the requirement of IEC overshoot can be calculated by using a test voltage
60060-1: 1989 were also compared. It is anticipated that factor k (or referred to as the k-factor in some early
the results can serve as a part of the basis for a new papers) which varies from 1 to 0 depending on the
procedure for determination of lightning impulse oscillation frequency:
parameters in the revised IEC 60060-1. U t = U mp + kβ , (1)
where β = U e − U mp
1 INTRODUCTION .

The result of the study was that in theory it should

The current version of IEC 60060-1 (1989) and be possible to find a better evaluation method that gives
IEEE Std 4 define that oscillations on lightning more coherent results, which are more in agreement
impulses with frequencies above 500 kHz should be with the actual breakdown behaviour of the materials or
disregarded, and that a mean curve should be drawn for devices under test.
the parameter evaluation. This definition has not been
very practical, and it has led to ambiguity in the
parameter evaluation of commonly recorded lightning
impulses.
At the end of the 90’s, a European Project within the
STM Program was developed by KEMA (Netherlands),

1
2 PROPOSED NEW PARAMETER 3.2. Filtering
EVALUATION ALGORITHM Either FIR or IIR type digital filter can be used. Both
The procedure is an implementation of Equation 1, approaches described here lead to a linear-phase
and it is used for computer aided calculation of digitally response.
recorded impulses. The procedure is, in short, to obtain 2.3.1. FIR
a test voltage curve from which the impulse parameters Typical number of coefficients needed to implement
are calculated. the k-factor using linear phase FIR approach range from
The proposed LI parameter definitions for the new 500 to 4000. Again, FIR filter design functions are
IEC 60060-1 are listed in Appendix 1, and the algorithm available in many programming environments [8].
proposed to be added as a normative annex is repeated 2.3.2. IIR
in Appendix 2. To obtain linear phase response, a dual pass (zero-
A number of alternative evaluation algorithms were phase) filtering approach has been proposed [ 9]. In this
studied before starting this round-robin, and the most approach the attenuation of the filter is only half of what
promising one was selected for this study. Discussion is needed, but the data is passed through the filter twice,
on the alternative algorithms can be found on another first forward and then in reverse order.
paper [2]. It appears that in this case an IIR filter can be
A parallel round-robin has been carried out to study implemented by using only three coefficients:
the applicability of the algorithm of equation 1 for y (i ) = b0 x (i) + b1 x (i − 1) + a1 y (i − 1) , (3)
manual evaluation of impulse parameters [3]. where x (i ) is the ith input sample for the filter, and
y (i ) is its ith output sample. Here b0, b1 and a1 are the
3 GUIDANCE FOR PROGRAMMING filter coefficients.
This chapter summarizes our experience on The filter coefficients can be derived by a relatively
implementing residual filtering algorithms [4], [5]. Here simple forward calculation:
we concentrate on two steps of the algorithm that are πT  (4)
x = tan S  ,
not straightforward to implement: double exponential  a 
curve fitting and design of the k-factor filter. x (5)
b0 = b1 = , and
3.1. Double exponential fitting 1+ x
1− x (6)
The double exponential function should include four a1 = ,
free parameters (A, B, C and D): 1+ x
where a is the -3 dB point of the K-factor filter, and TS
 − (t − D ) −
(t − D ) 
is the sampling interval used for recording the signal.
u (t ) = A e B − e C  (2)
  For example, 10 ns sampling interval leads
 .
to a1 = − 0.9585113 and b0 = b1 = 0 .02074434 .
The initial guess for the fitting algorithm could be
The filtering is then performed twice (once in
for example:
forward and once in backward direction) using the
A: The extreme value of the curve
following difference equation:
B: 70µs
y (i ) = 0.02074434(x (i) + x (i − 1)) + 0.9585113 y (i − 1) , (7)
C: 0.4µs
D: True or virtual origin of the curve where y is the output data vector and x the input data
Levenberg-Marquardt algorithm [6] and its vector. In order to avoid numerical problems often
derivatives have been used successfully for finding the typical for IIR filters, large enough number (in this case
best estimate in least squares sense. This algorithm is preferably ≥ 7) of decimals has to be used for filter
readily available in many programming environments coefficient.
[7].
Some software packages do not use partial 4 RESULTS
derivatives during the optimization process. In the case
The evaluation method was tested on 52 impulse
of typical impulse voltage data, where the x- and y-
waveforms. Their sources are summarised in Table 1.
values are of different magnitudes, time values in the
The waveforms vary a lot, and many of them are not
order of 10 -6 and voltage values of 10 3, omitting the
considered as standard lightning impulses by the current
derivatives leads to very poor convergence, or failure to
IEC 60060-1.
converge at all. In this case normalizing the data (i.e.
scaling it so that both voltage and time scales span
approximately from 0 to 1) has been necessary.

2
Tab. 1: Origins of the 52 curves analysed by each round robin 5 CONCLUSIONS
participant.
A new method for evaluation of impulse voltage
Case Origin
parameters has been proposed for the next revision of
1, 3, 4, 6, 8 and 9 From TDG (IEC 61083-2:1996).
A1 to A17 Created analytically. IEC 60060-1. The method has been tested by evaluating
B1 to B4 From transformer test s. the impulse parameters of 52 impulses both by the new
C1 to C4 Measured smooth waveforms . proposed method, and by the method of the current
D1 to D6 Measured waveforms with front oscillations. version of the standard.
E1 to E15 Misc. non-standard waveforms .
Each participant of this comparison created their
own software according to the algorithm proposed for
Each laboratory created their own software for the next revision of IEC 60060-1. These software
evaluation of the parameters. All participants provided packages (N=7) agree with very low standard
results calculated by the new method, and in addition deviations, when the test voltage value (Ut), front time
four participants delivered parameters evaluated by a
(T1), time to half value (T2) or relative overshoot (β’) is
method fulfilling the requirements by the current
measured. The agreement is good even for the distorted
version of IEC 60060-1.
waveforms e.g. from transformer testing.
The difference between the old method (as described
The results obtained by the old method show,
in the current version of IEC 60060-1) and the proposed
especially for front time evaluation, a large scatter: it
method (as described in the Appendices of this paper)
has to be noticed, however, that several waveforms
are shown in Figs. 1 to 4. The relative overshoot values
being proposed had parameters outside the range which
evaluated for each waveform according to the proposed
could be considered acceptable: in fact, the largest
new method are shown in Fig. 5.
systematic differences between the old and proposed
The results are scaled so that the zero line of Figs. 1
new method were found while evaluating the analytical
to 3 (Ut/Up, T1, and T2) is the mean value obtained by
impulses with very large (typically 10%, and up to 30%)
the proposed new method. In Fig. 4 (β/β’) the zero line overshoot.
means no overshoot. The error bars on the zero line of It seems that the scatter of the old method is larger
each figure show one standard deviation of the results than the systematic changes the introduction of the new
(N=7) obtained by the proposed new method. The dots method would cause to the parameter values.
on these graphs show the results obtained by different
software packages (N=4) fulfilling the requirements of
the current version of IEC 60060-1.
6%

U t, U p
4%
Difference of the old evaluation from the new method

2%

0%

-2 %

-4 %

-6 %
4 - LIFA
6 - LIFO

A1
A2
A3
A4
A5
A6
A7
A8
A9
A10
A11
A12
A13
A14
A15
A16
A17
B1
B2
B3
B4

E1
E2
E3
E4
E5
E6
E7
E8
E9
E10
E11
E12
E13
E14
E15
1 - LI
3 - LISL

8 - LILO
9 - LISO

C1
C2
C3
C4
D1
D2
D3
D4
D5
D6

Curve

Fig. 1 Peak value calculated according to current IEC 60060-1 (Up) compared with the test voltage value ( Ut) of the new proposal. The error bars
around the zero line show one standard deviation of the test voltage value evaluations (N=7). Values obtained by software packages (N=4) according
to the current IEC 60060-1 method is shown by dots.

3
 10 %

T1
5%
Difference of the old evaluation from the new method

0%

-5 %

-10 %

-15 %

-20 %
1 - LI

C1
C2
C3
C4
D1
D2
D3
D4
D5
D6
A10
A11
A12
A13
A14
A15
A16
A17

E10
E11
E12
E13
E14
E15
9 - LISO
8 - LILO
3 - LISL

A1
A2
A3
A4
A5
A6
A7
A8
A9

B1
B2
B3
B4

E1
E2
E3
E4
E5
E6
E7
E8
E9
6 - LIFO
4 - LIFA

Curve

Fig. 2 Front time (T1) calculated according to current IEC 60060-1 compared with the value according to the new proposal. The error bars around the
zero line show one standard deviation (N=7) of the front time evaluations with the proposed new method. Values obtained by software packages
(N=4) according to the current IEC 60060-1 method is shown by dots.

8%

T2
6%
Difference of the old evaluation from the new method

4%

2%

0%

-2 %

-4 %

-6 %

-8 %

-10 %
1 - LI

6 - LIFO
4 - LIFA

C1
C2
C3
C4
D1
D2
D3
D4
D5
D6
A10
A11
A12
A13
A14
A15
A16
A17

E10
E11
E12
E13
E14
E15
9 - LISO
8 - LILO
3 - LISL

A1
A2
A3
A4
A5
A6
A7
A8
A9

B1
B2
B3
B4

E1
E2
E3
E4
E5
E6
E7
E8
E9

Curve

Fig. 3 Time to half value ( T2) calculated according to current IEC 60060-1 compared with the value according to the new proposal. The error bars
around the zero line show the standard deviation (N=7) of the front time evaluations with the proposed new method. Values obtained by software
packages (N=4) according to the current IEC 60060-1 method is shown by dots.

4
 10

β, β '
5
Difference of the old evaluation from the new method [%]

-5

-10

-15

-20

-25
1 - LI

6 - LIFO
4 - LIFA

C1
C2
C3
C4
D1
D2
D3
D4
D5
D6
A10
A11
A12
A13
A14
A15
A16
A17

E10
E11
E12
E13
E14
E15
9 - LISO
8 - LILO
3 - LISL

A1
A2
A3
A4
A5
A6
A7
A8
A9

B1
B2
B3
B4

E1
E2
E3
E4
E5
E6
E7
E8
E9
Curve

Fig. 4 Overshoot ( β) calculated according to current IEC 60060-1 compared with the relative overshoot (β’) according to the new proposal. The error
bars around the zero line show the standard deviation (N=7) of the relative overshoot values. Values obtaine d by the software packages (N=4)
according to the current IEC 60060-1 method is shown by dots.

25 %

β'
20 %

15 %
Relative overshoot magnitude

10 %

5%

0%

-5 %
4 - LIFA

9 - LISO
1 - LI

C1
C2
C3
C4
D1
D2
D3
D4
D5
D6
8 - LILO
3 - LISL

A10
A11
A12
A13
A14
A15
A16
A17

E10
E11
E12
E13
E14
E15
6 - LIFO

A1
A2
A3
A4
A5
A6
A7
A8
A9

B1
B2
B3
B4

E1
E2
E3
E4
E5
E6
E7
E8
E9

Fig. 5 Relative overshoot magnitude evaluated using the proposed new method. Results of 7 evaluations are shown for each waveforms.

5
6 APPENDIX 1, DEFINITIONS l) apply the digital filter to the residual curve R(t) to
obtain the filtered waveform Rf (t);
Definitions in IEC60060-1 draft used for this work:
m) add the filtered residual curve Rf (t) to the base
- test voltage value, Ut curve Um(t) to obtain the test voltage curve , Ut(t);
peak value of the test voltage curve
n) calculate the value of the test voltage Ut and time
- overshoot magnitude, β parameters from the test voltage curve;
difference in peak values between the recorded
o) find the peak value Ump of the base curve Um(t);
curve and the base curve
p) calculate the relative overshoot magnitude,
- front time, T1
of a lightning impulse voltage is a virtual parameter U e − U mp
β ′ = 100 ⋅ %;
defined as 1 0, 6 times the interval T between the Ue
instants when the impulse is 30 % and 90 % of the q) display the recorded curve U(t) and the test voltage
peak value on the test voltage curve. curve Ut(t).
- time to half-value, T2 r) report the value of the test voltage, front time, time
of a lightning impulse voltage is a virtual parameter to half value, and relative overshoot magnitude, β′.
defined as the time interval between the virtual
origin, O1, and the instant when the test voltage 8 REFERENCES
curve has decreased to half the peak value
[1] J.E. Clem, J.R. Meador, W.J. Rudge, and A.H.
7 APPENDIX 2, EVALUATION Powell: “Proposed basic impulse insulation levels
PROCEDURE for high-voltage systems”, AIEE transactions, pp
953-963, Vol. 69, 1950.
The guidance given to the comparison participants
[2] A. Nilsson, A. Mannikoff, J. Hällström and Y. Li:
for creating their evaluation procedure:
“New Procedures for Determination of Parameters
a) remove voltage offset from the recorded curve of Lightning Impulse Voltage Waveforms”,
U(t), and use that curve for the remaining steps; ISH´07.
b) find the extreme value Ue of the recorded curve U(t); [3] S.M. Berlijn, M.Gamlin, T.Steiner, E. Pultrum,
c) find the base level of the recorded curve by J. Rickmann, W. Lick: “Manual Evaluation of
calculating the mean of the voltage values from the Lightning Impulses according to the new IEC
flat part in the beginning of the record; 60060-1 (measurement uncertainty and method)”,
d) find the last sample on front having a voltage value ISH’07.
less than 0,2 times the extreme value Ue; [4] J. Hällström, S. Berlijn, M. Gamlin, F. Garnacho,
e) discard data up to and including that sample; E. Gockenbach, T. Kato, Y. Li and J. Rungis,
”Applicability of different implementations of K-
f) find the last sample on the tail having a voltage factor filtering schemes for the revision of
value larger than 0,4 times the extreme value Ue;
IEC60060-1 and -2” in Proceedings of the XIVth
g) discard data after that sample; International Symposium on High Voltage
h) fit the following double exponential function to the Engineering, 2005, paper B-32.
remaining data: [5] Y Li and J Rungis, “Evaluation of Parameters of
 − (t − D ) −
(t − D )  Lightning Impulse with Overshoot”, in
u d (t ) = A e B − e C  . Proceedings of 13 th International Symposium on
 
  High-Voltage Engineering, Delft, The
Here t is time, ud(t) is the double exponential voltage Netherlands, 25 – 29 August 2003.
function, and A, B, C and D are the parameters to be [6] Marquardt, D., "An Algorithm for Least-Squares
found by fitting; Estimation of Nonlinear Parameters," SIAM
Journal Applied Math., Vol. 11, pp. 431-441,
i) construct the base curve Um(t) of the waveform, by
1963.
using the base level of the recorded curve for sample
[7] e.g. Matlab (lsqcurvefit), Octave (leasqr), Labview
points up to time D (as defined in step h)) and values (Nonlinear Curve Fit) and LabWindows
of ud(t) for sample points from time D up to the
(NonLinearFit).
instant of the last sample defined in step f). [8] e.g. Matlab (fir2), Octave (fir2).
j) subtract the base curve U m(t) from the recorded [9] Lewin, Tran, Swaffield and Hällström: “Zero
curve U (t) to obtain the residual curve Phase Filtering for Lightning Impulse Evaluation:
R(t) = U(t) – Um(t); A K -factor Filter for the Revision of IEC60060-1
k) create a digital filter with its transfer function H(f) and -2”, accepted for publication in IEEE Power
equal to that defined by the test voltage factor Delivery, see also CIGRE D1.33 WGTF1 (Larnaca
function k(f); 2006) IWD27.