Application of Genetic Algorithm to Estimation of Function Parameters in Lightning Currents Approximations

Genetic algorithm (GA) is applied for the estimation of two-peaked analytically extended function (2P-AEF) parameters in this paper. 2P-AEF is used for approximation of measured and typical lightning discharge currents. Lightning discharge channel is often modeled as thin-wire vertical antenna at perfectly conducting ground. Engineering lightning stroke models assume that the current along that channel is related to the channel-base current which may be measured at the instrumented tall towers and in triggered lightning experiments. Mathematical modeling of lightning currents is important in verification of lightning strokes models based on simultaneously measured electromagnetic fields at various distances, so as in lightning protection studies, computation of lightning induced effects and simulation of overvoltages in power systems. Typical lightning discharge currents of the first positive, first negative, and subsequent negative strokes are defined by IEC 62305 Standard based on comprehensive measurements. Parameters of 2P-AEF’s approximation of the typical negative first stroke current are determined by GA and compared to approximations obtained by other functions. Measured currents at Monte San Salvatore in Switzerland, at Morro de Cachimbo Station in Brazil, and in rocket-triggered lightning experiments at Camp Blanding in Florida are approximated by 2P-AEFs, and good agreement with experimentally measured waveshapes is obtained.


Introduction
Besides cloud-to-ground lightning discharges, intracloud and discharges between thunderclouds, there are also great transient luminous events above troposphere, such as blue jets, gigantic jets, red sprites, elves, and halos.A small percentage of powerful thunderstorms, mostly with positive cloud-to-ground flashes and lightning discharge currents of the order of hundreds of kA, trigger such phenomena at tens of kilometers above thunderclouds.However, the majority of cloud-to-ground lightning discharges (about 90%) transfer negative charges from thunderclouds to the ground.Typical waveshapes of the first positive, first negative, and subsequent negative strokes currents are obtained from comprehensive measurements of Berger et al. [1,2] and Lin et al. [3] and measurements at tall instrumented towers [4] and in rockettriggered lightning experiments [5,6].These currents, so as some estimated parameters of lightning discharges, are often used in modeling of lightning strokes, in research and lightning protection studies [2,7], and in the IEC 62305 Standard [8].Engineering models [9] assume that the current pulse propagates with distortion and attenuation along the channel between the ground and the center of charges in thunderclouds.Spatial and temporal current distribution (or charge-density distribution) along the channel radiates lightning electromagnetic field over the conducting ground.Based on lightning strokes parameters (measured channel-base currents and their derivatives, estimated speed of propagating front, and return stroke speed evaluated from the channel luminosity), model-predicted electromagnetic field results are compared to the measured data at various distances from the channel in order to validate such models.Electromagnetic models assume that the lightning channel is a lossy thin-wire monopole antenna at perfectly conducting ground.Current distribution along the channel is determined by solving Maxwell's equations, for the satisfied boundary conditions, and afterwards, electric and magnetic fields are computed and compared to the measured ones.The influence of finite ground conductivity may be also taken into account in these models.Approximation of various channel-base currents waveshapes by analytical functions is very important, so as evaluation of their parameters.
There are many possible applications of evolutionary algorithms in lightning research and modeling of lightning strokes.Particle Swarm Optimization (PSO) is used for evaluation of lightning return stroke models' parameters (return stroke speed, channel height, and the constant of exponential current attenuation factor), based on electromagnetic theory relations and measured magnetic flux density in [10], whereas remote electric fields are used for estimation of lightning current waveforms and return stroke velocity profiles in [11,12].Genetic algorithm (GA) is proposed in [13] and Powell algorithm combined with PSO in [14] for identifying channel-base current parameters in the sum of two Heidler's functions [15], but only for single-peaked waveshapes of typical negative subsequent strokes.GA is also used for evaluation of parameters in the same function based on measured triggered lightning currents, and for the estimation of reflection coefficients [16] based on lightning discharge currents measured at tall instrumented towers, but also for the single-peaked channel-base current.An inversion method based on time series neural network and backpropagation neural network is presented in [17].
Approximation of negative first stroke currents is more difficult than of negative subsequent stroke currents, due to multiple peaks in waveshapes.For approximations of such waveshapes, linear combinations of a few functions from literature are usually used, as double-exponential function [18], Pulse function [19], Heidler's function, and others.A sum of seven Heidler's functions is used in [20] to approximate first negative strokes currents as measured in Brazil [21].Such waveshapes are approximated, for example, in order to calculate overvoltages induced in transmission lines [22].There are in total 28 parameters needed for that sum.Two-peaked analytically extended function (2P-AEF) uses just 7 parameters for the same waveshape [23].The accuracy is estimated in [20] from the current's peaks values obtained by the sum of seven Heidler's functions, whereas in the case of 2P-AEF, these peaks are chosen as exact values at the corresponding time moments.Besides, Npeaked AEF is a linear combination itself, so its number of peaks and parameters determine the accuracy of the obtained approximation.It is suitable for a variety of multipeaked waveshapes.
Application of GA to estimation of 2P-AEF parameters is presented in this paper.This function is based on the new channel-base current function [24], used to approximate the IEC 62305 Standard currents and other typical lightning currents [25].Marquardt least squares method is used for estimation of its parameters in [26].In this paper, 2P-AEF parameters are obtained by GA and given for an assumed negative first stroke current [27], but also for the measured currents waveshapes at Monte San Salvatore in Switzerland [1] and at Morro de Cachimbo Station in Brazil [21].Parameters of 2P-AEF are determined also for the typical first negative stroke current [1] and compared to results obtained by other functions.Triggered lightning current waveshape measured in artificially induced lightning discharge [28] is also approximated by 2P-AEF and its parameters are optimized by GA. 2P-AEF is described in Section 2 and application of GA to estimation of its parameters for the assumed two-peaked current waveshape in Section 3. Measured negative stroke currents are approximated in Section 4, the typical negative stroke current in Section 5, and triggered lightning current in Section 6.In order to validate the results of the obtained approximations, Feature Selective Validation (FSV) method [29,30] is applied in this paper, adopted by the IEEE Standard 1597.1 [31] for validation of computational electromagnetics results, computer modeling, and simulations.GA proved to have better results than least squares method used for the Two-rise function (TRF) [23].

