Machine Learning-Aided Monte Carlo Simulation and Subset Simulation

Deterministic Analysis

Over time, several techniques, including Fellenius (Ordinary), Bishop’s, Janbu’s, Spencer, and Morgenstern–Price, have been devised for the stability problem of the slope. Each of these methods are totally based on the different statics of equilibrium (i.e., force equilibrium [FE] and moment equilibrium [ME]) and their assumptions. However, the OMS is completely based on the statics of ME and is an exception as it completely ignores the interslice forces (i.e., interslice normal and interslice shear). In this investigation, OMS has been adopted for deterministic analysis of slope.

In the OMS, the weight of each slice is transformed into forces perpendicular and parallel to the slice base.

The available shear strength for a slope is determined by the base normal force, which is perpendicular to the slice base, whereas the driving force is calculated using the weight component parallel to the slice base. Figure 1 illustrates the dividing the sliding mass vertically into many slices, and also the free body diagram of a particular slice is shown. The FS is determined by dividing the total shear strength available along the slip surface by the sum of the mobilized shear driving forces. The FS is calculated using OMS considering pore water pressure ratio, as illustrated in Equation 1.

FS = ∑ [ c Δ l + N tan ∅ + r u N ] ∑ W sin α ( 1 + k h k v ) = ∑ S resistance ∑ S mobilized

(1)

where c is cohesion, Δ l is slice base length, N stands for base normal (i.e., W cos α ) , r u stands for pore water pressure ratio, ∅ represent angle of internal friction, W represents slice weight, k h is horizontal seismic coefficient, k v is the vertical seismic value, and α is the slice base inclination. Conventionally, it is assumed that when FS is less than 1 slope is regarded as a fail.

OPEN IN VIEWER

Figure 1.

Illustrative typical slope geometry.

Deterministic studies assume that soil attributes have fixed and known values, whereas probabilistic analyses take into account the uncertainty associated with soil properties within a slope. To account for this uncertainty, MCS and SS methods have been utilized in the current study for PA. The following subsection provides an illustration of the use of these methods in PA.

Random Field Modeling in Inherent Variability of Soil

Soil properties can be viewed as stochastic variables used to gauge their variability or uncertainty. Despite their apparent uniformity, natural soil formations often reveal regional disparities attributable to the complex geological processes and historical loading events they have experienced ( 58 ). The random field, showing statistical homogeneity, has been discretized using different methods developed over time. In this work, the variability of soil within a statistically homogeneous layer is depicted by a one-dimensional lognormal field named F(X). This choice is made because of the unambiguous positivity of the soil parameter ( 59 ). The adoption of statistically homogeneous random fields is driven by the necessity to create a practical framework for numerical analysis, given the challenge of accurately representing the true physical homogeneity within natural soil layers. Through adjusting the random field parameters, such as correlation length and variance, in line with empirical soil data, we are able to incorporate the observed physical heterogeneity into our model. This ensures that the model accurately mirrors the real-world variability of soil characteristics. The process of calibrating these parameters, which is confirmed through comparison with independent datasets, significantly improves the model’s ability to predict within the parameters of the assumed framework. In this context, x is a lognormal random variable that exhibits mean value μ x and standard deviation (SD) σ x in a vertical direction which can be simulated using the expression as shown in Equation 2 ( 9 , 38 ).

F ( X ) = exp ( μ lnX s + σ lnX L χ )

(2)

where μ l nX = ln μ x − 1 2 σ l nx 2 and σ X = ln ( 1 + COV X 2 )

are the mean and SD of ln(X), respectively; CO V X is the coefficient of variation representing the ratio of σ x to the μ x value; s  = a vector with unity components; χ is an n-dimensional standard Gaussian vector; L is a lower triangular matrix obtained through Cholesky decomposition, ensuring ρ = [ ρ ij ] = L L T . Here, ρ is the Pearson correlation matrix at different depths, calculated as ρ ij = e ( − 2 | z i − z j | λ ) , where z i and z j denotes the depth at two distinct locations and λ signifies the correlation length. Moreover, the PDF of a lognormal distribution is established using Equation 3.

f ( x ) = 1 x σ ln X 2 π exp [ − 1 2 ( ln x − μ ln X σ ln X ) 2 ]

(3)

Probabilistic Analysis

Several factors can introduce inaccuracy into a slope stability analysis, such as the random spatial distribution of soil properties, the difficulty of accurately determining stratigraphy below the surface, and mistakes in the computer models ( 60 ). Accounting for the inherent uncertainties associated with soil factors is a process that can be systematically addressed using probability theory and statistics. Within this study, the stochastic analysis of slope stability for road or rail embankments is performed, utilizing both MCS and SS techniques. These analyses are aided by the utilization of the “UPSS 3.0 Add-ins” tool embedded within the MS-Excel platform.

Monte Carlo Simulation

