Application of Harmony Search Algorithm to Slope Stability Analysis

: Slope stability analysis is undoubtedly one of the most complex problems in geotechnical engineering and its study plays a paramount role in mitigating the risk associated with the occurrence of a landslide. This problem is commonly tackled by using limit equilibrium methods or advanced numerical techniques to assess the slope safety factor or, sometimes, even the displacement ﬁeld of the slope. In this study, as an alternative approach, an attempt to assess the stability condition of homogeneous slopes was made using a machine learning (ML) technique. Speciﬁcally, a metaheuristic algorithm (Harmony Search (HS) algorithm) and K-means algorithm were employed to perform a clustering analysis by considering two different classes, depending on whether a slope was unstable or stable. To achieve the purpose of this study, a database made up of 19 case studies with 6 model inputs including unit weight, intercept cohesion, angle of shearing resistance, slope angle, slope height and pore pressure ratio and one output (i.e., the slope safety factor) was established. Referring to this database, 17 out of 19 slopes were categorized correctly. Moreover, the obtained results showed that, referring to the considered database, the intercept cohesion was the most signiﬁcant parameter in deﬁning the class of each slope, whereas the unit weight had the smallest inﬂuence. Finally, the obtained results showed that the Harmony Search algorithm is an efﬁcient approach for training K-means algorithms.

In another study, Michalowski (2002) introduced stability charts for uniform slopes which, according to the kinematic method of limit analysis, lead to a strict lower bound on stability number c/γH or an upper bound on the safety factor. In the proposed charts, pore water pressure and seismic forces were considered for slopes [30]. However, some applications in the field of slope stability have also been recently presented using artificial intelligence approaches [31][32][33][34][35]. Mishra et al. (2019) conducted slope stability assessment based on several recent meta-heuristic techniques. They compared the results with the antlion optimizer and made some recommendations for slope stability analysis [36]. Yuan and Moayedi (2019) investigated slope stability analysis and failure recognition using several artificial intelligence techniques. They found that the combination of multilayer perceptron with biogeography-based optimization can provide a higher performance capacity in the stability assessment compared to other intelligence methods [37]. In another study, Mishra et al. (2020) carried out a slope stability assessment using the multiverse optimization (MVO) algorithm. The obtained results indicated that this algorithm can be applied as a powerful tool for investigating the factor of safety [38].
In the present study, machine learning (ML) was employed as an alternative approach for a preliminary evaluation of the stability condition of slopes. In particular, the harmony search (HS) algorithm was employed. Traditional A.I. techniques, such as linear programming, non-linear programming, and dynamic programming, have been commonly employed for analyzing optimization problems. These techniques can guarantee global optima in simple and ideal models. However, they present some drawbacks: in linear programming, considerable losses occur when a linear ideal model from a non-linear real world problem is developed; in dynamic programming, an increase in the number of variables would exponentially increase the number of evaluations of the recursive functions and tax the core memory; in non-linear programming, if the functions used in computation are not differentiable, the solving algorithm may not find the optimum. Careful attention is also required in the selection of the initial values in order to guarantee convergence to the global optimum and not to the local optima. The HS algorithm is employed in the present study for the purpose of overcoming the above mentioned limitations. It should be mentioned that the utilization of the clustering method for the evaluation of slope stability using the HS algorithm has not been used in previous studies. Additionally, in this study, the HS algorithm was used due to the smaller number of control parameters and also because it is easier to perform than many optimization algorithms. Even though this approach cannot overtake the capabilities of the methods commonly employed in geotechnical engineering (such as limit equilibrium methods or the strength reduction method), it could be taken into account to provide a primary evaluation of the slope safety factor.
The remainder of this paper is organized as follows: Section 2 describes the harmony search (HS) algorithm. Section 3 explains the reasons for selecting the main geotechnical parameters and geometrical characteristics for this analysis and introduces Case study. In Section 4, the necessary analyses are performed to achieve the optimal model and the efficiency of the model is discussed. Section 5 concludes the paper and recommends future studies.

