Locating Critical Circular and Unconstrained Failure Surface in Slope Stability Analysis with Tailored Genetic Algorithm

Abstract This article presents an efficient search method for representative circular and unconstrained slip surfaces with the use of the tailored genetic algorithm. Searches for unconstrained slip planes with rigid equilibrium methods are yet uncommon in engineering practice, and little publications regarding truly free slip planes exist. The proposed method presents an effective procedure being the result of the right combination of initial population type, selection, crossover and mutation method. The procedure needs little computational effort to find the optimum, unconstrained slip plane. The methodology described in this paper is implemented using Mathematica. The implementation, along with further explanations, is fully presented so the results can be reproduced. Sample slope stability calculations are performed for four cases, along with a detailed result interpretation. Two cases are compared with analyses described in earlier publications. The remaining two are practical cases of slope stability analyses of dikes in Netherlands. These four cases show the benefits of analyzing slope stability with a rigid equilibrium method combined with a genetic algorithm. The paper concludes by describing possibilities and limitations of using the genetic algorithm in the context of the slope stability problem.


INTRODUCTION
Slope stability analysis remains to this day one of the fundamental problems in geotechnical engineering.The most commonly used computational methods are the limit equilibrium method (LEM) and the finite element method (FEM) as described in [12], [14], [15].The LEM with its method of slices is used by engineers to assess the stability of slopes because of its reliability, speed and widespread experience with the method.The determination of the potential slip surface becomes particularly important in the vicinity of infrastructure, development housing, or when an embankment protects flood plains.The focus in the majority of cases is on determining the geotechnical parameters for calculations [18], [21] but, as studies presented in this article show, the right location of potential slip curves significantly affects the factor of safety (Fs).
The problem of searching for the critical slip curve is a typical optimisation task, consisting in finding the minimum of the Fs function, depending on geometrical data of the slip curve: {X 1 , Y 1 , R 1 } (Fig. 7) for circular and {{X 1 , Y 1 }, ..., {X n , Y n }} for unconstrained slip surface (Fig. 8).Publication [19] presents the formulation of a function of the factor of safety Fs according to Bishop, depending on variables {X, Y, R}, what resembles the basis of function stability formulation and the method of slices.In addition, in monograph [10] formulas and derivations as well as the origins of the development of method of slices have been presented in detail.
Nowadays, most computer programmes used for slope stability analysis on the market exploit a method for searching for a critical circular slip curve which involves sequential calculation of all combinations of grid points and tangent lines.The next step is to report the lowest registered value.In the context of developing methods of heuristic algorithms used in geotechnical engineering [3], [4], [5], [6], [7], [9], [13], [19], [23], [25], [26], this approach is not only time consuming but also imprecise and virtually impossible in the case of an unconstrained slip surface.Moreover heuristic algorithms are particularly useful in the case of evaluation slope stability with variational method [20].
The genetic algorithm is one of many heuristic algorithms used for optimisation tasks.In monograph [8], the basis of the algorithm is presented in a comprehensive manner and is explained using examples created in the Turbo Pascal software.In the Mathematica software [24], the GA procedure, based on the example of finding the minimum of the function of one variable with description of specific sub-programmes, is presented by Bengtsson [1].The pseudo-code presented in publication [11] provides a good illustration of the GA.The GA procedure for specific optimisation tasks presented in this article can be seen in more detail in Figs. 4 and 9.As the effectiveness of the genetic algorithm cannot be proven in purely theoretical manner, its speed and accuracy is best tested experimentally, as in this article.
The researches considering implementation GA to slope stability analysis have been carried out by many authors.The most relevant papers are: Zolfaghari et al. [26], Sengupta et al. [19], Li et al. [9] and Zhu et al. [25].Zolfaghari et al. [26] presented a simple genetic algorithm combined with Morgenstern-Prices for noncircular and Bishop for circular slip surfaces.The authors presented an advantage of GA over Simplex method.Sengupta et al. [19] proved that GA with Bishop method used less CPU memory than Monte-Carlo or grid method.Li et al. [9] showed an interesting real-coded genetic algorithm for unconstrained slip surface.In their work, comparison of the number of slip surfaces needed to obtain the minimum factor of safety Fs to other authors is made.Zhu et al. [25] introduced a study of hybridization GA with Tabu Search Algorithm (HGA) in order to fill GA deficiency of local searching.Authors presented examples where Hybrid GA searching procedure for circular and noncircular slip surfaces was applied.
Another contribution was made by Cheng [3] who implemented Modified Particle Swarm Optimization (MPSO) algorithm to Spencer method.
In the present article, a tailored genetic algorithm proved to be more effective and faster because of mitigation of the number of computed slip surfaces.In HGA or MPSO more slip surfaces need to be considered to get similar value of Fs: 14 000 for MPSO and 15 000 to even over 150 000 for HGA for circular and unconstrained slip surfaces respectively.The algo-rithm presented in this work needs 2500 slip surfaces at most to compute reliable factor of safety Fs (Table 1).This is due to the tailored properties of the algorithm, designed for an unconstrained slip surface search.In practical use the resolution of two decimal places of Fs is sufficient.
This paper shows the advantages of using a genetic algorithm combined with a rigid LEM with a truly unconstrained slip plane search by presenting a relatively simple GA to optimize an analytic solution and a simple slip circle (chapter two).The complexity in this chapter increases until a representative unconstrained slip plane can be found.After presentation of the tailored GA, chapter 3 is dedicated to four examples that show the effectiveness of the method.The GA is tailored in such manner that it finds the optimum free slip plane with relatively little computational effort (Table 1).
The presented method is relatively precise and very quick in determining the safety factor.It gives reliable answers under such conditions where Bishops method is not applicable, like layered soils and uplift conditions.The method will calculate the safety with the precision of a FEM method and the ease of a LEM method.The paper aims to give insight in the genetic algorithm and prove, using different examples, that this is a robust and reliable method.