The MCS technique was initially developed in the 1940s by John von Neumann and Stanislaw Ulam ( 61 ), and is a numerical procedure for iteratively computing a mathematical or empirical operator with random or uncertain variables and predetermined probability distributions ( 62 ). A slope fails when it slides along any slip surface (i.e., a critical slip surface). Determining the overall failure probability (also known as the system p f ) along each potential slip surface is a challenging mathematical task ( 63 ). MCS is a robust technique to determine the desired level of p f on a slope in an efficient way compared with deterministic analysis. To ensure the desired level of accuracy in the estimation of p f , the minimum number of samples required in direct MCS should be at least ten times higher than the inverse of the desired accuracy level (i.e., 10/ p f ) ( 64 – 66 ). For instance, to attain an accuracy level of 0.001 using MCS, a minimum of 10,000 random samples should be generated. Let s 1 , s 2 , s 3 … s n denote the N potential slip surface that are considered in the LEM analysis of slope stability. The slope failure occurs when any component (i.e., s 1 , s 2 , s 3 … s n ) fails (i.e., FS < 1); then p f is frequently calculated by adding the total number of failed samples to the total number of random samples generated, as demonstrated in Equation 4. Subsequently, β for the p f can be found using Equation 5.

p f = N h N

(4)

β = φ − 1 ( 1 − p f )

(5)

where N h is the total number of cases where failure has occurred, N represents the total number of generated random samples and φ stands for the standard normal cumulative distribution function.

Subset Simulation

The SS is a more sophisticated variant of the MCS that utilizes conditional probability and the Markov chain Monte Carlo simulation (MCMCS) technique to accurately estimate small tail probabilities ( 67 ). It is a technique for quickly and effectively producing samples that match specified levels of failure probabilities in a progressive manner. To achieve the desired failure domain, this approach employs uniquely crafted Markov chains that generate conditional failure samples at intermediate events. This technique expresses the low level of probability for a rare failure event F into sequence of intermediate failure events ( F 1 , F 2 , F 3 , … , F m ) with relatively higher conditional probability. To exemplify, let us examine the slope stability problem, where p f denotes the probability that the performance function T (i.e., the probability of FS less than fs or P [FS < fs]) falls below a specified threshold value t, or p f = P ( T < t ) denoted as t = t m < t m − 1 < , … , < t 2 < t 1 . These threshold values represent critical levels of the performance function T, where exceeding them indicates a degree of failure in slope stability. For instance, t m could represent the most critical threshold, while t 1 could denote a less severe failure condition. In this research, the target FS was considered to be less than 1.1 and threshold value greater than 0.909 (i.e., t = 1/FS), which was considered a failure sample. Therefore, if all samples have a threshold greater than 0.909, it is considered as a fail.

The SS methodology involves m + 1 steps, where the initial step employs a direct MCS to produce “N” unconditional samples. Subsequently, the remaining steps, consisting of m steps of MCMCS to generate conditional samples, can be calculated as N ( 1 − p o ) , where N signifies the number of random sample produced in each step and p o is the conditional probability, which has been reported as 0.1 based on literature review ( 9 , 67 , 68 ). The sum of conditional samples produced by SS in “m” number of simulation levels is { m × N − ( m − 1 ) p 0 N } . Notably, the produced samples fall into m + 1 mutually exclusive and collectively exhaustive subgroups such as { Ψ k , k = 0 , 1 , 2 , … , m } , determined by the intermediate threshold value. More specifically, the subgroup { Ψ o , Ψ 1 , Ψ 2 , … , Ψ m − 1 } includes N ( 1 − p o ) samples out of N samples. However, the last subgroup Ψ m consist of all “N” samples. The p f is then calculated as shown in Equation 6.

p f = ∑ k = 0 m P ( F Ψ k ) P ( Ψ k )

(6)

where Ψ o = { T ≥ t 1 } , k = 0; Ψ k = { t k + 1 ≤ T < t k } , k = 1, 2, …, m − 1; Ψ m = { T < t m } ; P ( Ψ k ) the occurrence probability of Ψ k and it is considered as p 0 k − p 0 k + 1 for k = 0,1, …, m − 1 and p 0 m for k = m ( 69 ); P ( F Ψ k ) represents the conditional failure probability of Ψ k .

Details of Employed Computational Model

The intrinsic spatial variability of soil parameters, uncertainty in the underlying stratigraphy, and errors in the modeling are only a few examples of the elements that could potentially influence the results of a slope stability analysis. The amount of uncertainty associated with the soil parameter can be logically incorporated into slope analysis through the employment of ML. ML is a data-driven, empirical technique where a computer program “learns” from a dataset without being explicitly programmed with a problem and solution ( 70 – 72 ). Following the implementation of MCS and SS in the spreadsheet, three neural network-based models, RNN, LSTM, and BNN, have been investigated for the optimization of the worst failure condition. The theoretical background of the developed models is discussed in the following subsection.

Recurrent Neural Network

David Rumelhart’s research from 1986 served as the foundation for the RNN, which is a specialized type of neural network ( 73 ). RNNs are neural networks that are adept at modeling sequence data to comprehend the temporal dependencies and patterns within the data, making them particularly effective for tasks such as natural language processing, speech recognition, and time-series prediction. The RNN is different from the traditional neural network, also known as feed-forward neural network, in that it has an input layer, a hidden layer, and an output layer. It has the ability to utilize previous information to influence subsequent steps. In contrast, the looping mechanism of an RNN serves as a bridge to facilitate the transfer of information between steps. This information is encapsulated in the hidden state, serving as a representation of the preceding input and enabling the model to incorporate temporal dependencies in its learning process.

