Ant colony optimization for slope stability analysis applied to an embankment failure in eastern India

The safety of road embankments is mainly assessed in the engineering practice by limit equilibrium analyses. Locating the critical slip surface of embankments and calculating the corresponding factor of safety is a crucial task. In this paper, the continuous ant colony optimisation algorithm is used to analyse the stability of slopes with non-circular slip surfaces. To illustrate the proposed procedure, one example from published literature and another engineering case study of a landslide at a road construction site in India is analysed. This latter study is supplemented by the results of geotechnical investigations performed before and after the failure of the embankment. The results demonstrate that the proposed technique identifies correctly the critical slip surface and can thus be used for engineering applications.

Cheng et al. [30] compared six metaheuristic optimisation methods and reported that no single method can outperform other methods in all aspects, each one having its own merits and limitations. In all these methods, it is paramount first to identify the critical slip surface. This critical slip surface is found by comparing the factors of safety obtained from several feasible trial surfaces and choosing the one with the minimum value. This may sometimes be erroneous due to the presence of local minima points. Many existing methods fail to identify the global minimum and may converge at a local minima even for simple homogenous slopes [31]. Neural networks on the other hand are site specific, require reliable training data and apply for particular cases only [25,32,33]. Discrete ACO used previously for FS calculation takes too many fitness function evaluations to compute optimum FS and has more uncertainty in FS estimation. Therefore several global optimisation algorithms are available to locate critical slip surfaces but are either too complex to implement for practitioners based on an excessive number of fitness function evaluations. The use of finite element techniques bypasses the need to assume the shape and location of the trial slip surfaces, but also introduces other complexities including defining appropriate soil constitutive models.
Most of the existing software packages available for slope stability calculations [34-36] implement a method of slip surface search called point grid approach. The location of the centre point of the slip surface, that in this case is assumed as circular, is shifted during the search until it corresponds to the lowest factor of safety. However, this method does not assure that the slip surface retrieved corresponds the global minimum. Some programs [37] implements cuckoo search [38] and simulated annealing [39] for noncircular failure surfaces. Monte-Carlo search techniques are also implemented in limit equilibrium slope stability programs (SAS-MCT 4.0) [40]. For evolutionary computation, codes such as Slope SGA [29,41] incorporate a Genetic Algorithm (GA) for both circular and non-circular slip surfaces. The commercially available software Geostudio [35] uses an optimising function which is based on search techniques used by Greco [16] and Malkawi et al. [40] to search for critical slip surface. Apart from that, to the authors' knowledge, there is no other commercial program using intelligent algorithms optimisation techniques for slope stability analysis. Published research [17,18,30] demonstrates that there is a growing interest in ant-based algorithm metaheuristic for slope stability assessment.
Most of the intelligent algorithms for slope stability calculations are relatively new, difficult to implement and not familiar to geotechnical engineers. Sengupta et al. [27] used a circular surface with the genetic algorithm to find the location of the critical slip surface and the corresponding factor of safety using Bishop's method of slices [1]. The factor of safety was defined as a fitness function expressed in terms of coordinates of the centre of the slip surface and its radius (three control variables). For a non-circular slip surface, the control variables will be more than three. Kahatadeniya et al. [17] used discrete ant colony optimisation and Kashani et al. used [42] Imperialistic Competitive Algorithm (ICA) to determine the optimal curve yielding minimum factor of safety using Morgenstern-Price method [3]. The ant colony optimisation technique used by Kahatadeniya et al. [17] was in discrete search space and required careful calibration of parameters to arrive to correct results. Furthermore by using discrete ACO the uncertainty in the estimation of FS was large and many fitness function evaluations were needed to arrive to an optimum value. In ACO, the number of fitness evaluations equals to product of population size with the total number of iterations as in each iteration the fitness function is evaluated only once. As different algorithms use different amounts of consumption of fitness function evaluations per iteration so a comparison made only on maximum number of iterations as a stop condition is not justified [43,44].
In this study, an extended form of ACO, known as continuous ant colony optimisation ( ACO R ), capable to solve problems in the continuous domain, is used for computing the factor of safety for slope stability calculations. The minimum factor of safety is identified using ACO R with control parameters varying continuously, differently than previously used approaches of ACO working with discrete intervals. Furthermore, the ACO R is validated against theoretical benchmark example. Additionally, an actual field failure is back-analysed using this technique. We will show that results obtained by the proposed technique are in agreement with those of other researches and with those obtained by standard software.
Apart from being applied in slope stability, ACO and its variants have been used in other civil engineering optimisation problems, e.g. identifying structural damages from changes in natural frequencies [45][46][47], allocation of railway platforms [48], design of concrete retaining walls [49], damage detection in structures [50], guiding vehicles to less congested paths [51,52], forecasting river flow [53] to mention a few of them. These works demonstrate the potential of ACO techniques in solving civil engineering optimisation problems.