MINIMUM OF RASTRIGIN'S FUNCTION
A standard task portraying the quality of a search algorithm is to find the minimum of the Rastragin's function f (X, Y, R) (Fig. 1).The set of solutions is located between the upper and the lower plane.Visually, it is easy to see that the minimum is f (0, 0, 0) = 0.In a computational procedure, it is difficult to find the optimum, though, as this function has many local maxima and minima.Figure 4 shows the procedure of optimization using the genetic algorithm for both the analytical tests and the Bishop analysis.
The red numbers in this figure correspond to the chapters of the programme placed in the supplementary data [16].Presenting the entire code allows for easy tracing of the computational process of the programme.
Locating critical circular and unconstrained failure surface in slope stability analysis with tailored genetic algorithm 89 Calculations are conducted for the following set of data: 1. pop = 50 -number of population; 2. gen = 100 -number of generations; 3. nrofLeaders = 2; 4. mutationRate = 0.2; 5. aa = 0.6 -parameter of roulette wheel selection; In the initial population, the elite is the value f (0.183295, -0.226707, 6.37283) = 7.56758 corresponding to the red point marked in Fig. 2 after 100 generations and the results f (0.00131752, -0.0148293, 0.00622) = 0.00943332 are obtained corresponding to the red point marked in Fig. 3.
As can be seen, the final result differs slightly from the exact solution f (0, 0, 0) = 0.The inconvenience of the algorithm lies in the fact that there is no clear stopping criterion.In addition, the algorithm is a chaotic process which each time gives a different but similar result for the same set of parameters (Fig. 5).The presented GA finds the position of global minimum quickly and rather precisely with a very small risk of getting stuck in a local minimum.These are the main traits of this algorithm compared to other methods of function optimization.The effectiveness of a GA for optimization tasks can be traced by ob- serving the behaviour of population of different generations in time, using the intermediate results.Figure 6 presents the average of the population fitness and the matching of all the individuals (strings) to the target function.This type of graph can be applied to assess the parameters used and to optimize the whole procedure to obtain a more accurate result in shorter time of computation.