Here, x_t stands for input, S_t and y_t stand for H L and outer layer outputs at time t, respectively, and W, U, V indicates shared parameters of weight of inputs, weight of inputs at current state and weight of outputs at time t, respectively, as shown in Figure 2. The modeling techniques are defined by Equations 7 and 8.

S t = f ( U x t + W S t − 1 + b h )

(7)

y t = f ( V S t + b o )

(8)

where b h and b o are the respective bias vectors for the hidden and outer layers, and f(.) is an activation function like sigmoid function. As the time in travel increases, it can be seen that the H L in the front of the RNN expands, causing information loss and error accumulation in the H L in the back ( 74 ). In practice, when the time interval is prolonged and there is a vanishing gradient problem, the performance of RNN in maintaining memory can be compromised during training using the backpropagation through time. To overcome this issue, several variations of RNN have been developed.

OPEN IN VIEWER

Figure 2.

Basic architecture of simple recurrent neural network ( 74 ).

Long Short-Term Memory

In 1997, Hochreiter and Schmidhuber introduced a novel variant of RNN called LSTM. The LSTM architecture was specifically designed to address the issue of vanishing gradients in traditional RNNs, enabling it to better handle long-term dependencies and to foretell the future events with greater accuracy ( 75 ). The LSTM architecture is designed to acquire knowledge of long-term dependencies among sequential data and time steps. Three different types of gates, namely, forget, input, and output, as well as a memory cell, are added by LSTM neurons. Additionally, it includes the sigmoid layer, the tanh layer, and the point-wise multiplication operation. The basic LSTM unit is shown in Figure 3. Notably, the σ 1 , σ 2 , and σ 3 are the three sigmoid activation functions for forget, input, and output gates, respectively, and ∅ 1 and ∅ 2 are the hyperbolic tangent activation functions (tanh) as shown in Figure 3. Mathematically, the LSTM unit can be parameterized using Equations 9 to 13.

f t = σ 1 ( W f . [ h t − 1 , x t ] + b f )

(9)

i t = σ 2 ( W i . [ h t − 1 , x t ] + b i )

(10)

c t = tan h ( W c . [ h t − 1 , x t ] + b c )

(11)

o t = σ 3 ( W o . [ h t − 1 , x t ] + b o )

(12)

h t = o t × tanh ( c t )

(13)

where x t , h t , c t , f t , i t and o t illustrate the input, final output, cell state, forget gate, input gate and output gate vector, respectively, whereas W f , b f , W i , b i , W c , b c , W o and b o are the weight and bias of the forget gate, input gate, cell state and output gate, respectively.

OPEN IN VIEWER

Figure 3.

Basic unit of long short-term memory.

Bayesian Neural Network

A BNN is a stochastic ANN that utilizes Bayesian inference in its development ( 76 – 78 ). Researchers came up with the BNN algorithm so they could take the best aspects of neural networks and probabilistic modeling and merge them ( 79 ). The purpose of this initiative is to make use of probabilistic modeling’s boundaries and guarantees, as well as the capability of neural networks to approximate universal functions. BNNs incorporate a probabilistic approach, whereas conventional neural networks use point estimates for their weights and biases. Instead, BNNs represent these parameters using posterior probability distributions ( 80 , 81 ). The universal continuous function approximation capabilities of BNNs are similar to those of conventional neural networks. The learner in the Bayesian framework draws inferences about the posterior distribution p ( w | k ) across the parameters w of the model after analyzing the data k. The posterior distribution is given by Baye’s rule: p ( w | k ) < ⋉ p ( k | w ) p ( w ) , where p ( w | k ) is the likelihood of k given by the model with parameter w, and p(w) is the prior distribution over the parameters. The prediction of the model of the new set of data is then given by Bayesian model average ( 82 ). The basic architecture of BNN is presented in Figure 4.

OPEN IN VIEWER

Figure 4.

Bayesian neural network architecture.

Illustrative Example

As an illustration, this section applies the PA of soil slope by MCS and SS with the use of MS-Excel “UPSS 3.0 Add-ins” to evaluate p f and β . As illustrated in Figure 5, the embankment has a total height of 11.693 m, which comprises (a) a 1000 mm-thick prepared sub-grade and (b) an embankment fill that reaches a height of 10.693 m. The embankment has a slope angle of 26.6°, corresponding to an inclination ratio of 2H:1V. In this study, the soil layers situated below the natural ground surface have been divided into two distinct zones: sub-soil layer-1, with a depth of 2.0 m, and sub-soil layer-2, which lies beneath sub-soil layer-1 and has a depth greater than 2.0 m. It is noted that a berm width 1.5 m has been provided at a height of 5.4 m from the top of the embankment. In accordance with the recommendations made by the Research Designs & Standards Organisation (RDSO), the width of the berm should be at least 1.5 m on either side of the embankment for every maximum height of 6.0 m (measured from the top of the embankment) ( 83 ). It has been presumed that sub-soil layer-1 is predominately made up of silty sand and clay that ranges from moderate to low compressibility. On the contrary, the lower stratum is made up of silty sand that is combined with well-graded sand and silt. Table 1 illustrates the summary of soil parameters for the intended embankment. Particularly, in the sub-grade and embankment layer, the parameters c 1 , ∅ 1 , and γ 1 are explicitly modeled using a 1-D lognormal random field, and stability of slope is examined by OMS. On the other hand, the parameters for the sub-soil layer-1 and sub-soil layer-2 are presumed to be deterministic. It should be emphasized that the horizontal seismic coefficient (k_h) and the vertical seismic coefficient (k_v) have been included in this analysis. For seismic condition, k_h were set to 0.12 for zone-III and 0.14 zone-IV. However, in each zone k_v = 0.5 × k_h is taken into consideration. This line has been planned in accordance with the zoning map of India, as recommended furnished by RDSO ( 83 ). It is noted that influence of pore water pressure ratio (viz., r u = 0 , r u = 0 . 05 , and r u = 0 . 10 ) has also been addressed in this illustrative problem.

