A Review: Evolutionary Computations (GA and PSO) in Geotechnical Engineering

This study briefly reviews the application of Genetic Algorithm (GA) and Particle Swarm Optimization (PSO) in geotechnical engineering since GA and PSO are widely used in civil engineering. The application of GA and PSO is studied in three popular families of geotechnical problems including unconfined seepage analysis, slope stability analysis, and foundation design. In each category the available results from different studies are reviewed and compared. The comparison of results shows the desirable accuracy in the predicting of optimal values in the process of analysis and design. The presented methods perform successfully in the reviewed problems. However, PSO predicts the optimum values in fewer numbers of iterations, which suggests higher performance in term of implementation/application.


Introduction
In geotechnical engineering the optimization techniques can be categorized as: 1) mathematical programming; 2) optimality criteria methods; and 3) heuristic search algorithms. In mathematical programming, the major characteristic includes the linear and nonlinear programming. In the linear programming, the objective function and the associated constraints are represented in a linear combination of design variables. The linearization of objective function or constraints is not always easy and if the linearization techniques are used, the error in the representation of linearized problem is inevitable [1] [2] [3]. On the other hand, the nonlinear programming is introduced which was developed for un-constrained nonlinear problems. The incidence of the optimal solution is shown by the proof of Kuhan-Tucker conditions (KT) [4] which are the necessary conditions for the justification of the optimal solution. However, the application of KT conditions is enormously difficult for the most of engineering problems [5].
The optimality criteria methods (OCM) are developed based on the combination of Kuhan-Tucker conditions from nonlinear mathematical programming and Lagrangian multipliers. In this approach, KT conditions support the necessary desires for the optimal solution while Lagrangian multiplier guarantees the satisfaction of constrains in the optimization problem. OCMs are used vastly in engineering problems including the continuous and discontinuous design variables [6] [7]. OCMs basically describe continuous design variables. However, by considering two major steps, this procedure is compatible for discrete variables, as well. First, the optimization problem is defined based on the assumption of continuous solution. Second, a set of discrete values are estimated based on the solution from first step [8]. However, two problems are associated with this approach; for instance, if OCMs are used in structural design while a single cross sectional property is chosen as a design variable and the other properties are expressed as a function of the design variable. First, the relationship between the design variable and cross sectional properties is not unique. Second, the selected discrete variables may lead to a different structural response which does not fit the prescribed constraints [5].
The techniques of the third group in the optimization problems are heuristic search algorithms which are relatively new optimization methods in geotechnical engineering. Evolutionary methods do not require an explicit relationship between the objective functions and constraints. Moreover, these methods are independent of objective function or constraint gradient. Instead, the set of random values for objective function is adjusted thorough an evolutionary procedure and subsequently, any violation of constraints is reflected in each iteration.
Regarding the gradient based methods, considerable part of computational cost is served in sensitivity analysis phase whereas applications of evolutionary techniques which are based on probabilistic searching algorithm save computational time [9]; this is because of evolutionary methods which do not require the gradient information and avoid sensitivity analysis. This property enables engineers to implement these techniques in the problems including discontinuity in the design variables space. In general, evolutionary methods work concurrently with a population of design points in the space of design variables instead of a single design point. These methods also can be easily implemented in discrete, continuous, and mixed optimization problems with minor adaptations. In addition, the open format for constraint statements and the possibility of defining multiple scenarios in the optimization process fascinated many researchers to implement heuristic search algorithms in engineering design optimization [10].
The safety in designs is highly respected in geotechnical engineering. On the other hand, contractors and private sectors prefer economical design whit minimum labor effort and construction cost. Subsequently, implementation of ro-bust optimization techniques which establish a trade-off between safety and total cost of geotechnical projects is necessary in practice.
In this study, we reviewed the two most implemented evolutionary methods in geotechnical engineering, Genetic Algorithm (GA) and Particle Swarm Optimization (PSO). The applications of GA and PSO in three common types of geotechnical problems including unconfined seepage analysis, slope stability analysis, and foundation designs are studied and subsequently, for each problem, the results from literature are regenerated, compared and discussed.

