Search of Non-circular Slip Surface Based on Improved FOA

: The determination of the most critical non-circular slip surface can be attributed to the optimization of complex non-linear multi-peak function due to its numerous control variables and large amount of calculation. It is a trend in recent years to apply intelligent optimization algorithm into slope stability analysis. Considering that the standard Fruit fly Optimization Algorithm (FOA) is prone to fall into local extremum, the improved Fruit fly Optimization Algorithm is obtained by incorporating the standard FOA with the simulated annealing idea. In order to improve the search efficiency, a fixed step size is adjusted to an adaptive step size, and a double tier search strategy is proposed to be applied: the potential non-circular slip surface is obtained from the outer layer, and the factor of safety along the potential slip surface is calculated step by step from the inner layer. The improved FOA is applied to a slope with weak interlayer. The feasibility, superiority and efficiency of the improved algorithm are proved by comparing its answers to the judges'. Different inter-slice force functions, various initial values of Fs and λ are assumed, and the results show that these parameters could hardly affect final solution for safety factor.


Introduction
Slope stability analysis has always been a research subject of great value in geotechnical engineering. Generally, a series of sliding surfaces are assumed, and the safety factor of a slope is calculated based on limit equilibrium method or strength reduction method, then the minimum safety factor is found out to judge the stability of the slope, and the corresponding sliding surface is named as the critical slip surface. The shape of the most dangerous slip surface could be summarized as circular or non-circular. The research on circular slip surface search has become mature [1], while the search of non-circular slip surface remains a problem, which could be attributed to an optimization problem of complex nonlinear multi-peak function due to its numerous control variables and large amount of calculation. Scholars have also put forward various methods to solve this problem. With the development of computer technology and the penetration of artificial intelligence, it has become a research trend in recent years to apply intelligent optimization algorithm into search of the critical slip surface of a slope.
Cao Ping [2] used an improved genetic optimization algorithm to search for any sliding surface in a slope, and got a result close to the actual situation. Gao Wei [3] developed the encounter ant colony algorithm based on ants' forward and backward foraging meeting which formed a complete path, and applied it to search for a non-circular sliding surface. Wu Wenyuan [4] put forward a joint search method for a circular and broken line sliding surface based on alternate search method and quadratic interpolation theory. Wei Hainan [5] combined genetic algorithm with MP method to establish a global optimization method for any sliding surface search. Based on critical slip field theory, Jiang Zefeng [6] established a critical slip field analysis and calculation system considering different hydraulic conditions and reinforcement conditions. Wei Changyi [7] reduced multiple indexes of data space to a few principal components through linear transformation, and established a slope stability analysis model based on genetic neural network. Shahrokhabadi [8] combined genetic algorithm with natural element method to determine a most dangerous non-circular sliding surface, and the algorithm was proved to have good robustness through an example solution.
Khajehzadeh [9] applied adaptive maximum speed constraint to universal gravitation algorithm, and the new improved algorithm exhibits good convergence in search of a most dangerous non-circular slip surface and calculation of the slope reliability index. Intelligent optimization algorithm has achieved great advantages in slope stability analysis, but there are still some shortcomings, such as slow convergence rate of the simulated annealing algorithm [10], low search efficiency of the ant algorithm [3], poor parameter universality of the genetic algorithm [10]. In this paper, Fruit fly Optimization Algorithm and Simulated Annealing algorithm were combined to get an improved optimization algorithm, and Improved MP (Morgenstern and Price) method was applied to calculate the safety factor of a slope. Meanwhile a double-tier search strategy for non-circular sliding surface was proposed and verified by an example slope with a soft interlayer.

Improved MP Method
Morgenstern and Price method (MP method) is a strict slice method based on limit equilibrium theory, which satisfies the requirement of force balance in both horizontal and vertical directions and overall moment balance, and it is suitable for safety factor calculation on a slip surface of any shape. The force analysis of one soil slice is shown in Figure 1, where W is gravity, natural unit weight and saturated unit weight are taken respectively while above or below the infiltration line;Q is load on slope surface, ηW is horizontal seismic force; E is horizontal thrust acted on the side of the soil slice, X is shear force acted on the side of the soil slice; N is normal force on the bottom of the soil slice; U is pore water pressure; S , , , / is anti-sliding force on the bottom of the soil slice, , and , are effective internal friction angle and effective cohesion, respectively, and is the safety factor. , and , are shear strength parameters derived from strength reduction method. In the original MP method reasonable initial values for λ and Fs were firstly assumed, and then integral solutions for 、 " were obtained one by one, until the results of # 、" # are zero, which is not conducive to programming. According to the improved MP method proposed by Zhu Dayong [11], the external forces on the soil slice are decomposed along the direction perpendicular and parallel to the sliding surface, and the following formula can be obtained: , . % /. ' 0 % Where,is the sum of shear resistance on the soil slice except for the inter-slice force , and , is the sum of sliding force. Assume, Considering that the inter-slice force on the two ends of the slope is zero, e.g. E # 0, safety factor can be expressed as: Considering moment balance of a single soil slice and no bending moment at two ends of the slope, e.g. M 0 =M n = 0, the expression of the scale factor is obtained as follows: In the improved MP method, the expressions of 、λ are implicit equations, which are conducive to iterative solution with MATLAB.