OPEN IN VIEWER

Figure 5.

Slope geometry of 11.693 m-high embankment (all dimensions are in mm).

OPEN IN VIEWER

Table 1.

Summary of Soil Parameters

Soil parameters	Mean	Unit
Sub-grade and embankment
c 1	10.00	kN/m2
∅ 1	26.00	°
γ 1	17.00	kN/m3
Sub-soil layer-1
c 2	20.00	kN/m2
∅ 2	26.00	°
γ 2	18.70	kN/m3
Sub-soil layer-2
c 3	0.00	kN/m2
∅ 3	30.50	°
γ 3	16.70	kN/m3

Deterministic Calculation

Deterministic analysis is a strategy to predict the behavior of system under conditions where input parameters remain constant or where uncertainty is not taken into account. In the current work, a deterministic analysis of a proposed case study of soil embankment has been reported using an Excel spreadsheet. To find the F S min , OMS has been investigated, assuming circular slip surface with a defined center ( x c , y c ) and radius (r). The calculation of F S min for the given problem involved dividing the sliding mass into multiple slices, each characterized by its weight ( W i ) , segment length ( Δ l i ) , soil property (i.e., c , ∅ , and γ ), as well as an angle α i relative to the horizontal plane. Importantly, numerous trial-and-error approaches with all distinct combinations of ( x c , y c ) and r for all potential slip surface by applying Equation 1 have been adopted for finding the F S min .

Figure 6 illustrates the constructed deterministic spreadsheet, showcasing the determined value F S min = 1 . 731 for an 11.693 m-high slope embankment. It is important to note that through numerous iterations the optimized ( x c , y c ) and radius r of the proposed embankment were reported as ( x c = 27 . 0 , y c = 23 . 0 ) and r = 23 . 0 , respectively, defining a critical slip surface. It is essential to highlight that the presented worksheet has been developed under non-seismic condition ( viz . , k h = 0 . 0 and k v = 0 . 0 ) with r u = 0 . 0 , and without considering soil uncertainty. However, the mean values slope soil parameters c 1 , ∅ 1 , and γ 1 have been taken from Table 1. The output results (i.e., lowest FS) of the generated deterministic analysis are then validated with the Slope/W module in GeoStudio 2016, as can be seen in Appendix Figure A1(A). Subsequently, the results have since been validated, and the different seismic conditions (viz., k_h = 0.12 and k_h = 0.14) and r u conditions (viz., 0.0, 0.05, and 0.10) of the illustrative problem were addressed.

OPEN IN VIEWER

Figure 6.

Generated deterministic Excel sheet.

In this study, a total thirty slices were divided in a deterministic calculation for the illustrative problem. As FS is dependent on the slip surface and soil property, the soil property is the most important input variable, as indicated in the section below.

Uncertainty Modeling

The uncertainty model (UM) spreadsheet is designed to generate random variables based on user-specified parameters. These parameters include factors such as the type of distribution, the deterministic correlation length, and the correlation structure. In this particular illustration, three parameters, namely c 1 , ∅ 1 , and γ 1 , have been treated as random variables along with the depth and considered to be fully corelated. Uniform IID (independent and identically distributed) random samples are generated using the built-in “RAND ()” function in Excel, and then transformed into random samples of the desired distribution type (e.g., normal distribution or lognormal distribution). A homogeneous lognormal random field that has an exponentially falling correlation structure is used to describe the variability with depth. In this illustrative problem, the proposed embankment which is 11.693 m high was assumed to be divided into ten layers, each having a thickness 1.1693 m; a total of eleven uniformly IID random variables are required at eleven depths. Furthermore, the generated random variables are converted into standard Gaussian variables, and then PDFs using function “NORMINV ()” and “NORMDIST ()” in MS-Excel, respectively. For this particular analysis, a one-dimensional lognormal random field was created using Equation 2. The random variables used in the analysis are correlated with one another, and the correlation matrix or covariance matrix is illustrated in Figure 7 which is subjected to Cholesky factorization to create a lower triangular matrix (L). Once the “L” has been constructed, random variables are correlated using Equation 2 and reported in Figure 8 as a UM worksheet. Of note, the correlation length ( λ ) is set as 1.0 m, whereas the construction of “L” is performed using written MATLAB code. Additional details on the process of generating random samples can be found elsewhere ( 84 ). All details are comprehensively presented in the UM worksheet, including input parameters, random sample generation, the creation of a lognormal random field for c 1 , ∅ 1 , and γ 1 , and the Cholesky matrix, shown in Figure 8.

