ASA Algorithm Combined with Current Simulation Method (CSM) for the Magnetic Induction under HV Power Lines in 3D Analysis Model

In this paper, a new hybrid approach to modeling the magnetic induction produced by HV overhead power line which combines the current simulation method (CSM) and adaptive simulated annealing algorithm (ASA) is discussed. The aim of the ASA algorithm is to find the optimal position and number of current loops used in bundles conductors for an accurate magnetic induction. Several parameters affecting the magnetic induction have been studied; it is observed that, taking into account the effect of conductor sag is much more interesting particularly at the mid-span length where the magnetic induction becomes very significant, the results also indicated that the maximum magnetic induction levels are less than the limits recommended by the International Commission on Non-Ionizing Radiation Protection (ICNIRP) standard for general public and occupational exposure. The calculated results are compared with those obtained from the COMSOL 4.3a Multiphysics software. A good agreement has been reached.


Introduction
Electricity is an essential factor to stimulate economic growth needed to improve and facilitate the daily quality of life of all human population, the continued rise in the world's population and with the tendency to concentrate in large cities. This has created an unprecedented demand in the consumption of electrical energy and has accelerated the installation of new high voltage power transmission lines near or in the populated areas. If there is no question of challenging the huge benefits made by electricity, the public is concerned more and more, for over thirty years about the potential effects of exposure to extremely low frequency electric and magnetic fields on the human health and the environment, this exposure is becoming increasingly important as technology advances and new electrical applications multiply.
However, the public awareness of the effects of electric and magnetic fields due to exposure to electric and magnetic fields from HV power lines became a major concern; the human population must meet the safety limits of exposure established by the international Standard such as the International Commission on Non-Ionizing Radiation Protection (ICNIRP). At extremely low frequency of 50 Hz, the limits of exposure set by this international standard, for general public exposure are 200 μT (magnetic field) and 5 kV/m (electric field). Respectively, for occupational exposure are 1 mT and 10 kV/m [1][2][3]. Therefore, it is important to accurately analyze the distribution of electric and magnetic fields in the vicinity of high voltage overhead power lines, in order to improve protection of human health and the environment [4][5][6].
In the last years, several publications have been made for the calculation and measurement of very low frequency electric and magnetic fields (ELF) created by the HV power lines. Most assume that power transmission lines are straight horizontal conductors, parallel to a flat ground, and parallel to each other. The sag due to the line weight is neglected or modeled by introducing an average height for the horizontal line [7][8][9][10][11].
In view of the above, this study aims to determine the levels of the magnetic induction created in the vicinity of High Voltage (HV) power transmission lines, using the current simulation method (CSM). However, this technique has rarely been used in the literature and the results reported with this method are almost very limited.
Successful application of this simulation method requires an appropriate number and position coordinates of the filamentary line currents, the values of filamentary currents are determined by satisfying the conditions of limits to a selected number of points in the sub-conductors. Once the values and positions of simulation currents are known, the distribution of the magnetic field in any region can be calculated easily [12][13][14][15].
In order to overcome difficulties and to improve the calculation accuracy, it is desirable to use an optimization technique of the adaptive simulated annealing (ASA) algorithm in combination with this method. ASA is a powerful global optimization algorithm, especially practical for nonlinear and / or stochastic systems. This algorithm is developed to statistically find the best global fitting of a nonlinear and non-convex cost-function over a Ddimensional space [16,17]. It should be noted that this calculation takes into account the effects of the catenary line, where the sag depends on the individual characteristics of the electrical line and environmental conditions. The catenary geometry effect is rarely considered in the literature. The performance of the proposed hybrid method will be verified with the COMSOL 4.3a Multiphysics software based on the Finite Element Method (FEM).

Model of Overhead Power Lines
The conductors of an overhead electrical line are not at all points at the same height along the span length (longitudinal axis). In fact, they regularly describe a catenary, where the sag depends on the individual characteristics of the electrical line and environmental conditions. Figure 1 depicts the basic catenary geometry for a single-conductor line [18][19][20]. Where, L: The distance between the two points of suspension (span); S: The sag of the conductors; hmin: The height of the conductor at mid span; hmax: The height of the conductor at the end of the line (foot of the pylon); have: The average height of the conductor.
The equation of the catenary shape of conductor placed in the (yz) plane is given by [21].
(1) z is the longitudinal position of the conductor about z axis, for a symmetrical line, you normally choose z = 0 at the midspan; δs is the solution of the equation below: (2) To calculate the height of the electrical line along z axis in a span length, the following equation can be used [22].
Some researchers, in the magnetic induction calculation in the vicinity of power lines assumes that the conductors are horizontal of infinite length, parallel to a flat ground and parallel with each other, and the sagging due to the weight of conductors is neglected (2D numerical model), taking into account an average height (see Figure 1) [22], this average height is given by the following expression. (4)