MINIMUM FS FOR THE CIRCULAR FAILURE SURFACE
The circular slip surface is a simple shape used for analysing slope stability with the use of LEM, especially the Bishop method.The formulation of the genome is rather easy.An circle is defined using three variables {X, Y, R}.These variables correspond to the centre point and the radius of the slip circle (Fig. 7).The same genetic operators as in the case of Rastrigin's function are used (Fig. 4).The main difference between the two functions is the way the results are presented.The solution space of the Rastrigin's function is easily shown as it is a simple function.The solution space of the slope stability analysis cannot be constructed quickly in an analogous manner, since it would require calculating the factor of safety Fs for all points {X, Y, R} of the area considered, which would take an enormous computational effort.The wheel is scaled with parameter "aa".If its value equals 0, the selection is completely random, and if its value equals 1, the chance of selection is based upon the relative fitness of the individual (parent); 3. Creating a child through crossover crosses randomly at 1 or 2 points on the genome; 4. In each generation a defined number of strongest members, the so called elite "nrofLeaders -nL", is directly passed through to the next generation; 5.The evaluation of the fitness of individuals is conducted using the D-Geo Stability software, after having exported data (the defined slip surfaces) from Mathematica.The calculated safety factor is imported from the D-Geo Stability output (see Figs. 4 and 9 for cells 3 and 4, respectively).It has been decided to use an existing programme to calculate slope stability, because writing a new equilibrium method from scratch requires an enormous amount of time and is not the subject of the research conducted; 6. Mutation causes a random gene with a given probability "mutationRate -mR" to change to a random value within a given range.

UNCONSTRAINED FAILURE SURFACE
The formulation of a free slip surface is more complicated than for a circular surface.The genome becomes more complex and breeding strategies become more important.The number of transversal lines "ln" is defined (Fig. 8).The first transversal line must follow the surface line where the slip plane enters the soil body.The last line must follow the surface where the slip plane exits the soil body.For each transversal line we adopt a length equal to 1.This length is measured from the slope into the ground (Fig. 8).In this case, the individual defining a random slip surface consists of an array of value ranged from 0 to 1.The genome is defined in the following manner: {l 1 , l 2 , l 3 , l 4 , l 5 , l 6 , l 7 , l 8 , l 9 }.Knowing the positions of the ends of "l n ", one can easily match them with the values of coordinates {x n , y n } corresponding to a given value "l n " (red line in Fig. 8).The number of possibilities of defining the individuals for any free-slip surface is practically infinite.However, several factors, decisive for the effectiveness of the algorithm, have to be taken into account, including: 1. Simplicity; 2. Freedom in creating a shape of a free-slip surface -the fewer constraints the better; 3. Generated slip surfaces have to be realistic in majority.If there are too many unrealistic slip surfaces, the procedure becomes relatively ineffective; 4.This definition has a constant genome length even if the size of the slip plane increases or decreases.This keeps the breeding process more simple.Figure 9 presents the flowchart of the genetic algorithm to find the representative unconstrained slip surface.Each number marked in red corresponds to a chapter of the programme in Mathematica software attached in the supplementary data [16].For clarification, the following remarks have to be made by the flowchart.
The first 5 points of the flowchart for the unconstrained slip plane are identical to the flowchart that belongs to the Rastragin's function and circular slip surface.Thereafter, the following adaptations are made: 1. Generating the initial population takes place by drawing parabolas in the l n space.This causes chaotic but somewhat realistic initial slip surfaces for the initial population; 2. Three kinds of mutations of individuals are defined: jump, creep and inversion: (a) Jump consists of changing a randomly selected l n by a random value in the interval (0,1) (b) Creep involves a random increase or decrease of the value of the randomly selected l n by a given percentage "pC" (c) An inversion changes the value of all l n by 1 − l n Mutations are controlled by the mutation probability parameter "mutationRate -mR", which is further divided into parameters of specific kinds of mutation "jProb, cProb, iProb"; 3.After the genetic operators are applied on a population, individuals (genomes) are adapted to ensure Fig. 9. Flowchart explaining the genetic algorithm for locating the unconstrained failure surface.Red numbers correspond to the numbers in the supplementary data [16] Locating critical circular and unconstrained failure surface in slope stability analysis with tailored genetic algorithm 93 that the value of variable X of the subsequent point increases, relatively to the previous ones, thus eliminating unrealistic slip planes.[3], [9], [26].In these two examples, the algorithm described in this paper is compared in terms of precision and efficiency with the examples from presented in literature.The two subsequent examples are practical calculations of stability of anti-flood defences in the Netherlands in complex soil and water conditions.Additionally, Table 2 presents parameters of the genetic algorithm used in the subsequent calculations, corresponding to resulting values of factor of safety Fs.