OPEN IN VIEWER

Figure 7.

Covariance matrix Excel worksheet.

OPEN IN VIEWER

Figure 8.

Uncertainty modeling Excel worksheet.

MCS-Based p_f

For the purpose of estimating p f , MCS and SS have been applied to this illustrative case. The following procedures need to be followed to accomplish the direct MCS through the UPSS dropdown toolbar: (a) click the add-ins in Excel; (b) select the direct MCS by clicking the “Simulation” UPSS dropdown menu in Excel; (c) select the no. of runs per sample (1 by default); (d) choose the sample size (1,000 by default); (e) select the standard normal IID variables from UM sheet; (f) select the system response, y (which can be defined as the reciprocal of FS, i.e., y = 1/FS) from UPM sheet; (g) to record the variable, select the “Advanced” option from the popup box; (h) to initiate the MCS simulation, click the Run option on the popup menu. Figure 9 depicts the direct MCS popup window box corresponding to general inputs and advanced options.

OPEN IN VIEWER

Figure 9.

Direct Monte Carlo simulation (MCS) popup window box: (a) general inputs and (b) advance inputs.

In this illustrative case, four different combinations of COVs (i.e., case-1, case-2, case-3, and case-4) were investigated, as shown in Table 2. For the purpose of PA, 10,000 random samples have been generated for each scenario for the direct MCS simulation with λ = 1 . 0 m in this instance. The p f and β is computed by utilizing Equations 4 and 5, respectively.

Table 3 shows the outcomes of the PA of the employed direct MCS for different values of the seismic coefficient, with the p f value showing the degree of failure probability of the illustrated problem in percentage. In each case, the values of r u were fixed to 0.0, 0.05, and 0.10, whereas the seismic coefficient was configured to be 0.12 and 0.14, and they were segmented as Set-1 and Set-2, respectively. In this context, it is important to note that the p f of the soil slope, determined based on the target FS value for seismic conditions, was less than 1.1 ( 51 ).

OPEN IN VIEWER

Table 3.

Results of Probabilistic Analysis (Monte Carlo Simulation) in Relation to p f and β

Pore pressure ratio, r u	Cases	Number of samples generated	Number of samples having factor of safety < 1.1	p f ( % )	Reliability, β
Set-1
0.0	Case-1	10000	0	0.00	NA
	Case-2	10000	0	0.00	NA
	Case-3	10000	8	0.08	3.15
	Case-4	10000	92	0.92	2.35
0.05	Case-1	10000	0	0.00	NA
	Case-2	10000	6	0.06	3.24
	Case-3	10000	118	1.18	2.17
	Case-4	10000	544	5.44	1.60
0.10	Case-1	10000	1	0.01	3.72
	Case-2	10000	461	4.61	1.68
	Case-3	10000	1126	11.26	1.21
	Case-4	10000	1812	18.12	0.91
Set-2
0.0	Case-1	10000	0	0.00	NA
	Case-2	10000	0	0.00	NA
	Case-3	10000	84	0.84	2.39
	Case-4	10000	379	3.79	1.78
0.05	Case-1	10000	0	0.00	NA
	Case-2	10000	167	1.67	2.13
	Case-3	10000	752	7.52	1.44
	Case-4	10000	1424	14.24	1.07
0.10	Case-1	10000	606	6.06	1.55
	Case-2	10000	2137	21.37	0.79
	Case-3	10000	2927	29.27	0.55
	Case-4	10000	3571	35.71	0.37

Note: NA = not available.

Table 3 presents the results obtained from the MCS of the proposed illustrative slope embankment. Set-1 (i.e., k_h = 0.12 and k_v = 0.06) p f lies within the range of 0% to 18.12%, and Set-2 (i.e., k_h = 0.14 and k_v = 0.07) lies between 0% and 35.71%. This indicates that at higher seismic values, the chances of failure increase. Furthermore, it can be observed from Table 3 that higher β values correspond to lower p f values. It is observed that, p f value increases significantly from 0% to 0.92% for r u = 0 . 0 , 0% to 1.18% for r u = 0 . 05 , and 0.01% to 18.12% for r u = 0 . 10 when k_h was set to 0.12. In a similar fashion, p f value increases from 0% to 3.79% for r u = 0 . 0 , 0% to 14.24% for r u = 0 . 05 , and 6.06% to 35.71% for r u = 0 . 10 when k_h was set to 0.14. This is primarily because of the execution of a higher COV value and r u . Furthermore, a dramatic increase in p f value was observed when comparing Case-1 and Case-4 of Set-1 for r u = 0 . 10 , where the p f value increased from 0.01 to 18.12 (about 1,812 times). These findings suggest that an increase in COV leads to higher p f and lower β. In addition, the complementary cumulative distribution function (CCDF) plot for Set-1 and Set-2 with different COVs and the various r u conditions (i.e., r u = 0 . 0 , r u = 0 . 05 and r u = 0 . 10 ) has been given separately in the Appendix Figure A1(B) and Figure A1(C). This CCDF plot shows a relationship between P(Y > 1/FS) and driving variable (i.e., y = 1/FS). Notably, slope failure occurs when the driving variable y > 0.91, considering the seismic effect. It can be observed from the trend of the CCDF plot that when the COV value of the soil parameter increases, the likelihood of failure increases dramatically. Additionally, it can also be deduced from Figures A1(B) and A1(C), the p f increases correspond to an increase in seismic value.

