Optimization of substation grounding grid design for horizontal and vertical multilayer and uniform soil condition using Simulated Annealing method

Grounding systems are critical in safeguarding people and equipment from power system failures. A grounding system’s principal goal is to offer the lowest impedance path for undesired fault current. Optimization of the grounding grid designs is important in satisfying the minimum cost of the grounding system and safeguarding those people who work in the surrounding area of the grounded installations. Currently, there is no systematic guidance or standard for grounding grid designs that include two-layer soil and its effects on grounding grid systems, particularly vertically layered soil. Furthermore, while numerous studies have been conducted on optimization, relatively limited study has been done on the problem of optimizing the grounding grid in two-layer soil, particularly in vertical soil structures. This paper presents the results of optimization for substation grounding systems using the Simulated Annealing (SA) algorithm in different soil conditions which conforms to the safety requirements of the grounding system. Practical features of grounding grids in various soil conditions discussed in this paper (uniform soil, two-layer horizontal soil, and two-layer vertical soil) are considered during problem formulation and solution algorithm. The proposed algorithm’s results show that the number of grid conductors in the X and Y directions (Nx and Ny), as well as vertical rods (Nr), can be optimized from initial numbers of 35% for uniform soil, 57% for horizontal two-layer soil for ρ1> ρ2, and 33% for horizontal two-layer soil for ρ1< ρ2, and 29% for vertical two-layer soil structure. In other words, the proposed technique would be able to utilize square and rectangle-shaped grounding grids with a number of grid conductors and vertical rods to be implemented in uniform, two-layer horizontal and vertical soil structure, depending on the resistivity of the soil layer.

Introduction A well-designed grounding system is an essential factor in ensuring a safe operation of electrical systems which plays a major role in ensuring electromagnetic affinity, human protection and wellbeing, and the reliability of devices. The goal of an ideal grounding system is to reduce construction costs and time involved while satisfying the safety threshold parameters. The safety threshold parameters which consist of the grid impedance, step, and touch voltages are used to measure the behavior and safety level of a grounding system. These safety threshold parameters are computed using the equations from IEEE 80 [1]. For a grounding system to be declared safe, the magnitudes of the safety parameters must be below the threshold values. The safety threshold parameters are affected by the grounding grid geometry, grounding design parameters, and the electrical properties of the soil on which the grounding system is installed. Computerized grounding analysis in uniform and two-layer soil categories has become extensive, owing to the enhanced precision, speed, and adaptability provided by computer use. To address the optimization challenge, the traditional trial-error approaches become extremely time-consuming. Therefore, some non-traditional probabilistic search algorithms like the Genetic Algorithm (GA) [2][3][4], Particle Swarm Optimization (PSO) algorithm [4,5], and Harmony Search (HS) algorithm [6], have been used to optimize grounding systems in recent years. Table 1 shows the various algorithms that are built to optimize grounding impedance, increase grounding effectiveness, and reduce engineering costs. Most of the optimization was conducted in uniform and horizontal two-layer soil structure but none was analyzed in vertical soil structure due to the complex nature of the computation involved.