EXAMPLES 1 AND 2 (ACADEMIC)
Figures 10 and 11 show the positioning of the subsoil layers as well as the representative slip surface divided into slices.The adopted parameters of soil and the geometry of slopes conform to the corresponding publications [3], [9], [26].From the calculations conducted one can see that the genetic algorithm has similar precision to the examples from literature.A clear advantage of the obtained results, compared to those of other authors calculations, is the number of slip surfaces under consideration, which directly shortens the time used for calculations (Table 1).Figures 12 and 13 present the initial and the final population for given conditions of searching for the slip surface (blue lines).The red colour has been used to mark the critical slip surface.Examples 1 and 2 from this chapter are simple examples with one global minimum, which does not constitute a major difficulty for a search algorithm.The differences in results obtained by various authors arise from matching the slip surface to the privileged surface in between layers boundaries with weak and strong properties along them with and a smooth entrance into the surface.It is obvious that the variation in reaching the function minimum depends on the adopted parameters of the algorithm and on the stochastic process itself.

EXAMPLE 3 (REAL CASE)
Example 3 regards a cross-section of the embankment along the trajectory Kinderdijk Schoonhovenseveer.This cross section is chosen as an example, because there is a vast amount of data readily avail-   Figure 16 presents the last population in the search area for the unconstrained slip surface.The critical curve resembles a circular one.This makes sense, as the water pressures in this area of the cross section are hydrostatic and the soil is relatively homogenous.As a circle is the representative failure mode, a Bishop analysis with three different GA input parameters was conducted for comparison (Fig. 17). Figure 18 shows the change of the factor pop -population count, gen -number of generations, aa -parameter of scaled roulette wheel selection, nL -elite size, mR -mutation rate, pC -percentage move in "Creep" mutation, cProb -probability of "Creep", jProb -probability of "Jump", iProb -probability of "Inversion".

EXAMPLE 4 (REAL CASE)
Example 4 was originally set up to test a regional levee along the "Gravelandse Boezem".This example was chosen because this cross section has two possible slip planes with almost the same safety factor.In a search algorithm, it is interesting to see how the procedure behaves when there are two local minima with almost the same safety factor.Figure 19 shows the representative slip plane, but the failure mode where the berm fails has only a slightly higher safety factor.
Figure 20 shows transversal lines defining the search area for the slip surface for Example 4b as well as the population of the last generation with the elite (red colour).Observing the population for each of the generations, it can be noticed that two critical slip surface exists.A high, circular one, and a long and deep slip plane.Clearly, there are two local minima (Table 1, Example 4c, d).The deep slip plane is caused by the effect known as Uplift [22].This effect is caused by high pressure of the ground water in the permeable sand layer, caused by a high water level in the neighbouring river.As a result, the normal effective stress under the foundation of the dike can decrease even to zero, thus causing instability of the slope.
The analysis shows, though, that the higher circular slip plane is representative as it has the lowest factor of safety.This slip surface resembles a circular one.Therefore, further calculations can be conducted for a circular slip surface using Bishops method (Table 1, Example 4c, d).