Genetic Algorithm (GA)
Genetic Algorithm which is based on the process of natural evolution is one the best-known heuristic search algorithms. GA has received substantial attention in engineering design optimization in the recent decades. The first time that GA was introduced in computer science goes back to sixties when a team of biologist tried to implement the process of evolution in nature in a computer code [11].
GA refers to any population-based algorithm which uses selection, crossover, and mutation over chromosomes to find the optimal solution. In general, a member of the population is called chromosome/genotype which is binary or real valued string. Following Barricelli [11] several types of GA have been introduced in optimization studies. However, one can simply define a given optimization problem in GA considering three main steps including [12]: Step 0: Initialization.
Generation of initial chromosomes in step 0 (genotypes) is the first step in GA algorithm. In the most practical problems genotypes are generated randomly and goodness of each chromosome is evaluated via objective function and associated constraints.
Step 1: the selection operator is introduced and applied to the current population to create an intermediate one. The initial and intermediate populations are same in the first step (step 0). However, in the subsequent iterations, the imposition of selection operator forms the intermediate population.
Step 2: crossover-mutation operators are implemented to the results from Step 1 to create next population. The term crossover is used when the generator (operator) forms a chromosome by combining the properties of each of two parental chromosomes. However, the term mutation is applicable when a new chromosome is formed by introducing small alterations to the properties of single parent chromosome [13].
There are many choices to design/encoding an optimization program based on the mentioned steps and the decision on design/encoding of GA depends on the nature of the problem [14]. The algorithm design depends on experience, the model specifications and associating the results from experiments with different heuristic search algorithms. However, a common design for GA with standard genetic operators can be defined as [15]: 1) Randomly generate an initial population including M chromosomes.
2) Calculate the fitness, ( ) F m , of all chromosomes in M.
3) Create the evolved population based on: a) Using proportional fitness selection for selected two chromosomes, m1 and m2. b) Apply crossover function to m1 and m2 to produce a new chromosome (m3).
c) Apply mutation function to m3 to produce m′ . d) Add m′ to the next population.

Particle Swarm Optimization (PSO)
Particle Swarm Optimization is a metaheuristic search algorithm introduced by Eberhart and Kennedy [16] to find optimal solution in engineering design optimization. PSO is based on the concept social models and swarm theories. The swarm consists of individual particles which mutually try to find the solution in the search space. In each iteration, individuals (particles) move toward the best solution which is experienced by them (Personal best) and concurrently to the best solution which is obtained by the other particles (Global best). PSO is established based on few or even no assumption on the search space; this feature enables PSO to search the optimum solution in a wide search space. In addition, PSO can be used in the optimization problems which are irregular, noisy, or dynamic [17].
Regarding the PSO algorithm, the swarm of particles walks through the n-dimensional search space of a problem including n degree of freedoms. and Global best solutions in the m th iteration [18]. In general, these parameters are determined based on tuning and trial-error observations.
With respect to Equation [1] and aforementioned concepts, PSO algorithm covers the following steps to find the optimum solution: 1) Random generation the initial position and velocity vector including 0 i x and 0 i v , respectively. 2) Evaluate the objective function with respect to predefined constraints.
3) Comparing the fitness value of each particle to its local best. If the current value is better than the previous one, then it is replaced with the new value. 4) Comparing a particle's fitness value to the population's best. If the current value is better than the previous one, it is replaced as the global best. In each iteration, the predefined maximum velocity and position control the changes in velocity and position of each particle.

GA and PSO in Seepage Analysis
Determination of the seepage path through the earth dams is one the most difficult problems in geotechnical engineering. This phenomenon is recognized as unconfined seepage problem since the location of phreatic line is not determined. Subsequently, it represents a boundary value problem where at least one of the boundaries is not determined. In order to find the location of phreatic line (free water surface) an iterative procedure is needed [19]. Varity of optimization techniques are used in unconfined seepage problems to find the phreatic line in the earth dams [14] [20] [21]).
In the first step, the governing equation of unconfined seepage problems is is the normal vector on the impermeable boundary).
2) Equipotential lines: Hydrostatic pressure is applied on this boundaries and the total head is constant along these boundaries. u φ and D φ are the total To summarize, the mathematical definition of governing equation and boundary conditions for a 2D seepage problem is defined: where j k shows the hydraulic permeability in saturated porous media, can be defined as [23]: Min , Subsequently, the mathematical constraints are introduced in the optimization procedure as follow: Fenton and Griffiths [24] analyzed the unconfined seepage problem using mesh deforming method based on nonlinear FEM. Their study represents an iterative procedure which adjusts the height of the nodes in the mesh geometry. However, Ouria and Toufigh [23] used linear FEM associated with Nelder-Mead simplex method. They approximated the phreatic line with a 4 th degree polynomial and minimized the goal function based on the conditions on the phreatic line. Shahrokhabadi and Toufigh [25] used the same procedure but they used Natural Element Method (NEM) and GA in their solution procedure. In fact, they used a mesh deformation technique in NEM which is not sensitive to mesh geometry and GA to avoid the local optimum solution. Moreover, they showed that the average error in optimization and number of iterations are noticeably affected by the initial guess in gradient based optimization algorithms while GA is not sensitive to initial guess. Figure 2 represents the comparison among the results using gradient based optimization method and FEM [23], nonlinear FEM analysis [24], and NEM-GA [25]. It is shown that linear FEM with Nelder-Mead simplex and NEM-GA lead to similar results for the seepage path. However, the obtained exit point based on nonlinear FEM and NEM-GA is identical and it is reported as 3.63 m.
Precise estimation of the Exit point (see Figure 1) is a challenging problem in unconfined seepage problems and considered as a significant criterion that de- In order to present more accurate study, Table 1 shows the exit points reported by different numerical and analytical studies. It is observed that MFS-PSO-TCF represents closer results to analytical solution in comparison with other numerical-optimal solutions. Regarding the presented results in Table 1,