The Harmony Search (HS) Algorithm
In the field of artificial intelligence, there is a wide range of optimization algorithms and techniques used in various fields of industry and academia [39][40][41][42]. In computational science, meta-heuristic algorithms play a key role in solving complex problems [43][44][45][46][47]. The harmony search (HS) algorithm is a meta-heuristic algorithm and was introduced by Geem [48][49][50]. The HS algorithm is an efficient algorithm in unpredicted and uncertain systems that can provide a higher performance capacity for dealing with linear and nonlinear data sets. In addition, the HS algorithm has been successfully applied in the engineering and academic sectors [51][52][53][54][55][56][57]. The working principles of optimization of HS algorithm are similar to the composition of songs by musicians using various notes. In fact, this algorithm is inspired by the process of producing a new song. In this process, the composer produces a new song by combining pieces and music notes, and in the next song, according to the composer's limited memory, the best notes and pieces are considered [58]. The basic HS algorithm comprises four main steps, as follows: Step 1: Initialize the harmony search memory (HM). A number of randomly generated solutions is produced according to the conditions of the problem. If there is an n-dimension Land 2021, 10, 1250 3 of 12 problem based on Equation (1), a HM with a size equal to the HMS can be introduced as follows:  HMS) is considered an answer candidate.
Step 2: Create a new solution from the HM.
Step 3: The new solution from the last step is evaluated and then the HM is updated. If the new solution is better than the worst solution stored in the HM, the two solutions are swapped.
Step 4: The algorithm process and calculations are repeated. In fact, the loop of steps 2 and 3 repeats until the stop conditions are achieved, such as the desired precision level or maximum iteration [58].
The flowchart of the HS algorithm is shown in Figure 1 [49].  HMS) is considered an answer candi Step 2: Create a new solution from the HM.
Step 3: The new solution from the last step is evaluated and then the H If the new solution is better than the worst solution stored in the HM, the are swapped.
Step 4: The algorithm process and calculations are repeated. In fact, the 2 and 3 repeats until the stop conditions are achieved, such as the desired p or maximum iteration [58].
The flowchart of the HS algorithm is shown in Figure 1 [49].

Case Study
In this study, the described HS algorithm was applied in the field of slope stability. The novelty of this research work consisted of considering the nature of complex and uncertain systems and their full compliance with meta-heuristic algorithm concepts. Various parameters have been used by many researchers to analyze the stability of slopes in artificial intelligence methods. Table 1 shows some of the most effective research works in this area. Figure 2 demonstrates the frequency of the parameters affecting the stability of the slope used in these research works. Hence, according to the importance of the parameters used in Figure 1, slope height (H), slope angle (α), cohesion (c'), angle of shearing resistance (ϕ'), unit weight (γ), and pore pressure ratio (r u ) were considered for modelling in this study. Table 1. Literature review of machine learning techniques applied to slope stability problems. Referring to the slope whose scheme is reported in Figure 3, the following input data were considered in the present study: slope height (H), slope angle (α), intercept cohesion (c'), angle of shearing resistance (ϕ'), soil unit weight (γ) and pore pressure ratio (r u ). A dataset consisting of 19 case studies was considered randomly, referring to the circular slip surface. For the considered slopes, the above-mentioned parameters were available, as well as the factor of safety (Table 2) [80].