Optimization method Soil structure References
IEEE & Finite Element Analysis (FEM) method using ETAP 12.60 software Not mentioned [15] Gravitational Search Algorithm (GSA), Particle Swarm Algorithm (PSO) Uniform [5] IEEE method using ETAP software Horizontal twolayer [16] Optimal Compression Ratio (OCR) & fitting function Horizontal twolayer [17] Quasi-static composite image method & Green function, dual-port model Uniform [10] Hybrid method using chemical electrolyte ground rods, auxiliary wire mats, and ground enhancing material Horizontal twolayer [24] Simulated Annealing (SA) Uniform [23] Particle swarm optimization (PSO), space variable technique & harmony search algorithm (HS) using MATLAB Uniform [6] Genetic algorithm (GA), multi-objective particle swarm optimization (MOPSO) algorithm, Bees, IEEE, and Schwarz's equation Uniform [4] Optimum Compression Ratio (OCR), Particle Swarm Optimization (PSO) Horizontal twolayer [22] Steepest Descent Method using CDEGS Horizontal twolayer [25] Variable spacing techniques created OPTIMA function in MATLAB Horizontal twolayer [18] Economical Substation Grounding System Designer using MATLAB Uniform [26] ATP-EMTP & genetic algorithm Uniform [19] Genetic algorithm optimization Uniform [20] Simulated Annealing, Optimized Finite Element Method Uniform [9] https://doi.org/10.1371/journal.pone.0256298.t001 most common strategies employed. Table 2 shows the pros and cons of various optimization techniques. Analytical techniques in [15][16][17][18] are presented to address optimal grounding system configuration issues. Analytical techniques have the disadvantage of being unable to solve multi-objective optimization issues. To address the grounding system challenge, several researchers employed the GA method [2,4,19,20]. Wesley et al. [2] proposed a method by using GA to acquire soil characteristics that may be expressed in a multilayer structure. The method is based on an apparent resistivity curve generated from soil resistivity measurements. However, GA is poor in tackling the problem of a restricted local search space, is easily overestimated, and relies on the initial population for its solution. Besides, the grounding system optimization problem in [6] is solved using the HS algorithm. However, it has a low convergence rate and a low precision of optimization. The HS algorithm could not directly deal with constraints, thus a variety of constraint-handling approaches should be employed to help with the optimization process [13,21]. Furthermore, the PSO method [5,6,22], is utilized to optimize the design parameters of the grounding system, however the method has difficulty with discrete optimizations. The PSO algorithm has a weakness since it converges prematurely after being trapped in local optima, which it mistook for global optima. When the PSO algorithm is used to a complicated multi-dimensional situation, it becomes difficult to escape the local optima and achieve the global optima due to specific restrictions. The SA method is used to solve the grounding design parameters issue in [9,23]. However, it has a slow convergence speed and parallel computing is challenging. From all of the optimization methods that have been compared, the SA algorithm is proposed in this paper due to the discrete nature of grounding design parameters and the complexity of the problem which will increase with complex soil properties in which the grounding grid is installed. As a result, the optimization process requires a strong global search capability, with the SA algorithm being able to avoid early convergence to local optima and diverge the solution when compared to other optimization methods. The proposed algorithm is validated against results simulated using the Current Distribution, Electromagnetic Fields, Grounding and Soil Structure Analysis (CDEGS) software and optimization techniques in Refs. [27] and [26].

Soil apparent resistivity
As previously stated, the ability of this optimization method for grounding design parameters to find an optimal solution in a variety of soil conditions is its main significance. The soil's apparent resistivity is one of the most critical soil characteristics that influence the performance of a grounding system under various soil conditions. The soil constants and empirical

Analytical method
The approach is easy and requires minimal work.
Only deals with optimization problems with a single-objective [15][16][17][18] Genetic Algorithm (GA) The convergence speed is fast and the versatility is strong Complex, poor at handling the problem of restricted local search space, easy to overestimate, and its solution is dependent on the initial population [2,4,19,20] Particle Swarm Optimization (PSO) Has a good memory, fast training speed, simple method It is easy to fall into local optimal, inadequate discrete optimization issue management [5,6,22], Harmony Search (HS) Reasonable execution time, simple The convergence that occurs too soon [13,21] Simulated Annealing (SA) Strong global search capability Convergence speed is slow, parallel computing is challenging [9,23] https://doi.org/10.1371/journal.pone.0256298.t002 formulas of the uniform soil, horizontal two-layer soil, and vertical two-layer soil model used to calculate the apparent resistivity are based on data available in [28,29], as seen in Tables 3  and 4. The calculated apparent resistivity will be utilized as the input for the optimization process of the grounding system in different soil conditions.

Simulated Annealing algorithm
The Simulated Annealing algorithm is an intelligent optimization tool, explicitly tailored towards the optimization problem and has a series of benefits relative to the conventional optimization approach. The significance of the results is a set of good solutions instead of a single solution. It offers an alternate benefit, so it is particularly ideal for coping with difficult nonlinear engineering optimization issues. Simulated Annealing is a probabilistic search algorithm for an iterative adaptive heuristic that could be used to solve various nonlinear issues and compensate for non-differentiable and sometimes intermittent features. The probability of Simulated Annealing algorithms on achieving optimal solution could be higher. The optimal global solution can be achieved, with increased precision, global convergence, conditional parallel processing, and broad adaptability, which can accommodate multiple types of variables in optimization without any secondary knowledge on the objective function and thus no specifications for the control function.
To begin with a simple grounding model, the design parameters in the analysis of grounding performance are usually chosen in the general order or based on the physical grounding system available at the site. As a general practice, the number of grid conductors in the X and Y directions (N x and N y ) as well as vertical rods (N r ), are increased to obtain a lower grounding impedance and below threshold values for step and touch voltages. A design engineer Table 3. Empirical formulas for apparent resistivity in different soil conditions [28].