Two-Peaked Analytically Extended Function
Analytical functions used for approximation of lightning currents with satisfying accuracy usually have a few parameters to be determined.The simplest function for representing pulse waveshapes is double-exponential (DEXP) function [18].First derivative of this function at  = 0 is not of zero value, as in realistic cases, which results in certain difficulties in lightning electromagnetic field computation.Its convex waveshape from zero to the maximum value is also not suitable for representing realistic lightning currents.
Heidler's function [15] has a concave rising part of the waveshape and first derivative of zero value at  = 0.It is given by for the current  0 , time constants  1 and  2 , degree  (usually chosen from 2 to 10), and the current peak correction factor  = exp[−( 1 / 2 )( 2 / 1 ) 1/ ].It is often used in the literature and adopted in the IEC 62305 Standard [8] for representation of typical lightning currents.A sum of two Heidler's functions is used for approximation of singlepeaked lightning currents, as suggested in [32], and used in [10][11][12][13][14].For the measured negative first stroke currents, a sum of seven Heidler's functions is needed [20].Linear combination of Heidler's function (for  = 2) and DEXP function is used in [33] for the negative subsequent stroke current and given by In order to approximate waveshapes with two distinct rise portions, a combination of power and exponential functions from [34] is simplified in [27], so the following form is obtained: for , parameters , , , , ,  1 ,  2 , current  0 , and peak correction factor .
Pulse function [19] is given by the following expression: for the current  0 , time constants  1 and  2 , degree , and the current peak correction factor 2P-AEF is given by the following expression: for the first peak  1 at  1 and the second peak  1 +  2 at  2 .There are 2 1 parameters to be determined ( 1 and  1 , for  = 1, . . .,  1 ) for the first interval (from zero to the first peak), 2 2 parameters ( 2 and  2 , for  = 1, . . .,  2 ) for the second interval (from the first to the second peak), and 2 3 parameters ( 3 and  3 , for  = 1, . . .,  3 ) for the last interval (after the second peak).The following three relations are valid: which results in  = 2( 1 +  2 +  3 ) − 3 unknown parameters to be determined for 2P-AEF( 1 ,  2 ,  3 ).For the two-peaked waveshapes, satisfying accuracy is obtained by 2P-AEF(2,1,2), for  1 = 2,  2 = 1,  3 = 2, and just  = 7 unknown parameters to be estimated by GA.Greater   increase the accuracy in i-th interval, at the expense of increased number of parameters to be determined.Different techniques may be employed for evaluation of these functions parameters.Among procedures based on least squares fitting are nonlinear curve fitting method (NLCF) [19] applied to Pulse function parameters and Marquardt least squares method [35] applied to DEXP [36], Heidler's function [37][38][39], Pulse function [39], and NP-AEFs [26].GA is applied in this paper for evaluation of 2P-AEF's parameters for the typical and measured waveshapes of negative first strokes.