References Geotechnical Properties and Geometry Slip Surface Number of Cases
In order to classify slope stability cases, the k-means algorithm (Lloyd's algorithm) was used as it is one of the most common and popular clustering methods optimized by the HS algorithm. Equation (2) shows the k-means algorithm [81]: where x i represents a member of a cluster with i = 1, 2, 3, . . . , n. In addition, m j and j are the center of each cluster and the number of clusters, respectively, with j = 1, 2, 3, . . . , k.
Moreover, d expresses the Euclidean distance between each member of a cluster (x i ) and the center of the same cluster (m j ). In fact, the fitness function can have two main purposes: maximizing the distance between the centers of the clusters and minimizing the distance between the members of a cluster. Usually, when the goal of optimization is minimization, the fitness function is also called cost function. Therefore, in this study, the purpose of optimization was to minimize the distances between each member of each cluster, and the value of the cost function was actually the minimum distance between each member and the center of each cluster. Generally, when the values of the cost function are close together across two consecutive iterations, it can indicate convergence of the algorithm, and when these values are equal, it means that the error rate is zero and the algorithm has reached convergence.
Land 2021, 10, x FOR PEER REVIEW 5 of 12 Figure 2. Frequency of geotechnical and geometrical parameters studied in the slope stability analysis research works documented in Table 1.   Figure 2. Frequency of geotechnical and geometrical parameters studied in the slope stability analysis research works documented in Table 1. nd 2021, 10, x FOR PEER REVIEW 5 Figure 2. Frequency of geotechnical and geometrical parameters studied in the slope stability a ysis research works documented in Table 1.

Modelling by the HS Algorithm
The current study focused on the applicability of a combination the HS algorithm with the K-means algorithm as an unsupervised clustering technique in order to cluster and evaluate slope stability according to compiled datasets from case studies of the circular critical failure mechanism.
In the first step of clustering, all data were normalized. For this purpose, all data belonging to each category were divided by the largest-value data of the same category [82]. Then, determining of control parameters of the algorithm was necessary. The control parameters of algorithms contributed greatly to finding the optimum results, and increased the convergent speed of the algorithm. There are often no special equations and rules, so, based on a trial-and-error procedure, expert opinions and data sets are considered [83]. After the initial analysis, the control parameters were determined, including maximum iterations equal to 300, and a pitch adjustment rate (PAR) equal to 0.1, i.e., the neighboring values were selected with a 10% probability. In addition, based on previous studies, the value of the harmony memory size (HMS) is normally considered to be within the range of 50 to 100. Specifically, a value of 60 was used in this application following a trial-and-error procedure [58,[83][84][85]. Two classes were considered: the class with the label "1" included unstable slopes (FS ≤ 1), whereas stable slopes (FS > 1) belonged to the class with the label "2". The process of optimization based on iterations, and the results of clustering for the 19 cases with two classes, are shown in Figure 4 and Table 3, respectively.  According to Figure 4, the best cost for the clustering process was reached after 55 iterations and equaled 8.209. Since the best cost represents the performance of the algorithm, Figure 4 indicates an appropriate convergence speed of the HS algorithm. It is worth mentioning that the best value is a dimensionless number, and that if the difference between the results of two consecutive iterations is less than the minimum acceptance precision, the process of optimization will be fixed. In the clustering analysis, whose results are summarized in Table 3, all cases were classified into two classes based on the results of the optimum partition. The optimum partition, determined the Euclidean distance of each case study from the center of each class and then the class of each case study, is indicated in the fourth column (Defined class). If the class of sample conformed to the initial assumption, the clustering was verified as "Satisfied", otherwise it was not verified, and labelled as "Not Satisfied". For instance, Case 3 had Euclidean distances equal to 0.474 from Class 1 and 0.352 from Class 2; hence, this case belonged to Class 2. Besides this, Case 15 had Euclidean distances equal to 0.392 from Class 1 and 0.751 from Class 2; hence, this case belonged to Class 1. In the fifth column, the values of the factor of safety for the considered slopes are indicated, and  According to Figure 4, the best cost for the clustering process was reached after 55 iterations and equaled 8.209. Since the best cost represents the performance of the algorithm, Figure 4 indicates an appropriate convergence speed of the HS algorithm. It is worth mentioning that the best value is a dimensionless number, and that if the difference between the results of two consecutive iterations is less than the minimum acceptance precision, the process of optimization will be fixed.
In the clustering analysis, whose results are summarized in Table 3, all cases were classified into two classes based on the results of the optimum partition. The optimum partition, determined the Euclidean distance of each case study from the center of each class and then the class of each case study, is indicated in the fourth column (Defined class). If the class of sample conformed to the initial assumption, the clustering was verified as "Satisfied", otherwise it was not verified, and labelled as "Not Satisfied". For instance, Case 3 had Euclidean distances equal to 0.474 from Class 1 and 0.352 from Class 2; hence, this case belonged to Class 2. Besides this, Case 15 had Euclidean distances equal to 0.392 from Class 1 and 0.751 from Class 2; hence, this case belonged to Class 1. In the fifth column, the values of the factor of safety for the considered slopes are indicated, and a comparison was made between the fourth and fifth columns. For instance, in Case 10, the Euclidean Land 2021, 10, 1250 8 of 12 distances were calculated as 1.053 and 1.158 from Class 1 and Class 2, respectively; hence, Case 10 should belong to Class 1, while it belonged to the Class 2 according to the value of the FS-therefore, it could be concluded that this case was wrongly classified. Afterwards, the accuracy of the results provided by the HS algorithm were verified, as shown in the last column. As can be seen, it can be concluded that the model was able to correctly classify 17 of 19 case studies, with a model accuracy equal to 89%.

Results and Discussion
In this study, the 19 case studies were classified into two classes by a combination of the HS algorithm and the K-means algorithm. As mentioned earlier, the accuracy of this modelling was 89%, with only two cases that were classified incorrectly. Consequently, it was found that the combination of HS algorithm and K-means can provide high performance capacity in clustering and evaluating slope stability. Figure 5 shows the result of clustering.
Land 2021, 10, x FOR PEER REVIEW 8 of 12 a comparison was made between the fourth and fifth columns. For instance, in Case 10, the Euclidean distances were calculated as 1.053 and 1.158 from Class 1 and Class 2, respectively; hence, Case 10 should belong to Class 1, while it belonged to the Class 2 according to the value of the FS-therefore, it could be concluded that this case was wrongly classified. Afterwards, the accuracy of the results provided by the HS algorithm were verified, as shown in the last column. As can be seen, it can be concluded that the model was able to correctly classify 17 of 19 case studies, with a model accuracy equal to 89%.

Results and Discussion
In this study, the 19 case studies were classified into two classes by a combination of the HS algorithm and the K-means algorithm. As mentioned earlier, the accuracy of this modelling was 89%, with only two cases that were classified incorrectly. Consequently, it was found that the combination of HS algorithm and K-means can provide high performance capacity in clustering and evaluating slope stability. Figure 5 shows the result of clustering. In order to evaluate the effect of each parameter on slope stability, referring to the dataset considered in this study, the Euclidean distance of each effective parameter from the centre of each class was determined and measured. The results are shown in Table 4. According to Table 4, for all cases in Class 1 (failed class), the cohesion (c') was the most significant parameter, with an Euclidean distance equal to 0.063, followed by the pore pressure ratio (ru), slope height (H), unit weight (γ) and slope angle (α), with Euclidean distances equal to, respectively, 0.714, 0.736, 0.753 and 0.812-contributing greatly to the existence of instability in the slopes. Finally, the angle of shearing resistance (φ') had the lowest impact on the stability of the analyzed cases. It is worth pointing out that this latter statement holds only for the considered dataset and does not have general validity. In Class 2 (the stable class), the cohesion (c') had the most significant effect on all the cases. It is worth mentioning that the cohesion (c') played a key role in the stability of slopes in both classes, including the failed and stable classes. In addition, the unit weight In order to evaluate the effect of each parameter on slope stability, referring to the dataset considered in this study, the Euclidean distance of each effective parameter from the centre of each class was determined and measured. The results are shown in Table 4. According to Table 4, for all cases in Class 1 (failed class), the cohesion (c') was the most significant parameter, with an Euclidean distance equal to 0.063, followed by the pore pressure ratio (r u ), slope height (H), unit weight (γ) and slope angle (α), with Euclidean distances equal to, respectively, 0.714, 0.736, 0.753 and 0.812-contributing greatly to the existence of instability in the slopes. Finally, the angle of shearing resistance (ϕ') had the lowest impact on the stability of the analyzed cases. It is worth pointing out that this latter statement holds only for the considered dataset and does not have general validity. In Class 2 (the stable class), the cohesion (c') had the most significant effect on all the cases. It is worth mentioning that the cohesion (c') played a key role in the stability of slopes in both classes, including the failed and stable classes. In addition, the unit weight (γ) had the lowest impact on all cases in the second class, with a Euclidean distance equal to 0.844. On the one hand, one of the most important strengths of clustering techniques as unsupervised methods is that it is possible to analyze problems using a limited amount of data. On other hand, in some geotechnical problems, it is difficult to prepare laboratory and field data; therefore, based on the results of this study, it can be concluded that a combination of the HS and k-means algorithms as a clustering approach can be used as a reliable system for modelling and evaluating some problems involved in geotechnical engineering. In addition, it should be noted that this model is unique, and it cannot be used directly in other case studies.

Conclusions
This paper presents an application of the machine learning technique to slope stability analysis by using a combination of the HS and K-means algorithms. For this propose, 19 case studies documented in the literature were considered, for which the most relevant geotechnical parameters and geometrical characteristics were collected. Specifically, these were the slope height (H), slope angle (α), intercept cohesion (c'), angle of shearing resistance (ϕ'), soil unit weight (γ) and pore pressure ratio (r u ). Besides this, all the considered slopes were characterized by a circular slip surface. A clustering analysis was conducted by defining two classes: the first class included unstable slopes (FS ≤ 1), whereas the second class concerned stable slopes (FS > 1). Then, a comparison was made between the results of the clustering analysis and the actual FS. The figures show that 17 of the 19 case studies were classified correctly into the two classes, with a model accuracy equal to 89%. In addition, for case studies considered in this work, the intercept cohesion (c') was the parameter that exerted the greatest influence on slope stability for both classes (failed and stable). The purpose of this manuscript was to demonstrate the capability of the HS algorithm of dealing with problems of slope stability. Even though this approach cannot outperform the methods currently commonly employed in geotechnical engineering (such as limit equilibrium methods or the strength reduction method), it could be taken into account for providing a primary evaluation of slope stability conditions. From a general point of view, the HS algorithm should be employed by using a much larger dataset to 'learn' how the factor of safety is affected by input parameters (geotechnical and geometrical). Afterwards, the algorithm can be used to perform a preliminary evaluation of the stability conditions of a slope that is not already included in the database. In this context, the dataset could be built referring to a set of slopes included in a particular geographical region.