SS-Based p_f

The SS is an advanced version of direct MCS, as has already been explained. SS was implemented for the same illustrative problem using the “Add-ins” option in Excel. The SS is performed by clicking the “simulation” button in the UPSS tool bar and choosing the “Subset simulation” option from the dropdown menu. Figure 10 depicts the general inputs as well as the advanced inputs section for the SS popup window, in which general inputs include the following seven keys: the number of runs (which has a default value of 1); the number of samples per level (N by default 200); the conditional probability ( p 0 which has a default value of 0.1); the number of simulation levels (m); the random variable; the PDF of the random variable; and the driving variable (Y, which is defined as the reciprocal of FS), whereas advance inputs options are for recording the variables.

OPEN IN VIEWER

Figure 10.

Subset simulation popup window: (a) general inputs and (b) advance options.

In this investigation, SS is performed in each case by following conditions: p 0  = 0.1, N = 500, m = 3. Keep in mind that for each simulation run, the total number of samples is m × N − ( m − 1 ) p 0 N . This gives a total of 3 × 500 − ( 3 − 1 ) × 0 . 1 × 500 = 1400 samples, which include 450 samples at level 0 and 1, respectively, and 500 samples at level 2. The p f and β can be determined through the use of Equations 6 and 5, respectively.

Referring to Table 4, the p f values for Set-1 and Set-2 fall somewhere in the range of (0%–18.2%) and (0%–35.8%), respectively. No failure is detected, when the value of r u is set to 0.0 in both sets of Case-1 and Case-2, and to 0.05 in both sets of Case-1. It is clearly seen from Table 4 for Set-1, the p f for r u = 0 . 0 lies in between 0% and 1%, and for r u = 0 . 05 lies in between 0% and 8.6% and for r u = 0 . 10 lies in between 0.02% and 18.2%. Furthermore, when k_h was set to 0.14, the range of p f values were found to be 0%–4.36% for r u = 0 . 0 , 0%–14.4% for r u = 0 . 05 and 7.72%–35.8% for r u = 0 . 10 .

OPEN IN VIEWER

Table 4.

Results of Probabilistic Analysis (Subset Simulation) in Relation to p f and β

Cases	r u	Number of samples generated	Number of samples having factor of safety < 1.1			Value of p f and β
			Level 0	Level 1	Level 2	p f ( % )	Reliability, β
Set-1
0.0	Case-1	1400	0	0	0	0.00	NA
	Case-2	1400	0	0	0	0.00	NA
	Case-3	1400	0	0	70	0.14	2.99
	Case-4	1400	0	0	500	1.00	2.33
0.05	Case-1	1400	0	0	0	0.00	NA
	Case-2	1400	0	0	30	0.06	3.24
	Case-3	1400	0	18	500	1.36	2.21
	Case-4	1400	0	380	500	8.60	1.36
0.10	Case-1	1400	0	0	10	0.02	3.54
	Case-2	1400	0	219	500	5.38	1.61
	Case-3	1400	8	450	500	11.60	1.20
	Case-4	1400	41	450	500	18.20	0.90
Set-2
0.0	Case-1	1400	0	0	0	0.00	NA
	Case-2	1400	0	0	0	0.00	NA
	Case-3	1400	0	14	500	1.28	2.23
	Case-4	1400	0	168	500	4.36	1.71
0.05	Case-1	1400	0	0	0	0.00	NA
	Case-2	1400	0	65	500	2.30	2.00
	Case-3	1400	0	326	500	7.52	1.44
	Case-4	1400	22	450	500	14.40	1.06
0.10	Case-1	1400	0	336	500	7.72	1.43
	Case-2	1400	63	450	500	22.60	0.75
	Case-3	1400	98	450	500	29.60	0.54
	Case-4	1400	129	450	500	35.80	0.36

Note: NA = not available.

Figure 11 illustrates the relationship between p f (%) and COVs for different values of r u (i.e., 0.0, 0.05, and 0.10) using SS at two distinct seismic values (i.e., k_h = 0.12 and k_h = 0.14). It is evident from the figure that as r u values increase at both seismic values, p f increases sharply. Moreover, for a specific combination of r u and seismic values, p f increases significantly with an increase in the COVs of soil parameter. For instance, in Case-1, the value of p f was found to be 0.02 for r u  = 0.10 and k_h = 0.12, whereas for the same case, p f increased to 7.72 for r u  = 0.10 and k_h = 0.14. This observation suggests that the increase in p f is caused by the rise in seismic effect.

OPEN IN VIEWER

Figure 11.