Application of GA to Estimation of 2P-AEF's Parameters
GA is an evolutionary algorithm which has been applied for a few decades already [40] to a great variety of problems in different domains.A problem plays the role of an environment with an initial population of individuals generated at random or heuristically.An individual represents a possible solution to the problem with a certain degree of adaptation to its environment.Biological principle of natural evolution is used for artificial selection based on the computed fitness function.The candidate solution is encoded in its "genome" and new genetic material is introduced in each generation as in natural evolutionary mechanisms.Like in nature, better solutions to the problem are obtained, and GA is possibly leading to the optimal ones.One of the possible applications is curve fitting [41].
GA is applied in this paper with a goal to optimize 2P-AEFs parameters, so that the resulting waveforms well approximate measured and typical negative first stroke currents.Fitness function ff is the L-infinity norm described as for steps  = 1, .An assumed two-peaked waveshape [27], based on (3), is given in Figure 1.It is chosen to be approximated by 2P-AEF(2,1,2), for the first peak  1 = 16.5 kA at  1 = 3.3 s and the second peak  1 +  2 = 30kA at  2 = 4.6 s.In order to do so, lower and upper bounds of all parameters being optimized are given in Table 3.Values of 2P-AEF parameters obtained by GA are denoted by GA 2P-AEF, whereas results obtained by least squares method for TRF [23] are denoted by LS TRF in Table 3. Parameters values [27] corresponding to the function given by (3) are also shown in Table 3.These three current functions are presented in Figure 2. Amplitude Difference Measure (ADM) obtained by 1D FSV Tool [42], which is applied to the assumed current from [27] and its approximation by 2P-AEF(2,1,2) using GA, shows the quality of this approximation.The corresponding Grade-Spread chart is presented in Figure 3. Grade and Spread values obtained by 1D FSV Tool may range from 1 (best quality) to 6 (worst quality).ADMc is probability density function showing the proportion of pointby-point analysis in Figure 4(a) in order to check the measure of approximation's confidence.ADMi is showing point-bypoint comparison of the amplitude differences in Figure 4(b).
In the -axis shown in Figure 4(a), EX stands for excellent, VG for very good, G for good, F for fair, P for poor, and VP for very poor approximated data.The EX score (for 92.5%) and VG score (for 7.5%) of this approximation are confirmed by 1D FSV Tool.For smaller values of Grade (rated from 1 to 6), the quality of approximation is better, and for smaller values of Spread (rated from 1 to 6), the reliability of results Figure 1: An assumed two-peaked current waveshape from [27] given by ( 3).[27] and its approximation by 2P-AEF(2,1,2) using GA for parameters evaluation.[27] and its approximation by 2P-AEF(2,1,2) using GA for parameters evaluation. is higher.Results obtained by 1D FSV Tool for GA 2P-AEF approximation are Grade value 1 and Spread value 1 (Figure 3) for the threshold set at 85% by default.These two values (1,1) stand for the best quality of approximation.