MFS-PSO-TCF obtains the closest results to analytical solution introduced by
Ozis [31].
In the last part of study on application of evolutionary methods in unconfined seepage problems, the comparison between results from MFS-PSO-TCF and experimental results is studied [32]. Fu and Jin [32] made a tank of armor plates    whereas an evolutionary optimization method is implemented to find the unknown boundary (phreatic line).

GA and PSO in Slope Stability Analysis
Slope stability is a major problem in geotechnical engineering that estimates the failure potential of geo-materials (soil/rock) covering the slopes. The stability is considered as the balance between shear stress and shear strength which are due to soil mass movement and material resistance, respectively. The stability of soil mass covering a slope depends on both slope geometry (slope angle, back slope angle, pore water pressure and etc.) and loading conditions (climatic events, loading and lateral pressure, and etc.).
Slope failure comes about in the failure zone (critical surface) in which the soil strength is less than shear stress. There are numerous techniques to determine the critical slope surface in the classical geotechnical engineering. In order to analyze and design slopes, there is a high demand to find either critical slip surface or optimum configuration [33] [34] [35]. where ds is the differential length of a given slip surface. If the slip surface is divided into n segments with discretized length j s ∆ , the evaluation of ( ) i FOS is described as: The combination of NEM-GA has been used to calculate FOS in slope stability analysis. In the proposed framework, GA is used to find the minimum FOS for circular slip surfaces while NEM is used to analyze the shear stress and strength on the given slip surface [37]. The critical slip surface which is defined by ( )  Table  2.
The results show FOS varies from 0.98 (Janbu simplified) to 1.09 (NEM-GA). In order to compare the results from classical limit equilibrium methods, FEM and NEM-GA, slip surface for above mentioned analysis is depicted in Figure 7. The proposed slip surface by FEM and NEM-GA are coarsely identical whereas the limit equilibrium methods find the critical slip surface in different coordinates. The difference in the results could be due to the failure criteria that have been introduced in different methods. Further details will be developed and explained in Discussions and recommendations.
The discussion in the previous methods is limited to circular slip surface while the heuristic algorithms can be expanded for finding the location of critical non-circular failure surfaces [39]. Zolfaghari et al. [40] introduced a simple genetic algorithm search to determine non-circular failure surface in slope stability analysis. They used Morgenstern-Price technique [41] to analysis the slopes. Despite of conventional studies, they considered the effect of internal-slice forces of adjacent slopes. This assumption significantly increases the computational efficiency in finding non-circular failure surfaces. This study suggested a robust algorithm to solve problems including finite and infinite slopes with heterogeneous materials.
In slope stability problems with heterogeneous geology the difference in the results between noncircular and circular slip surfaces is considerable. Figure 8 Figure 7. Critical slip surface proposed by classical limit equilibrium methods, FEM, and NEM-GA.
depicts the stability analysis in a finite slope including four different soil layers.
Zolfaghari et al. [40] showed that the FOS in the circular slip surface analysis based circular slip surface using Bishop and Morgenstern methods are 1.475 and 1.5 respectively. However, the analysis with noncircular failure surface shows that FOS is 1.24. Moreover, it shows that failure surface mostly includes the weaker layers which are expected in real problems. They extended this idea in infinite slope stability and Figure 9 represents the results regarding an infinite slope stability analysis. The problem includes three infinite layers with different material properties. In this problem the seismic load is applied to the problem in pseudo-static conditions while the coefficient of quake is 0.1 which is applied horizontally. Regarding the proposed framework, the FOS is different in static and seismic conditions from 1.14 to 0.95, respectively. However, the location of slip surface is identical under both conditions.
In the abovementioned study the control variables were defined within static bounds. However, Cheng [42] introduced the dynamic bounds which are controlled by the kinematic requirement of failure. This idea later followed by  Cheng et al. [43]. They introduced a Modified PSO procedure in coupling with slip surface generator introduced by Greco [44], Malkawi et al. [45], and Cheng [42]. The MPSO procedure allows obtaining the optimized solution with fewer evaluations in comparison with standard PSO. Since particles with better objec-tive function evaluation are permitted to search (fly) more within the given iteration step, subsequently MPSO suggested a faster approach in solution procedure. Moreover, in MPSO, the number of flies within the whole group of particles is limited. Figure 10 briefly depicts the MPSO algorithm by adding few simple steps to standard PSO algorithm.
In order to compare the results of presented methods, a natural slope with four soil layer that has been studied by Zolfaghari [40] is considered which implemented Spencer method within genetic algorithm to analyze the proposed slope under 4 loading conditions: Case I: No water pressure and no earthquake loading in the system.
Case II: Water pressure and no earthquake loading.
Case III: No water pressure but earthquake loading.
Case IV: Water pressure and earthquake loading.
Cheng et al. [43] used the same number of slices in order to establish a fair comparison between GA, PSO, and MPSO. In all cases the minimum FOS is reported which is considerably lower than GA solution (reported by Zolfaghari et al. [40]). Moreover, in all cases, the portion of slip surface which passes through Figure 10. Proposed MPSO algorithm for non-linear optimization problems (Cheng, 2007).
the weakest layer is more than the solution by Zolfaghari et al. [40]. Table 3 shows the summary of the results from GA, PSO, MPSO under 4 loading cases.
For illustration intentions, the slope stability analysis results are shown for Case I in Figure 11. It is shown that MPSO suggests more critical conditions in comparison with standard PSO and GA. In order to compare the results for other cases (II, III, IV), interested readers are referred to Cheng et al. [43].