DISCUSSION
In the supplementary data [16], the code of the algorithms used in Mathematica software is presented, which can easily be translated into any other programming language.This provides a good starting point for discussions and comments on the applied solutions and gives room for improvement.Flowcharts, which contain numbers corresponding to specific chapters and definitions in the presented code, make it easier to understand the algorithm.
The effectiveness of the GA depends to a great extent on the choice of its individual components, namely the generation of the initial population, selection, crossover and mutation methods, etc.That is why the genetic algorithm programmed by another author will always be unique.The correctness of the solutions used can be proved experimentally, the effectiveness depends upon the choice of the right components.
Other important factors determining the result include the definition of the search area and the parameters assumed (Table 2).It is important to have an "engineers intuition" supported by earlier calculations and evaluations of soil and water conditions.This enables understanding of the influence of the choices on the effectiveness of the GA in the particular application.
Using a genetic algorithm for optimization will always remain a challenge.A resulting slip plane will never be the exact solution, but one that lies closely to the exact solution.Based on the genetic algorithm one will not know how far the solution obtained is from the exact one.Pictures in Figs. 2 and 3 in chapter 2 are good illustrations of this issue.In the context of slope stability, the critical surface of the considered function cannot be seen and, consequently, there is no visual or numerical confirmation of how much the current solution differs from the exact one.This problem can be addressed by using a larger number of considered slip planes or by enlarging the population or the number of generations.Along with an increasing number of considered slip surfaces, computational time also increases.In this case, techniques such as parallel computing [20] are required in order to decrease time of computations.In author's opinion a few seconds with regular PC make algorithm useful enough.
Also, visual judgment of the slip plane in terms of its "smoothness" and its thrust-line while passing along potential privileged slip surfaces is valuable and, if interpreted correctly, will lead to better results of the factor of safety Fs.
Another way to obtain more accurate results is to make new calculations within a smaller search area, limited to the area of previously found critical curve.Then, the precision of the safety factor of the critical slip surface increases because the algorithm limits itself to one minimum (Example 3f).
A fundamentally different approach to reach a higher precision is to end the genetic algorithm with a hill climbing type of optimization.This will smoothen the slip plane and reduce the safety factor found.The advantage of using such a hybrid method is current subject of investigation.

CONCLUSIONS
The possibility to observe individual populations during generation is an undisputed advantage of the presented approach.This way other local minima can be noticed, which is of importance and can provide information about the location of equally dangerous slip surfaces (Example 4).It also aids in choosing the right parameters for the genetic algorithm to perform optimally.
The advantage of presented work, compared to previous ones, is its effectiveness in context of considered number of slip surfaces (small numbers of fitness function evaluation) and simplicity in creating the unconstrained slip surfaces.This effectiveness is reached through the right combination of initial population type, selection, crossover and mutation method.
Presented researches confirm that proposed solutions need small numbers of computed slip surfaces to get optimal value of FS (Table 1).

Fig. 4 .Fig. 5 .Fig. 6 .
Fig. 4. Flowchart of genetic algorithm for analytical test and locating critical circular failure surface.Red numbers correspond to the numbers in the supplementary data[16]

Fig. 7 .
Fig. 7. Formulation of task for circular slip surface To optimize Rastrigin's function and circular slip surface the following assumptions are made (based on the graphs presented in this chapter and on prior experience with genetic algorithms): 1.The genome is represented with real numbers for X, Y, R; 2. Scaled roulette wheel is used for selection of parents.The wheel is scaled with parameter "aa".If its value equals 0, the selection is completely random, and if its

Figure 14
Figure14presents the change in the minimum factor of safety Fs for each generation for all the calculations conducted in Example 1.It is obvious that the variation in reaching the function minimum depends on the adopted parameters of the algorithm and on the stochastic process itself.

Fig. 14 .
Fig. 14.Minimum value of Fs versus number of generations for Example 1a, b, c, d

Fig. 17 .
Fig. 17.Last population for Example 3f Fig. 18.Minimum value of Fs versus number of generations for Example 3d, e, f

Table 1
presents the results of calculations for the 4 examples.The first two examples are presented in literature

Table 2 .
Parameters of GA applied for calculation