An assumed current function
The first derivative of the current (3) for the assumed waveshape [27], TRF's derivative, and 2P-AEF's derivative for parameters obtained by GA are given in Figure 5.

Approximation of the Measured Negative First Stroke Currents by 2P-AEF
Previously described GA procedure is employed for optimizing 2P-AEFs parameters in order to approximate the measured negative first stroke currents.For all the examples, GA uses exact values for the current peaks and corresponding time moments as the input data (peaks are  1 at  1 and  1 +  2 at  2 ).As the first example, two-peaked waveshape of the negative first stroke current measured at Monte San Salvatore [1] is approximated by 2P-AEF(2,1,2) with 7 parameters estimated by GA.Parameters values for Example 1 are given in Table 4. Waveshape of the measured current denoted by MSS-FST#2peaks in [20], approximated by the sum of seven Heidler's functions (1), is the goal function for GA.Parameters of TRF are obtained by the least squares method [23] and denoted by LS TRF in Table 4.

International Journal of Antennas and Propagation
As the second example, the negative first stroke current measured at Morro do Cachimbo station in Brazil [21] is approximated by 2P-AEF(2,1,2) with 7 parameters estimated by GA.Parameters values for Example 2 are given in Table 4. Waveshape of the measured current denoted by MCS-FST#2peaks in [20], approximated by the sum of seven Heidler's functions (1), is the goal function for GA.Parameters of TRF are obtained by the least squares method [23] and denoted by LS TRF in Table 4. Results of these approximations are given in 100 s in Figure 8, and first derivatives of these currents are given in 20 s in Figure 9. Grade-Spread chart, obtained by FSV code which is applied to the sum of seven Heidler's functions representing MSS FST#2peaks results from [20] and its approximation by 2P-AEF(2,1,2) using GA, is presented in Figure 10.ADMc confidence histogram is given in Figure 11(a), ADMi point-by-point comparison of the amplitude differences in Figure 11(b).EX score is obtained (94%).For the approximation of amplitudes by GA 2P-AEF, the following scores are obtained: EX for 94%, VG for 5%, and G for 1% of results.Analysis by FSV code results in the Grade value 1 and Spread value 1 which stand for the best quality of  approximation.For MCS FST#2peaks current from [20] and its approximation obtained by 2P-AEF(2,1,2) using GA, Grade value is 1 and Spread value 1, as given in Figure 12.Amplitude approximation has EX (89%) and VG (11%) scores as given in Figure 13(a), and ADMi is as in Figure 13(b).

Approximation of the Typical Negative First Stroke Current by 2P-AEF
Typical waveshape representing 88 measured negative first stroke currents [1], normalized to the maximum current value, is approximated by and presented in Figure 14.As 50% of measured negative first strokes currents  [20] and its approximation by 2P-AEF(2,1,2) using GA for parameters evaluation.exceed 30 kA [1], it is often chosen as typical maximum current for this waveshape (which is the second peak   =  1 +  2 = 30 kA at  2 = 15.9 s).
DEXP function parameters are obtained in [36] by Marquart least squares method in order to approximate current waveshape of typical negative first stroke from [1].However, due to a very steep rise of DEXP function at  = 0, unlike slow rise of the approximated waveshape, this could not result in an accurate approximation.Certain improvement in the rising part is obtained by using Pulse function and is somewhat better by using Heidler's function, both optimized by Marquart least squares method [39].Such a waveshape may be approximated only by more complex functions or by a linear combination of usually used functions.However, 2P-AEF(2,1,2) is suitable for approximation of this waveshape with just 7 parameters to be determined, as in previous examples.2P-AEF's parameters are obtained by GA and their values are given in Table 4 for Example 3. Waveshapes of the normalized typical negative first stroke current [1] and its approximations by DEXP, Pulse, Heidler's function, and 2P-AEF(2,1,2), for parameters evaluated by GA, are presented in Figure 14.
Better approximation in the decaying part may be also obtained by 2P-AEF with more parameters in the third interval (after the second peak), for example, by 2P-AEF(2,1,3) with  = 9 parameters or 2P-AEF(2,1,4) with  = 11 parameters.A simple form of the function is preferable in many applications, if it is still able to correctly approximate the measured waveshape, so it was the reason to choose the smaller number of parameters.