Soil condition Soil Constant Equations
Vertical two-layer soil ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi  might use as many grid conductors and vertical rods of different lengths as possible to protect a grounding system. However, a good grounding system should not be not only efficient but also cost-effective. Therefore, an optimization algorithm to reduce the engineering costs of a grounding system, easy modelling of grounding system in any engineering software, and lesser computation time is developed using the Simulated Annealing (SA) method. The limitations of grounding grid design parameters and safety thresholds are referred to from relevant equations and standard data from IEEE 80 [1] for developing the optimization algorithm. The design parameters of the grounding grid consist of the X-direction number of horizontal conductors, (N x ) and Y-direction number of horizontal conductors (N y ), and the number of vertical rods, (N r ).

Objective function
The purpose of grounding design optimization is to reduce the objective function, which is material and installation costs while maintaining safety. The cost of the grounding grid system can be broken down into four major components, with the cost function computed using (1) [6].
1. Cost of grounding grid conductor (C 1 ), it is associated with the size and length of the grid conductor.
2. Cost of burying the conductor (C 2 ), it is associated with the depth of grid buried, length of grid conductor 3. Cost of welding (C 3 ), associated with the number of conductors on each side.

Cost of rods (C 4 ), it is associated with the number of vertical rods.
Cost:Func: Where, A is the total area enclosed by the ground grid, m 2 ; N x is the number of grid conductors in x-direction; N y is the number of grid conductors in y-direction; L x is the length of grid conductor in x-direction; L y is the length of grid conductor in y-direction; h is the depth of grid conductor; N r is the number of rods.

Safety constraints
Computation of safety threshold parameters. A stable grounding system should be capable of sustaining the actual mesh and step voltages under the equivalent tolerable level within a substation. Equivalent grounding system impedance should be sufficiently low to ensure the more fault current disperses across the grounding grid into the earth to allow the grounding systems to be more secure. The optimization algorithm presented in this paper can be implemented for both 50 kg and 70 kg of body weight, although the results presented in this paper only demonstrated the step and touch voltages calculated for 50 kg only as an example for a worst-case scenario. The main focus should be given to the detail wherein any circumstances, the actual mesh and step voltages must not surpass the threshold values given in (2) and (3) [1] which define the maximum value of voltages that human beings can withstand in the case of an accidental step voltage and an accidental touch voltage scenario when designing the grounding system. Where E step50thres = tolerable step voltages for a person weighing 50 kg; E touch50thres = tolerable touch voltages for a person weighing 50 kg; E s = calculated step voltages; E m = calculated mesh voltages For the grounding system to be protected, the measured voltages present on substation earth as in (4) and (5) and grounding resistance value calculated using (6) must be lower than the threshold values provided by (2) and (3) for voltages and 5Ω for grounding resistance. Fig  1 displays a flowchart explaining the IEEE Grounding System Design approach. Detailed instructions and related equations for the design of the grounding system are available from [1].
Where ρ s = soil surface resistivity (Ω.m) t s = duration of a shock for determining allowable body current (s) C s = surface layer de-rating factor ρ = soil resistivity (Ω.m) L t = total length of conductors in the grounding system (m) I g = maximum grid current (A) K i = correction factor of grid geometry K s = step voltage spacing factor L s = total length of conductor for a step voltage (m) K m = touch voltage spacing factor L m = total length of conductor for touch voltage (m) h = depth of grounding grid buried (m) A = area of the grounding system (m 2 ) D = spacing between conductors (m) n = geometric factors of grounding The limitations of R g, E step and E m are determined in (7)-(9) by the protection criterion of equipment and systems concerning the remote earth which usually considers the insulation level and the performance of protection devices and systems. In most substations, R g is in the range of 5 Ω. Based on the above descriptions, the safety constraints can be expressed as the following formulae: Solution algorithm. The use of the Simulated Annealing algorithm to lessen the design cost while ensuring grounding grid safety is shown in Fig 2. The following Steps 1-11 are the