Trends between p f (%) and COVs at different r u condition: (a) k_h = 0.12 and (b) k_h = 0.14.

It is evident from Tables 3 and 4 that the results obtained through SS are nearly comparable to those obtained through MCS. However, the p f values obtained by SS are either equal to or higher than those obtained through direct MCS. It can be inferred that the same level of accuracy (in relation to p f ) can be achieved with SS using a smaller number of random samples (i.e., 1,400) compared with MCS, which generates 10,000 random samples. Additionally, SS takes less computation time (i.e., 16.37 s) compared with direct MCS (which takes approximately 1 min) to achieve the same degree of p f . The Appendix displays the CCDF plot for p f versus P(Y > 1/FS) for the illustrative slope using SS in Figures A1(D) and A1(E), respectively. A comparison of the CCDF plots obtained from the SS and MCS analyses indicates that the SS analysis can offer lower p f values. MCS reported a lower p f value of up to 10 − 4 , whereas the SS analysis provides values up to 10 − 6 .

Referring to Table 4, the p f of Set-2 and the coefficient of r u  = 0.10 for Case-2 have been selected as an illustrative example for estimating the p f using SS. Table 5. summarizes the procedure for estimating the p f of an instructive example. In a total of 1,400 SS samples, 1,013 were found to have FS less than 1.1, making them failure samples. These 1,013 samples are broken down as follows: sixty-three samples represent Level 0 (i.e., sixty-three samples failed in Level 0), 450 represent Level 1(i.e., all samples failed in Level 1), and 500 represent Level 2 (i.e., all sample failed in Level 2) of the simulation. With the use of these failure samples, P ( F Ψ k ) for k = 0, 1, 2 is estimated as (63/450), (450/450), and (500/500) (see column 5). In addition, P ( Ψ k ) for k = 0, 1 and 2, the values are 0.9, 0.09, and 0.01 (see column 2), respectively. The p f can be estimated as ( 63 450 ) × 0 . 9 + ( 450 450 ) × 0 . 09 + ( 500 500 ) × 0 . 01 = 0 . 226 or 22.6% that has been calculated by using Equation 6. Based on the resulting p f , β can be computed as ϕ − 1 ( 22 . 6 % ) which is equal to 0.75.

OPEN IN VIEWER

Table 5.

Failure Probability from Subset Simulation of Instructive Example

Simulation level	P ( Ψ k )	Number of samples in bins, Ψ k	Number of samples having FS < 1.1	P ( F Ψ k )	P ( F Ψ k ) P ( Ψ k )	p f = ∑ k = 0 m P ( F Ψ k ) P ( Ψ k )
0	0.9	450	63	0.14	0.126	0.226
1	0.09	450	450	1	0.09
2	0.01	500	500	1	0.01

Computational Modeling

Following the completion of a PA using MCS and SS, the worst p f combination is chosen and subsequently optimized using ML. Remarkably, p f  = 35.8% of SS has been proven to be the poorest possible combination (viz., Set-2, Case-4, and r u = 0 . 10 ). The run sheet of the worst combination of SS with an input random variable (i.e., c_11.693 to c_0, ∅ _11.693 to ∅ _0 and γ _11.693 to γ _ 0 ) and FS are chosen to generate a total of 1,400 random data by using written MS visual basic code (VBA). Following the generation of 1,400 datasets at random, the employed models, namely, RNN, LSTM, and BNN were constructed by using all of the input and output parameters, followed by the random bifurcation of the training and testing datasets. To accomplish this, a total of 1,400 datasets are selected at random, normalized, and then divided in the following ratio: 70% (i.e., 980) are used for the training phase, and the remaining 30% (i.e., 420) are used for the testing phase. Of note, the bifurcation of datasets was done using written MATLAB code. It is also noted that the training datasets are used to create three neural network-based models, while the testing datasets are used to assess the model’s prediction precision in estimating the FS of an embankment that is 11.693 m high.

Following computational modeling, multiple evaluation metrics, namely determination coefficient (R²) ( 89 ), root mean square error (RMSE), Willmott’s index of agreement (WI), Legate and McCabes’s index (LMI), and Theli’s inequality coefficient (TIC) were determined, the mathematical expressions of which are given in Equations 14 to 18.

R 2 = ∑ i = 1 N ( y i − y avg ) 2 − ∑ i = 1 N ( y i − y i pre ) 2 ∑ i = 1 N ( y i − y avg ) 2

(14)

RMSE = 1 n ∑ i = 1 n ( y i − y ¯ i ) 2

(15)

WI = 1 − [ ∑ i = 1 n ( y i − y ¯ i ) 2 ∑ i = 1 n { | y ¯ i − y avg | + | y i − y avg | } 2 ]

(16)

LMI = 1 − [ ∑ i = 1 N | y i − y i pre ∑ i = 1 N | y i − y ¯ i pre ]

(17)

TIC = 1 N ∑ i = 1 N ( y i pre − y i ) 2 1 N ∑ i = 1 N y i pre 2 + 1 N ∑ i = 1 N y i 2

(18)