Simulation of Non-circular Slip Surface
The shape of non-circular sliding surface includes straight line, logarithmic spiral, quadratic parabola, polyline, etc.. In practical engineering, the shape of sliding surface is often irregular due to the complexity of engineering geological conditions, the difference of slope geometry and material parameters, thus it will be more extensive and objective to adopt arbitrary shape for a sliding surface. As the broken lines can approach any form of sliding surface infinitely, they'll be taken during non-circular sliding surface simulation process.
As shown in Figure 2, the broken line slip surface is composed of a series of discrete points connected by a straight line. If there are n vertices to be connected, and the coordinates of each control point are taken as variables, then there would be 2n control variables in a two-dimensional plane problem. Even if the longitudinal coordinates of the endpoints of a slip surface are given, there are still 2n-2 control variables. Obviously, the number of control variables increases with the number of discrete points, and the search workload increases exponentially with the number of control variables. Therefore, the first step is to determine a reasonable number of variables for slip surface simulation.

Figure 2. Sliding surface formed by broken lines.
Take the control variables of the broken line sliding surface as the coordinates (x, y) of each point, as shown in Figure 2, determine range of W X 、W Y , which respectively stands for horizontal coordinates of the exit and entrance point of the broken line sliding surface, and corresponding vertical coordinates Z X 、Z Y can be determined solely by the slope equation, W U 、 W [ … W #*+ are obtained by dividing the horizontal distance between W X and W Y evenly, thus any broken line sliding surface can be fitted by n variables in vector S. The generation of broken line sliding surface is random, which may produce some unreasonable sliding surfaces, such as "zigzag" and "non-concave" sliding surfaces shown in Figure 3. These unreasonable slip surfaces will not only waste resources of program and reduce search efficiency, but also cause the non-convergence problem during safety factor calculation. In order to avoid sliding surfaces of "nonconcave" or "zigzag", the control points for broken line sliding surface have to meet following requirements: 1) for abscissa of adjacent control points: W B+^W .
2) angle between connecting line of adjacent control points and horizontal direction α ∈ $&45°, 80°), and α B+^α . The steps to fit a polyline sliding surface are as follows,  Figure 4, DF is an extension of prior sliding surface, and vertical coordinate of point F can be used as the lower limit Z e # for vertical coordinates of node i, meanwhile the lower limit should ensure that broken line sliding surface does not go through the stiff soil layer which lies below. The intersection of DB and vertical line through node i is G, which can be used as the upper limit Z e=f for node i. 3) To calculate the angle α between the line connected by two adjacent control points and the horizontal plane, and "zigzag" and "non concave" cases should be avoided.

Standard Fruit Fly Optimization Algorithm (FOA)
Standard fruit fly optimization algorithm (FOA), developed by Pan Wenchao [12], based on the foraging behavior of fruit flies, is a new global optimization algorithm of swarm intelligence, which has been widely used in many fields to solve optimization problems in recent years. Zhang Shuiping summarized the design ideas of FOA, focusing on the existing improvement methods and related applications [13]. The mechanism of fruit fly optimization algorithm is based on the simulation of the foraging process of fruit fly population. The good smell and vision of a fruit fly is utilized, the corresponding sensory search link is designed, and the optimal solution is achieved through continuous iteration. The feeding iteration path of fruit fly population is shown in Figure 5.
FOA is exactly to internalize the foraging behavior of fruit fly into continuous iterations approach to the target value, thus to achieve the optimization of solution space. Therefore, the algorithm has advantages such as simple structure, easily understood and strong operability.

Improved FOA
The essence of standard fruit fly optimization algorithm is to copy all individuals in the population into the one with the current best adaptability, e.g. all individuals fly to the best location, and then take it as the starting point for iterative updating of the next generation. There is a natural shortage in the optimization mechanism, if the previous generation of fruit fly falls into local extreme value, the next generation which learns from the prior one will also fall into local optimum, and so on, eventually the fruit fly population falls into local extremum, that is, the optimization failed [14]. In addition, a fixed step size is generally used in the standard fruit fly optimization algorithm, which ignores the difference of the current position of individual fruit fly, and the overall search ability of fruit fly population has been limited. In order to overcome these problems, simulated annealing (SA) is integrated into the fruit fly optimization algorithm to increase the probability of jumping out of local extremum, and the fixed step size in the original algorithm is adjusted to adaptive step size according to the position difference among individuals. Thus an improved fruit fly optimization algorithm SA-FOA is proposed.