PLOS ONE
Optimization of substation grounding grid design using Simulated Annealing method details of the grounding optimization processes of the Simulated Annealing algorithm. These steps are similarly applied to calculate the touch voltage and grid impedance using relevant equations.
Step 1. Input system data Set the number of soil layers, soil resistivity (ρ) for uniform soil, soil resistivity (ρ 1 ) and (ρ 2 ) for two-layer soil, the total length of the grid and rod conductors (L t ), depth of grid (h), fault current (I g ), and so on.
Step 2. Calculate the apparent resistivity of the soil (ρ a ) using the equations in Tables 3 and 4 according to the soil structures where a grounding system will be placed.
Step 3. Set the initial configuration of a grounding grid.
X-direction number of horizontal conductors, N x (nrow) and Y-direction number of horizontal conductors, N y (ncol) of a grounding grid design. Set the lower [lb] and upper boundary [ub] for the nrow and ncol of the grounding grid.
Step 4. Calculate safety threshold parameters and initial cost function Step voltage threshold value (E step50thres ) is calculated using (3) while the step voltage of the current grounding system (E s ) using (5). The initial grounding cost (f c ) is calculated using (2).

Set 5. Generate a feasible solution
A new trial solution of ncol and nrow will be generated randomly, ncol and nrow are random numbers between [lb, ub].
Step 6. Calculate f cnew using new solutions of ncol and nrow.

Step 8. Check loop criterion
If Δf is � 0, the current solution is approved and the optimal solution estimation is changed to the adjusted solution. Otherwise, the existing solution is unchanged. 0 is the threshold of function tolerance which if it is exceeded, the iterations of the solver will stop and move on to the next step.
Step 9. Check safety constraints Step voltage of the current grounding system (E s ) is calculated using the new solution of nrow and ncol and compared with (E step50thres ). If the constraint is violated, go to Step 5. Otherwise, proceed to the next step.

Step 10. Check stop criterion
If E stepnew �E step50thres the current solution is approved and the optimal solution estimation is changed to the adjusted solution. Otherwise, the existing solution is unchanged. Repeat Steps 5-10 above, until loop iteration steps fulfill the criteria.
Step 11. Output the optimal configuration

Validation of the Simulated Annealing algorithm
Various parameters have many impacts on the efficiency of the grounding system. The Xdirection (N x ), and Y-direction (N y ) number of grid conductors, area of the grid (A), number of vertical rods attached to the main grid (N r ), and spacing between conductor's burial depth (h) were defined as the most prevalent and controllable parameters affecting the output of the grounding system [30]. A detailed explanation of the response of these parameters on the efficiency of the grounding system can be found in [31]. Up to two of these parameters can be the changing variables when deciding on the optimization issue of grounding system design. Optimization is implemented to make the design reliable and cost-effective in such a manner that the technological efficiency of the grounding system in terms of human protection and safety of equipment is never compromised and additional cost advantages are incurred by meeting all the essential electrical safety specification criteria.
To validate the results given by the SA algorithm, a simple design problem is solved using Matlab software and the MALZ module of CDEGS which is one of the grounding and electromagnetic analysis software packages available. The input data of the grounding system design parameters for uniform and two-layer soil structures are shown in Table 5 is adapted from [26]. Tables 6 and 7 provide further information on the soil characteristics of horizontal and vertical two-layer soil structures respectively based on assumptions, as Ref [26] only presented a uniform soil structure. As a consequence, the uniform soil structure findings in Table 8 are validated using both CDEGS and Ref [26] whereas the results in Tables 9-12 are validated  between SA and CDEGS. From the tables, there is a close agreement with the results of both the calculated and threshold values of step and touch voltages obtained by the SA algorithm and CDEGS. It can be seen that the percentage error between the SA algorithm and CDEGS is in the acceptable range of 2.97% to 10.7% for uniform soil, 1.09% to 5.86% for horizontal two-layer soil, and 2.35 to 8.16% for vertical two-layer soil structure. As for the comparison between the SA algorithm and the results presented in Ref [26] for uniform soil, the percentage error is 0% for threshold values of step and touch voltages. The percentage error for the calculated step voltage and grid resistance is 4.87% and 10.89% respectively. A major difference can be seen for calculated touch voltage at 52.67% for comparison between Ref [26]. This could be due to the difference in the equation used in Ref [26] as the percentage error for touch voltage between CDEGS and SA is only 10%. Thus, this validates the Simulated Annealing algorithm's performance developed using Matlab software as most of the percentage error falls within an acceptable range. Table 5. Grounding system design input for uniform and horizontal two-layer soil structure [26].