where N stands for the number of datasets, y avg and y ¯ i pre is the mean of the actual and anticipated slope safety factor, respectively. y i and y i pre represent the actual and modeled predicted values of the ith values. It is absolutely necessary to have a solid comprehension that the values of these indices should match to their ideal values; details are already mentioned in the literature ( 90 – 92 ). Notably, conventionally evolution metrices should have the value 1 for R², 0 for RMSE, 1 for WI, 1 for LMI and 0 for TIC ( 51 , 92–95).

Performance Assessment

In this section, the performance of the employed neural network-based model is presented and discussed. The construction of the model involved utilizing a total of 980 datasets, as described above. The number of hidden neurons ( N h ) was studied for the RNN utilizing a trial-and-error methodology throughout a range that extended from 1 to 40; the optimal value of N h was found to be 20. The optimized RNN architecture consists of a hidden layer with twenty nodes, which is a simple RNN. The input shape has thirty-three time steps and one input feature, and the activation function used is ReLU (rectified linear unit). This is followed by two dense layers and an output layer. In addition, the “mean square error” function was used for best fitting, and the “Adam” function was used for optimization in the optimized RNN architecture. Similarly, in the LSTM neural network, a random hit-and-trial approach was used to investigate the first layer with various input nodes to achieve optimal performance. The optimized LSTM model consists of one input layer with thirty-three nodes, a LSTM hidden layer with fifteen neurons, an input shape of (33,1), and an activation function “ReLU”. This is followed by a dense layer and an output layer. Additionally, the “mean square error” function was used for best fit, and the “Adam” function was used for model optimization. It should be noted that in both the RNN and LSTM neural networks, 500 epochs and a batch_size of sixteen were used. For constructing the BNN model, first layer with thirty-three input parameters, one hidden layer with ten N h and an output layer were used. Subsequently, for best fit of the model, fifteen batch_size and 500 epochs were considered. It should be noted that for each model, different epochs between 300 and 500 were trialed, and the optimum epoch was found to be 500.

The prediction outcomes of the ML model were built for optimization of worst condition (i.e., p f = 35.8%) for estimating the FS. Here, evaluation metrics for the proposed model are illustrated in Table 6.

OPEN IN VIEWER

Table 6.

Evaluation Metrics of Training Datasets

Model	R 2	RMSE	WI	LMI	TIC
RNN	0.9897	0.0092	0.9973	0.8940	0.0041
LSTM	0.9953	0.0062	0.9988	0.9314	0.0027
BNN	0.9510	0.0201	0.9876	0.7870	0.0089

Note: RNN = recurrent neural network; LSTM = long short-term memory; BNN = Bayesian neural network; RMSE = root mean square error; WI = Willmott’s index of agreement; LMI = Legate and McCabes’s index; TIC = Theli’s inequality coefficient.

Of note, evaluation metrics have been employed to determine the model’s performance. From R², the LSTM is 0.9953, which is higher than RNN (0.9857) and BNN (0.9510), showing better performance. From the RMSE and TIC, LSTM is 0.0062 and 0.0027 on the training phase, respectively, which is lower than 0.0092 and 0.0041 for the RNN, and 0.0201 and 0.0089 for the BNN model, respectively. Following a comparative analysis, it was concluded that the LSTM model demonstrated the most superior performance during the training phase. Additionally, it can be demonstrated that all the employed models obtained more than 95% accuracy (based on the R² value), demonstrating a good fit with the datasets that have been obtained. However, to provide a clearer illustration of the fit, the marginal histogram and error plot of the actual and predicted FS values for each of the models that were used are presented in Figures 12 and 13, respectively.

OPEN IN VIEWER

Figure 12.

Marginal histogram for training dataset: (a) RNN, (b) LSTM, and (c) BNN.

OPEN IN VIEWER

Figure 13.

Error plot of (a) RNN, (b) LSTM and (c) BNN (Training phase).

Following the building of the models, the unused datasets (a total of 420) were utilized during the testing phase of the models to determine their potential. Table 7 provides an instance of the evaluation metrics of the proposed-built models for testing datasets, and it can be concluded that throughout the testing phase, the suggested LSTM had the best precision in relation to R² = 0.9962 and RMSE = 0.0051. In contrast, during the testing phase, the BNN model was shown to be the least effective across all evaluation metrics. However, the marginal histogram and error plot of the actual and predicted FS values in the testing phase for each of the models are presented in Figures 14 and 15.

OPEN IN VIEWER

Table 7.

Evaluation Metrics of Testing Datasets

Model	R 2	RMSE	WI	LMI	TIC
RNN	0.9923	0.0072	1.0000	0.9946	0.0032
LSTM	0.9962	0.0051	1.0000	0.9964	0.0022
BNN	0.9488	0.0186	0.9998	0.9871	0.0083

Note: RNN = recurrent neural network; LSTM = long short-term memory; BNN = Bayesian neural network; RMSE = root mean square error; WI = Willmott’s index of agreement; LMI = Legate and McCabes’s index; TIC = Theli’s inequality coefficient.

OPEN IN VIEWER

Figure 14.

Marginal histogram for testing dataset: (a) RNN, (b) LSTM, and (c) BNN.

OPEN IN VIEWER

Figure 15.

Error plot of (a) RNN, (b) LSTM and (c) BNN (Testing phase).