Approximation of Rocket-Triggered Negative Stroke Currents by 2P-AEF
The approximation of rocket-triggered lightning current [28] was dealt with in [38,39].Heidler's, Pulse, and DEXP functions were optimized by using Marquart least squares method in order to represent this set of measured data.Heidler's function is better than the other two in representing the rising part of this current waveshape, but the decaying part is not represented well by any of these functions, due to its two-peaked waveshape.These waveshapes are given in Figure 15.However, the maximum current value and time to half of the maximum value are well approximated.2P-AEF(2,1,2) is applied for approximation of this waveshape and its parameters are obtained by GA.The parameters values are given in Table 4 (Example 4).This function is also presented in Figure 15 and denoted by 2P-AEF.
If compared to other functions, the waveshape obtained by 2P-AEF's approximation has better accuracy for just 7 parameters determined by GA.Accuracy of this approximation may be also improved by introducing more terms in the decaying part (for  3 ≥ 3).

Conclusion
There are many possible applications of evolutionary algorithms in lightning research [10-14, 16, 17, 43].Among these, genetic algorithm application to estimation of function parameters for curve fitting is presented in this paper.Measured multipeaked waveshapes of lightning currents and typical lightning current waveshapes of the negative first strokes are approximated in this paper by 2P-AEFs with parameters estimated by GA.If compared to other functions DEXP [36] Pulse function [39] Heidler's function [39] Time t (s) used for this purpose, better accuracy of approximation is obtained by 2P-AEF for a smaller number of parameters to be determined.Analytical expressions for 2P-AEF's first derivative, integral, integral of the square of this function, and its Fourier transform are given in [23][24][25] which makes 2P-AEF suitable for time and frequency domain calculations of lightning electromagnetic fields and lightning induced effects.Genetic algorithm is applied to the estimation of 2P-AEF's parameters in approximation of an assumed current waveshape [27], and the same procedure is used in cases of measured lightning currents of negative first strokes.For 2P-AEF and its parameters evaluated by GA, the quality of approximation is confirmed by FSV tools [29][30][31]42].There

2P-AEF
Measured data [28] DEXP [36] Pulse [39] Heidler [39]  is an excellent agreement of results obtained by 2P-AEFs and the goal functions.Currents derivatives waveshapes are also compared.Further research will aim at approximation of both currents and currents derivatives waveshapes by using GA.
In comparison to other functions from the literature, the best quality of 2P-AEF's approximation is obtained for the typical negative first stroke [1] which is useful in lightning models verification and lightning protection studies.Measured rocket-triggered lightning currents are also well represented by 2P-AEFs.
The main advantage of this function is that its N-peaked form may be used for approximation of various multipeaked waveshapes.For such waveshapes, N peaks at the corresponding time moments are chosen prior to application of International Journal of Antennas and Propagation GA.Greater numbers of terms in the intervals between the peaks improve the accuracy of approximation at the expense of increased number of parameters to be determined.However, a procedure applying GA for optimization of NP-AEF's parameters is generalized, so it may be useful in curve fitting of any waveshape.

Table 1 :
Parameter settings for the genetic algorithm.
. ., , whereas   are function's goal values and  AEF values obtained by 2P-AEF.The goal values may be obtained either by measurements or by survive to the next generation).The stopping criterion is the number of generations, that is, maximum number of iterations after which the procedure is stopped.In order to assist GA, upper and lower bounds (UB-upper bound; LB-lower bound) for each parameter being optimized are defined.These search regions are given in Table2.It should

Table 2 :
Search regions of the genome variables.

Table 3 :
(3)er and lower bounds of the genome variables and parameters values for 2P-AEF, TRF, and current function(3).