Design Parameters Values
Grid

Optimization of number of horizontal conductors in X and Y directions and number of vertical rods
The input parameters needed for the optimization of grounding system design in uniform soil are shown in Table 5. The soil resistivity for uniform soil is 40 Ω.m. Table 13 shows the corresponding results for optimizing the number of horizontal conductors (N x and N y ) and vertical Table 7. Soil properties of vertical two-layer soil structure.

Design Parameters Values
Distance between four electrodes (a) 5 m Angle between line where four electrodes are located and perpendicular line to the boundary between layers 1 and 2 (β) 60 0 Normal distance between the first electrode and boundary between layers 1 and 2 (d)  rods (N r ) in uniform soil structure. The cost of grounding is not included in this study since it is well known that using a less or optimal number of grid conductors and vertical rods can minimize grounding system implementation costs while ensuring system protection. The grounding design optimization in horizontal two-layer soil structure is conducted using the input data of the grounding design parameters shown in Table 14 adapted from [27]. For two-layer soil model, ρ 1 = 400 Ω.m; ρ 2 = 200 Ω.m for ρ 1> ρ 2 , and ρ 1 = 400 Ω.m; ρ 2 = 800 Ω.m for ρ 1< ρ 2 . The corresponding results for optimizing the number of horizontal conductors (N x and N y ) and vertical rods (N r ) in horizontal two-layer soil structure are presented in Tables 15 and 16. There is no research article on the optimization of the grounding system in the vertical soil layer available. Therefore, the results in Table 17 are given by comparing the results from SA and a CDEGS-simulated optimal grounding system.
The optimum number of grid conductors and vertical rods determined from the results would help in lowering the cost of grounding system implementation while maintaining the system's protection. Without jeopardizing the safety of the grounding system, the total amount of horizontal and vertical conductors has been reduced from initial numbers ranging from 9% to 35% for uniform soil, 57% for horizontal two-layer soil for ρ 1> ρ 2, and 33% for horizontal two-layer soil for ρ 1< ρ 2 , and 29% for vertical two-layer soil structure. These percentages refer to the number of horizontal conductors in the X and Y directions, as well as the number of vertical rods in the grounding system, being reduced from the initial number to the optimal number of conductors and rods. Besides, compared to the optimization results from Ref [26] for uniform soil, Ref [27] for horizontal two-layer soil, the optimization algorithm proposed in this paper gives better results in optimizing the number of horizontal conductors and vertical rods. This algorithm can optimize the design of the grounding system according to the updated system, even if the characteristics of the power system network are changed, such as the fault current or shock duration. As previously said, this algorithm employs the SA method to determine the best grounding design. Therefore, the algorithm would utilize the initial grounding design parameters to search for the best solution locally.

Conclusions
The overall cost of a grounding system was reduced by lowering the number of conductors in the X and Y directions, as well as the number of vertical rods connected to the grid, using the Simulated Annealing (SA) approach, without affecting the ground grid's protection. The proposed algorithm would be able to utilize square and rectangle-shaped grounding grids with a number of grid conductors and vertical rods to be implemented in uniform, two-layer horizontal and vertical soil structure which is yet to be done in any other research. The SA algorithm will diverge the parameters and delay premature convergence to the local optimum using its excellent local search capacity. The algorithm can successfully improve the optimization technique's convergence to prevent misleading the local optimum. Comparison has been made between the SA algorithm in Matlab and grounding analysis software (CDEGS) in different soil conditions, where the minimal differences between the two methods show the validity of the algorithm developed.