Setting of slope stability analysis
For setting up a slope stability analysis, the first step is to generate a kinematically feasible slip surface for the slope and calculate the factor of safety associated with it. The study considered a non-circular slip surface by linking straight lines using a method proposed by Cheng [54]. To start the generation of a slip surface, the failure mass is divided into N slices of equal width as shown in Fig. 1. The limits of variable x 1 and x N +1 are the entry ( x l , x u ) and exit points range ( x L , x U ) ( Fig. 1) which depend on the experience of the practitioner while all other control variables ( σ 1 ,…, σ N −1 ) vary between 0 and 1. The slip surface can hence be described as follows: In terms of coordinates it can be written as: The x co-ordinates of points lying between entry and exit ( x 2 ,…,x N ) are equally spaced as follows: For determining the y-coordinates of each vertex ( V 2 ,…,V N ) the upper and lower bounds ( y imax and y imin ) are considered, where i denotes the index of vertices. The process of determination of upper and lower bounds of y i is shown in Fig. 1. As can be seen y-coordinates of each vertex are different with each slip trial surface, so control variable σ is used which always varies between 0 to 1 irrespective of the considered slip surface. Hence, y i can be determined dynamically between the maximum and minimum values for y i by the following equation: The optimisation problem for searching a critical slip surface can be summarised as: where V can be obtained by the procedure explained above.
A generalised slip surface is thus defined by a set of N + 1 decision variables to optimise; x 1 , σ 1 ,…,σ N −1 , x N +1 . Each set of combinations of decision variables represents one trial slip surface.
After the generation of a trial slip surface, the second step is to determine the factor of safety (fitness function) associated with each slip surface. The original formulation of Morgenstern-Price method for calculating FS of a slope involves nonlinear equations which have to be solved numerically. The algorithm developed by Zhu et al. [55] based on the Morgenstern-Price method is used here to calculate the factor of safety. Finally, all the control variables [ the termination criteria is met using the ACO R algorithm to find the right combination leading to the minimum safety factor.

Ant colony optimisation (ACO)
Ant colony optimisation proposed by Dorigo and Stutzle [56] is a biologically inspired swarm optimisation method mimicking the foraging behaviour of ant colonies. They find the shortest possible path from their origin point (start) to the food source (finish) by releasing a chemical known as pheromone. When a lone ant randomly tries to find the food source, it can sense the pheromone trail left by other ants that have previously followed the same path. The pheromone trail stores local information of the path and is reinforced as each successive ant passes. It is a population based search technique for optimisation of problems having many inputs based on trail-laying and trail-following of its members to find the optimum solution. This foraging behaviour can be artificially simulated and can be used to solve engineering problems.
The research presented here explores the ACO optimisation algorithm in calculating the minimum factor of safety of the slope. In fact, the problem of finding the critical failure surface is very much similar to problems leading to a shortest optimised path. Figure 2 shows the permitted paths for a number of variables with their limits. The objective is to find the path which leads to a minimum factor of safety. The weighted connected graph shows the set of nodes and connections between them. In each case, a node represents the value of the variable to be chosen. The ants move from node to node to the right recording each node visited so as to update the pheromone trail. The lines between the nodes indicate the flow of ant path from one node to the next, chosen probabilistically as explained in the following point 1 on probabilistic transition and depending on the intensity of pheromone on that node. When ACO framework is implemented as an algorithm for solving either discrete or continuous optimisation problem, the following rules have to be followed by ants: 1. Probabilistic transition: Fig. 2 Architecture of the discrete ACO algorithm and construction of search space for m number of ants for optimisation of slope stability analysis ( x 1 , σ 1 , …,σ N−1 , x N+1 are control variables with their limits for FS computation as shown in Fig. 1) where p k ij is transition probability for ant k to choose node j when located at node i, τ i,j is the amount of pheromone on trail associated with node (i, j), N k i is the possible neighbourhood nodes of ant k when placed at node i and α > 0 is a parameter denoting the degree of importance of the pheromones. α is generally chosen as 1 [56] to allow diversification in solutions otherwise a higher value of α will influence ACO to intensify amount of pheromone to previously visited nodes leading to premature convergence.
The rule as shown in Fig. 3 explains how the ant will select the path with a higher level of pheromone concentration. The ant path is characterised by vector 5 as this path has maximum pheromone concentration at these nodes (marked by solid circles). In essence, calculation in this case is carried out starting from the first control variable ( x 1 ) and exhaustively identifying all possible paths until reaching the finish variable ( x N +1 ). The movement to successor node is determined by the stochastic decision rule given in Eq. 6 which is function of its local pheromone level and relative importance of pheromone.

Updating pheromone concentration:
where ρ is the evaporation rate varying from 0 to 1 of pheromone and controls the quantity of pheromone on each trail path and �τ k ij is the additional pheromone Fig. 3 Selection process of possible ant paths as a function of pheromone concentration level ( τ ij ). The solid nodes representing high pheromone concentration will have more probability of selection than their counterparts deposited by kth ant at node (i, j). τ up ij and τ prev ij denotes the updated and previous pheromone trail at node (i, j). A higher value of ρ will reduce the number of significant iterations, this decreasing the diversification of solutions leading to premature convergence. For discrete ant colony optimisation algorithm, a number of artificial ants are considered (m). Each ant builds a solution by visiting node to node on the graph (Fig. 3) in forward direction where each node has transition probability attached to it. An ant visits nodes representing variables of the optimisation problem (1 to N + 1). The initial intensity of pheromone on each node can be initialised to 1, so that all the paths have an equal probability of being chosen and no path is being left out. Hence, at first iteration, all the nodes have the same pheromone concentration that is τ ij = τ 0 = 1 (i = 1, 2, …, N + 1; j = 1, 2, …, m). The pheromone concentration at each node is globally updated dynamically in each iteration after all ants have successfully deposited the pheromone on their trail path. The pheromone updating aims to boost the pheromone trails in promising regions. The shorter path will have a higher quantity of pheromone level than the longer path owing to its lower evaporation. The iterations are carried out for a number of cycles to produce the optimum solution. At the end of iteration cycle, based on the quality of solution constructed by ants, more pheromone will be concentrated on nodes that construct the best solution.

ACO for continuous optimisation-ACO R
An extension of the original ACO algorithm developed for discrete domains is continuous ant colony optimisation ( ACO R ) which works for decision variables in continuous domains. The ACO R has been more practical to use as in practical applications decision variables lay in continuous domains. The continuous ACO algorithm was proposed by Socha et al. [57] inspired from discrete ACO and works without making any major changes in the original ACO framework. The fundamental idea of ACO R is to consider continuous probability density function (PDF) P(x) instead of discrete PDF according to Eq. 6 to construct solutions. As the domain of each variable in ACO R is continuous, an ant sampling is based on a continuous probability-density function P(x). The procedure of using ACO R meta-heuristics is outlined below.
In this study, slip surface is represented by a set of N + 1 decision variables; x 1 , σ 1 ,…, σ N −1 , x N +1 . For a N + 1 dimensional continuous problem, an archive of k solutions denoted by S (Eq. 8) is constructed. Each solution in the archive contains the value of N + 1 variables used to construct a generic trial slip surface and the factor of safety associated with it. This solution archive is equivalent to pheromone model used by discrete ACO to keep a memory of the search history. As pheromone cannot be updated directly, the equivalent step is updating the solution archive. The size of the archive k, which is also a controlling parameter of the algorithm, must be greater than the dimensions of variables being solved. For example, in our case we are using 30 slices (N) for slip surface discretisation having 31 ( N + 1 ) control variables, so a minimum archive size of 31 must be taken. The structure of solution archive in ACO R comprises of three matrices as follows.
The position of control variables defining slip surface outlined in "Setting of slope stability analysis" section are stored in the following matrix S of size k × N + 1 : The matrix S represents k trial slip surfaces used in ACO R solution archive. For a more general case matrix S of size k × N + 1 can be written as where s i j denotes the value of jth solution of the ith unknown variable. The movement of ants along the path leads to different slip surfaces, with a set of variables stored in matrix S.
f is the matrix storing the factor of safety of each slip surface: ω is the matrix saving the associated weights with each solution of the archive: To initialise the algorithm, random sets of all k solutions having N + 1 control variables are chosen within the specified parameter range to construct all the rows of the archive matrix S, which is equivalent to pheromone initialisation in discrete ACO algorithm. The solutions in the archive are always sorted according to fitness function ( where ω is defined as matrix of associated weights of each solution whose elements are computed by employing a gaussian function [57]. For the jth row of the archive associated weight ω j is defined as follows: where q is the selection pressure varying between 0 to 1 of the algorithm which controls the exploration and exploitation rate in ACO R . If q is small, the best-ranked solutions are preferred, while for larger values of q, the probability of choosing control variables becomes more uniform. The ant subsequently chooses the value s i j and samples its neighbourhood for a new value for variable i = 1, 2,…, N, N + 1 and reiterates this operation for all control variables using the same jth solution using a probability density function (PDF). Although different choices are possible, in this work a Gaussian function is again chosen [6]: where for the mean µ i j , the values of the ith variable for all the solutions deposited the archive S become the elements of the vector µ i j : Furthermore, the values of standard deviation σ i j can be calculated dynamically as the average distance from the chosen solution s i j to other solutions and then multiplying with the parameter ξ > 0 which determines the convergence speed of the algorithm as follows: After each iteration, the set of newly generated solutions will be refreshed in the archive, and then the archive for keeping the best k solutions probabilistically. The candidate solutions stored in archive are used to modify the probability distribution such that it is biased towards sampling high quality solutions. This process is repeated until the maximum number of iterations are reached.

Application of the proposed technique
For the purpose of demonstrating the efficacy of the ACO R algorithm to locate the critical slip surface and the associated minimum safety factor, two case studies are analysed. One case study is from published literature and the second one is from a failed road embankment located in eastern India. Since ACO is a non-deterministic method i.e. giving different solution for each computation, the FS values and standard deviation reported are based on a number of simulations equal to 10 to capture the uniformity in the solutions.
As described in "Theoretical background of ant colony optimisation (ACO)" section, in the ACO R introduction, the size of the archive (k), the number of ants (m), the selection pressure (q) and pheromone evaporation rate ( ξ ) play important role in convergence speed of the algorithm. To strike a balance between diversification and convergence of the search process, proper selection of these parameters is important. After carrying out parametric analysis and following the suggestions by Socha and Dorigo [57] and Majumdar et al. [58], the constant parameters of the proposed algorithm have been set up as follows: m = 50, q = 0.5, ξ = 1, k = 40 and maximum iterations = 200.

Literature example
The literature example with four soil layers has been taken from Zolfaghari et al. [29]. Geometry of the soil slope and soil properties are shown in Fig. 4. The soil layers are horizontal except for the first and second ones which are slightly inclined. The example is analysed using ACO R in dry conditions for two cases: (1) No earthquake loading (case 1), (2) Earthquake loading (pseudo-static horizontal coefficient of 0.1) (case 2). Table 1 summarises the results of both loading cases and compares them with those of optimisation techniques used by other researchers. Figure 4 reports the critical slip surface obtained and the one determined by the commercial software Geoslope [35]. The minimum FS obtained for case 1 i.e without earthquake loading is 1.293 with standard deviation of 0.008 and 1.010 in the presence of earthquake loading, with a standard deviation of 0.025. The proposed algorithm gives lower factors of safety than the GA, However, computational experiments indicate that, for practical purposes, the desired accuracy can be achieved after 200 iterations (i.e., after a number of 10,000 fitness function evaluations), and thus, it is computationally similar to ABC. The fitness function evaluations when compared to discrete ACO (150,000) [17] are significantly reduced, thereby leading to reduction in computation time. Furthermore, the computational time, determined by calculating the average time for 10 computations, was assessed on the basis of an Intel(R) Core(TM) i5-7500 CPU 3.40 GHz processor with 4.00 GB RAM. For one simulation, average computing time is approximately 40 s for the ACO R algorithm.

Engineering application
The proposed methodology is also applied to a case study of an embankment failure, located in eastern India. The National Highway no. 316 is a very important highway in eastern India and a part of it connects the cities of Bhubaneswar and Puri (Fig. 5a) is constructed over soft clay deposits. The existing two lane national highway is being augmented in capacity by increasing the number of lanes. Furthermore, a rail over bridge was proposed in the place of barrier crossing. To take the road to the elevation of the rail over bridge, an approach embankment was proposed with a maximum height of 11 m. The embankment is 27 m wide at the top and 35 m at the bottom with a side slope of 70°. During construction, signs of distress were observed on the left hand side of the embankment, when the embankment height reached its maximum value. After 2 months, the approach embankment failed (Fig. 5c, d) at Samajajpur at about 12.8 km from the city of Puri.

Geotechnical investigations
Five boreholes concentrated primarily on the right hand side were drilled for the initial investigations. After failure, 23 boreholes were drilled both on the left hand and right hand side (Fig. 5a). The construction was stopped as heaving was seen at about 16 m from the embankment (see Fig. 5b-d). The subsoil was throughly characterised for index properties and undrained shear strength of the foundation soil. Figure 6 shows the variation of undrained shear strength as obtained from initial and post failure investigations, as well as their respective average profiles. It can be observed that the left hand side soil profile is weaker than the soil profile in right hand side. Average values of soil strength parameters were used in the analysis for both the cases considered (first and second investigation) obtaining the results reported in Fig. 7.

Analysis with ACO R
In this case, the algorithm was modified as the entry point x 1 is not variable and is taken as the middle point of the road ( x 1 = 0). Indeed, the pictures of the site after the failure reveal that the entry of the failure surface was in the middle of the embankment. The factor of safety was calculated using Slope/w [35], a Limit Equilibrium Method software using Morgenstern Price method and a Finite Element Method using strength reduction module in RS2 [37]. Table 2 compares the results from the analysis, and it can be seen they closely match with each other. The ACO R analysis provides a FS range whose lower value is from 3 to 9% lower than the FS of the other analyses and whose upper value has a difference of 6% at the most from the other analyses.
It can be concluded that FS is more than 1 if we consider the results of the initial investigation (as designers did) with LEM, FEM and ACO whereas FS < 1 if we consider the results of the investigation which was carried out after the failure of the embankment. Figure 7 reports the critical slip surface obtained by the calculation carried out using ACO R , obtaining the safety factors reported in Table 2. The effect of geogrids and pore water pressure was not considered in the analysis for simplification but could be incorporated with minor modification in the Matlab sub-routines.

Summary and conclusions
In this paper, we presented a technique for slope stability assessment using continuous ant colony optimisation algorithm ( ACO R ). The algorithm is tested with literature benchmark example and with data of landslide failure coming from slope failure at a road embankment site. The parameters used in the ACO R algorithm are calibrated to avoid trapping into local minima leading to sub-optimal results. We also presented an  application of our technique applied to one heterogeneous slope having one weak soil layer. The results of the proposed technique are compared with reported results of GA, PSO, HS, ABC and discrete ACO for the chosen literature example. The performance of ACO R was found to be better than that of the other approaches in terms of computing the minimum FS and determining the critical slip surface. Such characteristics make the proposed algorithm suitable for practitioners to identify the critical non-circular slip surface in real slope stability problems.