Step Adjustment
In standard fruit fly optimization algorithm, the search step size for an individual is: Where Value is an initial step and is a fixed one, Random is a random number from interval [0,1], so range of search step size is [0, Value]. The problem here is that for individuals with smaller search step, convergence is slow and prone to fall into local extremes, while for those with larger step size, there will be a fluctuation effect in later stage of search. Generally speaking, the larger step in early stage of search can speed up convergence process, and the smaller one can improve convergence accuracy in later stage of search, which is beneficial to obtain an optimal solution. Based on this, the idea of inhomogeneous variation is applied into step setting, thus self-adaptive search step is determined according to current position difference of individuals, g ( * i1 & $ k & 1 /l Wk ) 6 m * %k 0. 5 & o (10) where N is the initial step, set according to maximum iteration number maxgen, g is current iteration number, is inhomogeneous variation factor, and sign is a symbolic function which returns positive or negative symbols to parameters. Self-adaptive step is adjusted from two aspects on the basis of the initial step, including inhomogeneous variation step and plus-n-minus. Inhomogeneous variation can ensure that step decreases gradually with increasing of iteration number, which meets the requirements for proper search step in different search stage in FOA, and plus-n-minus can guarantee FOA to jump out of local extremum and obtain global optimal solution.

FOA Based on Simulated Annealing (SA-FOA)
The idea of simulated annealing (SA) was first proposed by Metropolis in 1953 [15], and applied to combinational optimization by Kirkpatrick in 1983 [16]. SA algorithm is a global optimal method developed according to principle of solid annealing in physics. When solid is heated up, internal particles are in disorder, and internal energy gradually increases, while cooling down, particles arrangement tends to be stable, and internal energy decreases. Assume at state A, internal energy of the system is E(A), after external disturbance, at state C, internal energy becomes E(C). Based on Metropolis criterion, probability that particles tend to balance at temperature T is: p q 1, r s t ; q Wv w& P x *P X N y , r z t It can be seen from formula (11) that while solid annealing, internal energy decreases and system certainly accept the state, but if new state increases internal energy, then system accepts the state with a certain probability, that is, in the optimization process, optimal solution and poor solution of certain probability could be obtained. The internal energy function E(X) is evolved into target function, the annealing temperature T evolves into control parameter, and with the variation of T, the algorithm executes evolution of "new solution-judgment-abandonment or acceptance" until optimal solution is obtained.
The search process of SA algorithm is a heuristic random search based on Monte Carlo iterative solution, which is an extension of local search algorithm, and Metropolis criterion is key to its convergence to global optimal solution, so compared with FOA, SA algorithm hardly falls into local optimum. Theoretically, as long as there is enough time, it is sure for converging to global optimal solution with probability of 1. When simulated annealing idea is incorporated into standard fruit fly optimization algorithm, SA-FOA is obtained with following steps: 1) Initial setting of parameters related to fruit fly population: location X_axis, Y_axis, scale Sizepop, maximum iteration time Maxgen. 2) Update random distance (Random Value) and orientation during individual foraging.
3) Calculate distance between the individual and the origin, and determine taste concentration. accept the disturbance, and take the location (Xˊ, Yˊ) and taste concentration Smellˊ after disturbance as optimal value for the iteration, else reject the result produced by the disturbance, and take the former location and taste concentration as optimal value for the iteration, or wait until iteration number reaches set Maxgen to abort the program.

Double Tier Search Strategy
In non-circular slip surface search, there are still n control variables to be determined after simplification. If there are less soil slices, search workload can be reduced, but optimal solution may be missed, if there are more soil slices, search workload will be increased, also more freedom will be included, the objective function (safety factor) will have more local extremum in independent variable space. It is more difficult to obtain global extremum across many extremum. In order to solve the contradiction between quantity of soil slices and efficiency and accuracy of slip surface search, a double tier search strategy is put forward in this paper as shown in Figure 6. Since there are many control variables of non-circular slip surface, which includes W X , Z U , Z [ , Z … W Y , safety factor Fs and proportional coefficient λ , and the relationship between control points coordinates and Fs and λ belong to progressive level. It is inappropriate to optimize control points coordinates, Fs and λ at the same time. Thus in this paper, potential slip surface and its safety factor were proposed to be sought and determined through internal and external circulation. At outer layer through constraints and penalty function mechanism, SA-FOA is used to search different slip-in/out points, and intermediate points for obtaining multiple sets of potential slip surfaces, and outer circulation is aborted by maximum loop count. At inner layer Newton-Raphson numerical iteration is adopted to calculate Fs and λ with improved MP method for those determined slip surfaces. Both the differences between two Fs, λ from two adjacent iterations should be within allowable error ε, which determines whether the program should be ended or not.
With double-tier search strategy, the process looking for non-circular slip surface is as following, 1) Set geometric parameters such as slope height, slope angle, and mechanical parameters such as cohesion and inner friction angle. 2) According to experience set calculation accuracy, select control points coordinates to divide soil slices. 3) Calculate self-weight W , anti-sliding force R and slide force T for each soil slice. For non-homogeneous slope, the bottom of soil slice may go through different soil layers, when average value of parameter is calculated according to proportion of the slice bottom length located in different soil layers. The average parameters " , " can be defined as:

5) Update
( , λ λ ( , calculate 2 ,and ( , then calculate E , λ ( . 6) Observe whether the differences between new ( , λ ( and former corresponding , λ are within allowable error -+ and -U respectively. If they are, then exit from program, and return renewed ( and λ ( as the result from inner circulation. If they are not, then check if the maximum iteration number is reached in order to avoid "endless loop" of the program. 7) If the maximum iteration number is reached, then current slip surface is unreasonable, and its result of safety factor should not be included in the optimization scope. In this case, the maximum value set for Fs, such as 100, will be taken as return result from inner circulation based on the idea of "penalty function", then exit this optimization.

Example Solution
The example is selected from ACADS assessment question EX3, a slope containing a weak interlayer, whose profile features and material parameters are as shown in Figure 7, and the weak interlayer is above groundwater table, and pore water pressure is zero. The assessment question is divided into two parts, (a) to determine the most dangerous broken-line slip surface and its corresponding safety factor, (b) to determine safety factor of a specified slip surface, coordinates for the specified slip surface is shown in Table 1. The slope contains a weak interlayer of 0.5m thickness, while the most dangerous slip surface may not have to pass through the weak interlayer, when slip surface intersects the weak interlayer, the weak interlayer between entrance and exit of intersection is considered to be a section of the broken-line slip surface. and generated broken-line slip surface has to be above lower boundary of the weak interlayer.  For problem (a), the minimum safety factor is 1.177 by the proposed algorithm, the slip surface is shown in Figure 8, and the iterative process is shown in Figure 9. It can be seen that the improved fruit fly optimization algorithm gradually approaches the minimum value of safety factor by searching various slip surfaces. Due to implementation of double tier search strategy, with improved MP method to determine safety factor, it only takes 13 iterations to be close to optimal solution, and 30 iterations to converge to the optimal solution, which shows the application of double tier search strategy has greatly improved search efficiency. Compared to the recommended referee answer which is 1.26, the error from this paper is 6.6%. For problem (b), the improved algorithm is applied to calculate safety factor on the specified slip surface, and the result is 1.47. the iteration process is shown in Figure 10. Since control points coordinates of the slip surface are specified, only the inner layer search needs to be executed, i.e. safety factor is calculated using Newton-Raphson numerical iteration. Compared to the recommended referee's answer 1.34, its error is 9.7%.

Discuss on Effect of Initial Value
With different initial value of Fs, λ and the inter-slice force function f(x), safety factors of given slip surfaces are calculated, results are shown in Figure 11, Figure 12 and Figure 13. It only took about 15 iterations to converge to the same value for safety factor in each parameter assignment.   It can be seen that different values for initial Fs, λ and f(x) hardly affect the final convergence value for safety factor.

Conclusions
1) Broken-line is used to simulate non-circular slip surface of any shape, which is more conformable to reality, and control point coordinates are taken as independent variables. Combined with the improved MP method, the safety factor along some slip surface could be determined. 2) In order to overcome the shortcoming of standard FOA, the idea of simulation annealing is integrated with FOA, fixed step size is replaced by self-adaptive step, and double tier search strategy is proposed in the improved FOA: the potential non-circular slip surface is obtained from the outer layer, and the factor of safety along the potential slip surface is calculated step by step from the inner layer, which is proved to have high feasibility and efficiency in non-circular critical slip surface search. The improved FOA is also applied to a slope with weak interlayer, and the safety factors solution are very close to the referee's answers, which shows applicability and feasibility of the improved algorithm in non-circular slip surface search. 3) Refer to the improved algorithm, parameter assignment's effect was also been discussed. Different inter-slice force functions, various initial values of Fs (0.8, 1.0, 1.2) and λ (0.2, 0.4, 0.6, 0.8) are assumed, and the results show that these parameters could hardly affect final optimal solution for safety factor, which indicate the improved FOA has a good application prospect in slope stability analysis.