GA and PSO in Foundation Design
The last part of this study introduces a brief overview on application of evolutionary optimization techniques (GA and PSO) in foundation design. Recently, in the concepts of foundation design a new idea which is called Robust Geotechnical Design (RGD) is the matter of interest [46] [47] [48]. RGD is interested for obtaining a certain level of design robustness in addition to safety satisfaction and cost requirements. Subsequently, a single best design does not exist anymore and a trade-off among multi objectives may be required. In these conditions, Juang et al. [49] implemented GA in a multi objective optimization program to  Figure 11. Critical slip surfaces obtained from MPSO, PSO, and GA. find a Pareto Frontier which describes a trade-off between cost and robustness at a given safety level. Following Juang et al. [49], reliability-based RGD approach using GA is shown in six steps: Step 1: Classify the design parameters and noise factors in the system (foundation). The dimension of system (B: width; L: Length; and D: Depth) are the design parameters. Subsequently, the design space is finite including M finite designs. On the other hand, uncertain soil parameters ( φ : Friction angle; c : cohesion) could be remarked as noise factors.
Step 2: This step includes an inner loop ( Figure 12) which is used to compute variation of system response. In addition, coefficient of variation (COV) of noise Figure 12. Flowchart representing reliability-based RGD approach using GA. factors for each design is considered in this step. Note: The mean value of each noise is fixed while its COV is allocated based on point estimate method.
Step 3: It includes the reliability analysis for each design in the design space. The deterministic model for ultimate limit state (ULS) and serviceability limit state (SLS) for the system can be counted for probabilistic analysis. ULS and SLS probability may be calculated based on first order reliability method (FORM [48]).
Step 4: This step accomplishes the reliability analysis for N times and uses the results to calculate the mean and standard deviation of probability of failure.
Step 5: The outer loop in Figure 12 that using point estimation method (PEM [48]) for M times is accomplished in this step. Subsequently, the mean and standard deviation of the failure probability are obtained for all designs.
Step 6: The best design is determined in this step. Regarding the RGD, the multi-objective functions are defined as cost and robustness whereas safety is defined as a constraint. The cost of each design may be computed based on procedure introduced by Wang and Kulhawy [50] and the multi-objective optimization is done using non-dominated sorting genetic algorithm version II (NSGA-II) introduced by Pratp et al. [51].
RGD is based on satisfying three major factors: Safety, cost, and robustness.
Therefore, one can cast an optimization problem including robustness and cost as two objective functions which are subjected to safety as a constraint [49]. In order to find the Pareto Frontier, first we need to generate a random "parent population" in the design space which is similar to step 0 in standard GA (See Section 2). Here, we consider the number of initial "optimal" designs as n in the design space. Subsequently, a series of GA steps including selection, crossover, and mutation are executed to obtain the new "offspring population" (Step 1 in Section 2). Then the refinement of parent population is accomplished through an iterative process (Step 2 in Section 2). For a given generation, the new population is formed through the combination of parent and offspring population.
Next, the non-dominant sorting is applied over new population and produces different sets of design points based on the non-dominant ranking. For instance, the best class is recognized as F 1 ; the next is F 2 and so on. Subsequently, the best n points are nominated for next generation. Crowding distance technique can be used in the generation on new population with respect to keep diversity in the selected design points when the number of points in 1 2 , F F  exceed n [51].
This procedure continues until the convergence criterion is satisfied and stable solution obtained. Subsequently, the final solution in the proposed optimization process is determined as Pareto Frontier. Figure 13 shows the converged solution (Pareto Frontier) for a spread footing using NSGA-II optimization.
Pareto Frontier includes an optimal set of possible designs in RGD approach.
This set plays a trade-off role in relationship between cost and robustness against uncertainty factors. The designer may opt for a greater robustness at a higher cost or vice versa. On the other hand, if a threshold is predefined for cost then the designer could decide on the most robust design from the Pareto Frontier. In Figure 13. The Pareto frontier for the design of spread footing obtained using NSGA-II.
addition, based on feasibility robustness one can refine the decision for design objectives ( Figure 14). For further information regarding feasibility robustness, interested readers are encouraged to study was done by Juang et al. [49].
Application of evolutionary methods in foundation design is not limited just to spread footings. For instance, Lei [52] implemented PSO on the settlement fitting prediction of highway foundations. The settlement of highway foundations plays a major role in the deformation and stability of the roads. In general, there are three approaches to calculate settlement fittings: 1) Implementation of soil constitutive models to calculate the final settlement of road bases, this type of analysis is based on consolidation theories [53].
3) Employment of direct methods which are based on actual measured data (settlement) and presenting a mathematical model to predict settlement in road bases [55] [56] [57].
Lei et al. [52] used the direct method which is based PSO and time series analysis to present a subgrade settlement model for Nanjing-Hangzhou highway. Since settlement in foundation follows a time dependent trend, one can represent the settlement via a polynomial based time series. Time series analysis has been successfully applied in dynamic data processing phenomena (i.e. Prediction of dynamic deformation of landslides [58]: If m and jk a are determined with PSO, the above model turned to a predictive model which can be used to predict the future deformation prediction [52].
In the proposed study, parameters in the model are generally determined by empirical methods. In the process of finding fitting parameters for the model, the measured data has a considerable effect on the predictive model. Lei et al. [52] Regarding the predictive model presented in Equation [10], the learning value, predictive values are compared with the measured data in Figure 15. The comparison shows a good agreement in the analysis of measured data and predicted values.

Discussions and Recommendations
Seepage through earth dams or ground water flow represents a problem with an unknown boundary. In order to find a reasonable solution, this study suggests the simultaneous application of a robust technique for solving the governing Slope stability analysis is a challenging problem in geotechnical engineering and it needs the application of constitutive models which include elasto-plasticity [60]. In NEM-GA framework only elastic behavior of material is considered which can result in error in slope stability analysis. However, it is encouraged to extend the proposed framework to elasto-plastic materials. Subsequently, the application of evolutionary methods in the optimum design of deep foundations and soil reinforcements is strongly suggested for future studies.

Conclusions
In the presented review, we briefly introduced two well-known evolutionary methods (GA/PSO) and their application in geotechnical engineering. The implementation of GA/PSO in three major geotechnical problems is considered: 1) Unconfined seepage problems; 2) Slope stability analysis; and 3) Spread footing designs.
Regarding the unconfined seepage problems, at least one of the boundaries In slope stability analysis, determination of the failure surface location is a challenging problem in geotechnical engineering. In this study, we found that GA, PSO, and Modified PSO are successful in finding critical slip surface and subsequently estimation of FOS. In real life scenario, usually critical slip surface is not a semi-circular, especially in problems including heterogeneous soil layers.
In the presented study, MPSO performs adequately in the determination of critical slip surface with respect to number of iterations and computational cost.
In the last part of study, the application of GA/PSO in foundation design is investigated. It is shown that the combination of GA and RGD leads to find the trade-off between cost, safety, and design robustness. This criterion is depicted by a Pareto Frontier. Finally, the application of PSO in prediction of settlement in shallow foundations is presented by studying the spread footings in highways.
It depicts that PSO is successful to predict the settlement trend based on observed data.
The safety in designs is highly respected in geotechnical engineering. On the other hand, contractors and private sectors prefer economical design with minimum labor effort and construction cost. Subsequently, implementation of robust optimization techniques which establish a trade-off between safety and total cost of geotechnical projects is necessary in practice.