Current simulation method (CSM)
The idea of current simulation technique is very simple and is similar to the charge simulation method. This approach is adapted to the electrical conductor with bundles of several conductive elements called sub-conductors [12][13][14][15].
In this method, each sub-conductor current is simulated by a finite number n of filamentary line currents distributed on a fictitious cylindrical surface of radius Rj, the number and coordinates of the simulation currents are not arbitrary but based on the number of conductors and their arrangement in space.
The simulation currents for all the number of streams of filaments in all sub-conductors Ii, (i = 1,….., 3×n×m) must satisfy the following conditions [12][13][14][15]. To determine the unknown filamentary currents, a set of equations is formulated at a number of boundary points chosen on the sub-conductors surface to satisfy the boundary conditions as follows:   Pij is the magnetic normal field coefficient determined by the coordinates of the i th boundary point and the j th filamentary line current and is given by [13]. (7) Rij is the distance between simulation current point j and match point i; Rj is the fictitious radius of current filament simulation (see Figure 2). We find another expression [12,14] which uses Equation (8) to calculate the magnetic coefficient.
Once the set of equations (5) and (6) are solved for the unknown filamentary line currents, the deviation of the normal component of the magnetic field strength of the zero value is calculated at a set of checkpoints (match points) chosen on the sub-conductors surfaces, or we calculate the new magnetic potential coefficient corresponding to the chosen check points, the relative error is considered as a measure of simulation accuracy. When the simulation currents satisfy the boundary conditions at the checkpoints with a considerable accuracy, the magnetic field could be calculated at any point in space. On the other hand, we repeat this calculation by changing the place and / or the number of simulation current points and match points. The components of the magnetic induction at any point in the space around the power line can be calculated by the following equations [12][13][14][15].
Where, Aij is the magnetic vector generated by the threephase power line current, it is expressed as follows: For magnetic field calculation, the image of the filamentary current in a given sub-conductor is located at the depth De different from the sub-conductor height above ground given by [12,14].
Where, ρ is the electrical resistivity of the medium; f is the frequency of the source current. The resultant magnetic induction is calculated as [12,14].
In this analysis of magnetic induction produced by HV power line, the catenary form of the overhead power line conductor (sagging conductor) is taken into account as indicated in Figure 1 above; this 3D quasi-static analysis can supply to good magnetic induction evaluation.

Adaptive simulated annealing (ASA)
Simulated annealing is a probabilistic method proposed by Kirkpatrick, Gelett and Vecchi (1983) and Cerny (1985) to find the global minimum of a cost function, it based on the analogy of the physical process of annealing a metal. At high temperatures, the atoms are randomly distributed. With decreasing the temperature, they tend to arrange themselves in a crystalline state which minimizes their energy. Using this analogy, the algorithm generates randomly new configurations by sampling from probability distribution of the system. New configurations are generated by applying movements, which are randomly selected with a defined probability. A high temperature is given initially, when all movements are accepted even if they increase the energy. As temperature decreases, fewer movements are accepted, it shows that the algorithm converges to a global energy minimum if the temperature is reduced slowly enough .Adaptive simulated annealing (ASA) is a variant of simulated annealing (SA) algorithm, the difference lies in the fact that the algorithm parameters that control the temperature schedule and the selection of random step are automatically adjusted according to the progress of the algorithm [23][24][25][26][27].
The algorithm of the ASA is detailed as follows [28,29].
(1) -Initialization: An initial x is randomly generated, the initial temperature of the acceptance probability function, is set to the initial value of the cost function , the initial temperatures of the parameter generating probability functions, are set to 1. A control parameter c in annealing process is given, and the annealing times, for and , are all set to 0. (2) -Generating: The algorithm generates a new point in the parameter space with: (13) Where, Mi is a number determined by the difference between the upper and lower boundary of xi, qi is calculated as: (14) And vi a uniformly distributed random variable in [0, 1]. The value of the cost function is then evaluated and the acceptance probability function of is given by: A uniform random variable is generated in [0, 1]. If is accepted; otherwise it is rejected.
(3) -Reannealing: After every acceptance points, calculating the sensitivities: (16) Where, is the best point found so far, δ is a small step size, the n-dimensional vector ei has unit i th element, and the rest of elements of ei are all zeros. Let .
Each temperature Ti is scaled by a factor and the annealing time ki is reset: (17) Similarly, is reset to the value of the last accepted cost function is reset to , and the annealing time kc is rescaled accordingly: (4) -Annealing: After every Ngenera generated points, annealing takes place with: Otherwise, go to Step 2.
(5) -Finish: The algorithm is terminated if the parameters have remained unchanged for a few successive reannealings or a preset maximum number of cost function evaluations has been reached; otherwise, go to Step 2.
The objective function (cost function) used for optimization is given by the equation below: (20) Where, Aci is the magnetic potential calculated by the current points; Ami is the new magnetic potential estimated by the match simulation current points and nt is the total number of match points.

Finite Element Method (FEM)
The finite element method (FEM) is a numerical technique used to find approximate solutions of partial differential equations as well as of integral equations. The solution is based either on eliminating the differential equation completely, or rendering the partial deferential equations into approximating system of ordinary differential equations. The FEM analysis of any problem involves basically four steps [30][31][32][33].
(1) -To discretize the solution region into a finite number of sub-domains or elements; (2) -To derive the governing equations for a typical element; (3) -To assemble all the elements in the solution region; (4) -To solve the obtained equation system. In quasi static analysis, the distribution of the magnetic induction B is derived from the magnetic vector potential A [30][31][32][33][34][35][36].
When the magnetic vector potential is applied to the transmission lines wrapped space whose physical medium is air, there is no current density; this relationship reduces to the magnetostatic Laplace's equation as follows [30][31][32][33][34][35][36]: Since the calculation of magnetic fields in the transmission lines includes areas which are partially open and infinite, it is necessary to perform specific processing of the solution of equation (22), which is possible by using the FEM. From the variational principle, it can be proved that the minimization of the following energy functional gives the solution of equation (21) [ [36][37][38][39]: In FEM, the volume of the proposed region is divided into (m) small triangular elements where their sides form a grid with (n) nodes. Thus, the vector potential A at any point inside this element can be approximated according to the shape functions N [36][37][38][39].
The terms Nk (where k =i, j, m) are given by: Where, Δe represents the area of the triangular element, Ak is the potential at the k th node.
The solution of the magnetic field problem is obtained by minimizing the function of equation (22), it can be formally as follows: Through the solution of the nodal potential, the magnetic induction B can be determined from the calculated magnetic vector potential A, according to equations (21) and (24). Finally, in two dimensional problems (2D), equations in the x and y directions for magnetic induction can be written as [36][37][38][39].
Finally, according to the equation (28), we can calculate the resultant magnetic induction strength as below, The software package of COMSOL Multiphysics (version 4.3) based on the finite element method (FEM) can be used for modelling the main characteristics of electric and magnetic fields in many electrical applications (High Voltage overhead power lines). This simulation package solves linear and nonlinear equations in (1D, 2D, or 3D) models [32,[37][38][39][40].
The COMSOL Multiphysics steps may include the following [40,42]: -Creating the Geometry model and assign the material properties to the domain under study -Modeling Physics and Equations (Boundary and interface conditions) -Creating Mesh -Solution -Visualization of results In this study, we consider a three-phase HV overhead power line of 275 KV, with the arrangement and geometric details in the vicinity of the suspension pylon as shown in Figure 3, each phase of the power line consists of a bundle of two conductors separated by 30 cm with a radius 10 mm, the ground wires radius is 10 mm. the span length is 300 m, the maximum sag of the phase conductor is 8 m and 6 m for the ground wire. The three-phase currents have been assumed under a balanced operation with the magnitude of 500 A (with a nominal frequency of 50 Hz), the current in the ground wires is assumed to be zero; the earth is assumed to be homogeneous with a resistivity of 100Ω.m.

Results and discussion
First, one must look for the optimal values of the proposed method of calculating to have good simulation accuracy. This research is based on an optimization using the Adaptive simulated annealing (ASA) algorithm. The control parameters of the ASA algorithm and search intervals of the variables for the current simulation method are summarized in Table 1.    After multiple runs for the objective function optimization mentioned above in equation (20), once the algorithm terminates execution, the optimized function value and the optimal parameter values are obtained for the ASA algorithm, the solution converges to the minimum of the objective function, the ASA algorithm convergence with number of iterations is shown in Figure 6.
The search process at successive iterations with optimal solutions for the number and position of the filamentary line currents are represented in the Figures 7 and 8, respectively, where it becomes clear that this algorithm converge quickly to these values. The best solution obtained by the proposed ASA algorithm is incorporated into the proposed method and is summarized in Table 2.  Figure 9 shows the lateral distribution of the magnetic induction strength at mid-span point under the HV power line at a height of (1m) above ground level in different points along the right-of-way corridor, it can be seen that the maximum value of the magnetic induction is obtained under the middle phase conductor, and then decreases rapidly with increasing lateral distance in both sides of the right-of-way. The magnetic induction vertical component a less intense under the middle phase, and increases to a maximum value under the lateral conductor phase, then decreases progressively as one moves away from the power line center, in which it becomes very negligible. For the horizontal component the magnetic induction is maximum at the power line center and decreases to a much lower value in the immediate vicinity of the side conductor phase, from this distance; it increases slightly and then again decreases progressively with the lateral distance to reach much neglected values.

Position of filamentary line current [m]
Phase conductor Ground wire The instantaneous variation of the module and the magnetic induction argument with its horizontal and vertical components at mid-span length at a fixed point of observation (x = 0, y = 1) is shown in Figure 10. As time passes, the magnetic induction oscillates as a function of time along the axis (ox). The maximum value at this point is obtained at instant t = 0.0032 (s) for an argument θ = 0.4°. The contour lines of magnetic induction strength around the power line conductors at mid-span length in the (xy) plane are depicted in Figure 11. This difference in the levels shown in this figure (e.g., 5,10,15,20,25,30,50) is due to the variation of the coordinates of the calculation point. The intense lines are in close proximity to the phase conductors surface. The lateral distribution of the magnetic induction strength at mid-span point and around the pylon is presented in Figure  12. It can be easily noticed that the maximum value of the magnetic induction calculated at the mid-span is greater than the maximum value obtained at the pylon foot. This is certainly due to the height difference between the two adjacent pylons, the height of the conductors at mid-span point (ground clearance) is less than that of the conductors at the end of span length (pylon foot), the magnetic induction profile above ground level follows a nearly Gaussian shape.

Longitudinal Span [m] Magnetic induction [µT]
Sagging effect Average height Figure 13 illustrates the longitudinal profile of the magnetic induction under the middle phase conductor, at the point where the magnetic induction is maximum, It can be observed that, taking into account the effect of the sagging conductor of the power line, the magnetic induction strength around the pylon is much lower (3.93 μT), as one moves away from the pylon, the magnetic induction gradually increases to a maximum value under the mid-span point (8,85 μT), and then decreases again to the same lower value at the other end of the span length.
For the average height when the effect of conductor sag due to the line weight is neglected, the magnetic induction strength at 1 m above ground level is 7.33 μT; this value remains constant along the span length of the power line. This illustrates that taking into account the effect of conductor sag in the magnetic induction analysis is a very practical factor to model the real curvature of the power line; it plays an important role in determining the exact values of the magnetic induction. As a result, it can be seen that the maximum values obtained of the magnetic induction are below the limits set by the ICNIRP recommendations and therefore they do not present any appreciable hazard to human health [43,44]. Figure 14 shows the three-dimensional profile of the magnetic induction over an area equivalent to a right of way equal to 50 m either side of the HV power line center and a longitudinal span between the suspension pylons of 300 m. It will be noted that the concentrated level of the magnetic induction exists in a small area under and in the immediate vicinity of the conductors in the mid-span, and then it decreases slowly toward the pylons and even with the increasing of lateral distance from the power line center in both directions along the right-of-way corridor. The magnetic induction magnitudes are weaker near the pylons than between them, because the power lines are at their highest point above the ground level. .    As indicated above, the conductor sag is the principle factor affecting the magnetic induction strength. In what follows, it is worth recalling some other factors, which might influence the value of magnetic induction. The different effects of the variation of the conductor height above the ground and the spacing between the phases, also the observation point height of magnetic induction calculation are shown in Figure 15. From these results, the conductor height affects the magnetic induction strength under power line directly, when the conductor's height increases above ground level, the magnetic induction will decrease rapidly. Moreover, with increasing the spacing between two adjacent conductors, the magnetic induction increases slightly. On the other hand, an increase in the observation point height from the ground leads to increased in magnetic induction strength.
In order to know what type of circuit configuration creates a maximum magnetic induction, the profile given in Figure 16 illustrates the lateral distribution of the magnetic induction for different circuit configurations a single power line (see Figure 5). It is clearly shown that the horizontal configuration gives a higher magnetic induction than the other configurations because the three phase conductors are close to ground level, while the inverted triangular circuit produces lower magnetic induction to better cancellation effect of the three-phase currents. In addition, the vertical configuration gives almost similar magnetic induction profile to the inverted triangular circuit.
For double circuit lines shown in Figure 5 above (same phasing abc-abc), as can be seen in Figure 17, the vertical lines geometry causes higher magnetic induction values at all points over the lateral distance, on the other hand, the inverted triangular configuration represents smaller values of the magnetic induction under and close to the power line centre.
In double circuit lines, the phase arrangement has an important influence in the magnetic induction evaluation; it is very possible to organize the phase conductors to improve the cancellation. As an example, the magnetic induction for different phase arrangements in a vertical double circuit line with the same parameters is illustrated in Figure 18. As can be seen in this figure, the inverse phase arrangement (abccba) gives the lowest value of the magnetic induction for the different points around the power line centre, while the same phase arrangement (abc-abc) produces a higher magnetic induction than all other arrangements of phase conductors.
Finally, in order to verify the validity and reliability of the adopted method, the COMSOL Multiphysics (version 4.3a) can be used to ensure this analysis effectiveness. The magnetostatic module of this software was used to analyze the domain solution of this model in two-dimensional space. Meshing of the domain geometry with linear triangular elements generate by this software is shown in Figure 19. A fine mesh is required in the given region, in order to achieve the quality level required by a more precise simulation. Figures (20 and 21) show the lateral profiles of the magnetic inductions strength at mid-span length and pylon foot of the

Lateral distance [m] Magnetic induction [µT]
abc-abc abc-acb abc-bca abc-bac abc-cab abc-cba power line at 1 m above the ground level using COMSOL Multiphysics 4.3a software. It can be seen that the highest value of the magnetic induction is under the middle phase conductor, and then decreases rapidly with increasing in the lateral distance from the power line center. From these figures, it is evident that the magnetic induction strength immediately below the lowest point (mid-span) of the power line is significantly higher than in the proximity of the pylon and at some lateral distance from the power line. Therefore, the magnetic induction at 1 m above ground level is strongly affected by conductor sag since an increase in conductor sag will increase the magnetic induction strength, especially in mid-span where the phase conductors will have maximum sag. It's also important to remember that the proximity effect of the grounded tower which influences the magnetic induction distribution is thereby often ignored. A comparison of magnetic induction has been made between the adopted method and COMSOL Multiphysics 4.3a, based on the finite element method (FEM). Generally, the FEM constitutes the most powerful and accurate tool for multiphysic problems analysis, it is mathematically more rigorous, leading to more accurate results. On the other hand, the CSM uses a rather clever physical reasoning, but considerably very approximate. A visualization of the magnetic field profiles produced by these two approaches is given in Figure 22. A good agreement and a similarity between the results were found, which validates the accuracy of the proposed method. Concerning the comparison of the time execution between the two tools, it's important to mention that an Intel (R) Core (TM) i3-4005U processor with 2 physical cores and 4 logical processors has been used in this simulation, it was concluded that the CSM+ASA is approximately much faster 2 times than the COMSOL Multiphysics.

Conclusions
This paper presents a novel optimized technique for quasistatic modelling and analyzing the magnetic induction generated by HV power line based on the current simulation method (CSM) coupled with the adaptive simulated annealing (ASA). ASA algorithm is developed in order to determine the optimal arrangement (position and number) of filamentary current-loops to achieve a more accurate calculation. From the results, it is clear that the maximum value of the magnetic induction is located under the middle phase conductor then decreases by increasing the lateral distance from the power line center. The magnetic induction is strongly affected by several factors; the height of phase conductor and the spacing between the phase conductors, the geometric configuration of the phase conductors (single or double), the phase sequence arrangements for double circuit line can also affect the magnetic field distribution. The most Figure 19: Finite element discretization of the study domain given in Figure 3.   important parameter which affect significantly the strength of magnetic induction is the conductor sag, especially under the mid-span length where the magnetic induction becomes very important than that in the proximity of the pylon. It is noted from these obtained values that the ICNIRP standard is respected for the general public and occupational exposure. The results of the proposed method were compared with those obtained from the COMSOL 4.3a software. The results were virtually identical; this comparison is quite satisfactory and can correctly verify and validate the hybrid approach.