Exploring the Potential of Machine Learning in Stochastic Reliability Modelling for Reinforced Soil Foundations

Abstract

This study introduces a novel application of gene expression programming (GEP) for the reliability analysis (RA) of reinforced soil foundations (RSFs) based on settlement criteria, addressing a critical gap in sustainable construction practices. Based on the principles of probability and statistics, the soil uncertainties were mapped using the first-order second-moment (FOSM) approach. The historical data generated via a parametric study on a validated finite element numerical model were used to train and validate the GEP models. Among the ten developed GEP frameworks, the best-performing model, abbreviated as GEP-M9 ($R^2$ = 0.961 and RMSE = 0.049), in the testing phase was used to perform the RA of an RSF. This model’s effectiveness in RA was affirmed through a comprehensive evaluation, including parametric sensitivity analysis and validation against two independent case studies. The reliability index (β) and probability of failure (P_f) were determined across various coefficient of variation (COV) configurations, underscoring the model’s potential in civil engineering risk analysis. The newly developed GEP model has shown considerable potential for analyzing civil engineering construction risk, as shown by the experimental results of varying settlement values.

1. Introduction

The use of geosynthetic reinforcements is becoming popular for constructing safe and sustainable reinforced soil foundations (RSFs), especially for low-rise buildings. In order to ensure the safe design of RSFs, determining their settlement under a service load is of paramount importance for practitioners. In the literature, numerous analytical and in situ test (standard penetration, cone penetration, dilatometer, and plate load test)-based methods are available for predicting the immediate settlement of unreinforced sandy soil foundations. Lutenegger and DeGroot [1] rigorously reviewed and summarized the details of all the well-known methods for evaluating the settlement of shallow unreinforced foundations on cohesionless soils. However, the settlement estimation of RSFs require further investigations.

Since Binquet and Lee [2] investigated the behavior of soil reinforced with metal strips, numerous researchers have carried out experimental, numerical, and analytical works to study the load-settlement behavior of the footings resting on reinforced soil bed [3,4,5,6,7,8,9,10]. Omar et al. [6] studied the increase in the bearing capacity ratio (BCR) of the square and the strip footing resting on sand reinforced with geogrid layers. They reported that the maximum BCR occurred when the geogrid layers were placed at a depth of approximately 1.4 and 2 times the width of footing (B) for square and strip footings, respectively. Adams and Collin [7] investigated the behavior of geosynthetic-reinforced soil foundations on granular soil by conducting large-scale footing load tests. They attributed the increase in the BCR of the reinforced bed to the increase in the shear strength of the composite (sand–geosynthetic system) due to the inclusion of geosynthetic layers. Chen [8] carried out a comprehensive experimental study on geosynthetic-reinforced soil foundations using two types of geosynthetics, namely, geotextiles and geogrids. The results showed that the inclusion of geosynthetic layers helps in spreading the load over a wider area, thus significantly enhancing the bearing capacity and decreasing the settlement of RSFs under a particular load. Moreover, they also concluded that geogrid-reinforced foundations performed better than geotextile-reinforced foundations in settlement control tests. Abu-Farsakh et al. [4] investigated the behavior of a geogrid-reinforced soil foundation on sand using small-scale laboratory tests. They observed a decrease in the immediate settlement of approximately 20%, for all footing pressure levels, when the footings were reinforced with two or more geogrid layers. Mosallanezhad et al. [11] performed several footing load tests to investigate and compare the performance of simple geogrids (four layers) and grid-anchor reinforcing systems. They observed increases in the BCR of the geogrid and grid-anchor systems by 1.65 and 2.98, respectively. All this experimental work comprehensively proves the potential benefits of constructing RSFs. However, there is a need to develop a method for predicting the settlement of RSFs under service loading conditions in a realistic and rational way.

Uncertainty is a common factor in settlement prediction, as in many geotechnical engineering problems. Although an accurate settlement prediction is crucial for design, such uncertainty can lead to erroneous estimates of settlement. For an unreinforced soil foundation, the uncertainty in settlement prediction is directly linked to the uncertainties associated with the input variables and the model employed for settlement prediction [12,13]. Similarly, the settlement estimation of RSFs is also affected by those uncertainties.

Utilizing probability theory and statistical methods, geotechnical evaluations can effectively incorporate uncertainties related to soil properties [14]. In this approach, the performance of geotechnical structures is statistically assessed, often through the use of the reliability index (β) and/or the probability of failure (POF). The latter is defined as the chance that the structure will fail to meet established performance standards [14]. Over the years, various methodologies like direct simulation Monte Carlo (DSMC), the first-order reliability method (FORM), and first-order second-moment (FOSM) have been devised to conduct the reliability analysis (RA) of geotechnical structures. These methods analyze probabilistic estimates of soil features and sub-soil layers to calculate the β and/or POF for specific geotechnical design scenarios. Historically, such techniques have been utilized to analyze the RA concerning the settlement of shallow foundations [14,15,16].

In recent times, machine-learning (ML) techniques have become prevalent for solving engineering problems, demonstrating their wide-ranging applicability [17,18,19,20]. Capitalizing on these technological advances, there is a growing interest in using artificial intelligence (AI) to refine geotechnical reliability analysis (GRA) with ML, aiming at boosting the precision and efficiency of computations [21]. The core principle behind integrating ML into GRA involves learning from data, which mainly comprises random variable inputs or samples from random fields of geomaterial characteristics (like friction angle, cohesion, standard penetration test (SPT) value, etc.), along with the outcome of interest, often derived from geotechnical software (e.g., PLAXIS, GeoStudio, etc.). After sufficient training and validation, the ML-powered reliability analysis models achieve their intended performance levels, enabling them to predict the failure probability of geotechnical structures reliably and efficiently [21,22]. A review of the existing literature indicates that this innovative approach has been successfully implemented in evaluating the reliability of various geotechnical engineering challenges, such as the bearing capacity and settlement of shallow foundations [12,23,24,25,26]. However, to date, the RA of RSFs, in terms of the settlement prediction, has not received any attention in the literature.

Numerous studies have leveraged ML models like artificial neural networks (ANNs), support vector machines (SVMs), evolutionary polynomial regression (EPR), neuro-fuzzy systems, and hybrid ANNs to predict the settlement of unreinforced soil foundations [27,28,29,30,31,32]. However, the research related to the application of ML techniques in RSFs is very limited. A handful of researchers have employed data-driven methods to investigate the behavior of RSF [33,34,35,36,37,38,39,40]. Soleimanbeigi and Hataf [33,34] proposed a back-propagating neural network for predicting the ultimate bearing capacity and settlement at peak footing loads of RSFs. Similarly, Raja and Shukla [35,38] applied extreme-learning machines (ELMs) and multivariate adaptive regression splines (MARSs) to predict the load-bearing capacity and settlement of RSFs, demonstrating that their models could surpass many conventional methods in accuracy. Their proposed models successfully predicted the ultimate bearing capacity and settlement, and have outperformed several traditional approaches. A few researchers have employed non-linear regression and a hybrid ANN model to predict the settlement of geosynthetic-reinforced soil foundations [37,41]. These studies showed that the ML-based models performed well in predicting the settlement of footing resting on soil reinforced with geosynthetic layers.

However, there are certain limitations associated with these models, for example, the black-box nature of the various ML models such as the ANN or hybrid ANN, SVM, and ELM. Moreover, the local minima trap issue in ANNs may produce unfavourable outcomes during the validation process [42,43]. Despite these challenges, a comprehensive review of the existing literature reveals a gap: to date, there has been no effort to apply ML-based models specifically for the reliability analysis (RA) of an RSF.

This research aims at filling this research gap by proposing a gene expression programming (GEP)-based model to perform the RA of RSFs. The main advantage of GEP over other ML-based models is its ability to produce closed-form mathematical expressions for predicting the outcome. By utilizing a validated finite element model database, several GEP models were trained and tested, and the most optimal GEP model has been selected to perform the RA of the RSF.

The remainder of this manuscript is organized in the following manner: Section 2 outlines the methodology employed, encompassing the reliability analysis and the use of GEP. The experimental approach and numerical analysis, along with data processing and analytical techniques, are detailed in Section 3 and Section 4, respectively. Section 5 delves into the results and discussion, featuring an evaluation of performance metrics, an examination of parametric sensitivity, and the execution of independent validation procedures. The conclusions are presented in the final section of the manuscript.

2. Methodology

This section focuses on the development of GEP-based intelligent models for predicting the settlement of RSF, as well as incorporation of the RA. The theoretical background of RA and a brief overview of GEP are included in this section.

2.1. Reliability Analysis

Geotechnical analysis and design can logically account for uncertainties (such as the natural variability of soils) in soil parameters with the help of probability theory and statistics. The reliability of a geotechnical structure is defined as the extent to which it meets its design objectives over a specified time frame and under specified loads. In other words, it is the likelihood that the structure will not reach the limit over a specified time frame. One of the first steps in performing RA is defining the connections between the various inputs and outputs. This study uses a serviceability limit state criterion to assess the RSF’s dependability, and the FOSM for RA. FOSM estimates the mean and variance of the performance function using the first terms of a Taylor series expansion of the function. This section provides a short overview of the aforementioned procedure.

Let d_l represent the demand (anticipated loadings) placed on an engineering system, and c_r represents its capacity (available resistance). As d_l and c_r are unknown, the connected variables can be characterized by the mean, variance, and correlation. A limit state function, also known as a performance function (c_r d_l), can be used to define the system’s safety margin as $M_s=\left(c_r-d_l\right)$. Accordingly, the probability of failure, $P_f$, can be expressed as follows:

$$P_f=P\left[\left(c_r-d_l\right)\le 0\right]$$

(1)

Considering the above equations, the reliability index ($\beta$) and P_f are defined as follows:

$$\beta =\frac{{\mu }_c-{\mu }_d}{\sqrt{\left({\sigma }_c^2+{\sigma }_d^2\right)}}$$

(2)

$$P_f=1-\psi \left[\frac{{\mu }_c-{\mu }_d}{\sqrt{\left({\sigma }_c^2+{\sigma }_d^2\right)}}\right]$$

(3)

where ${\mu }_c$ and ${\mu }_d$ are means of the capacity and demand, respectively; ${\sigma }_c$ and ${\sigma }_d$ are the standard deviations of the capacity and demand, respectively; while $\psi$is the cumulative distribution function (standard normal) represented as N (0,1).

2.2. Overview of GEP

In comparison to more conventional methods like regression analysis and neural networks, GEP, a relatively recent supervised machine-learning technique, has shown to provide more accurate modeling of experimental outcomes [44]. Like genetic programming (GP) and genetic algorithm (GA), GEP is an evolutionary algorithm that uses samples of population, selects fitness from them, and includes genetic variability through the use of one or more genetic operators, all of which are important to the programming of GEP. It combines the linear, fixed-long chromosomes of genetic algorithms with the expression parses of several genetic-programming languages [45]. The main parts of GEP are the chromosomes and the expression trees (ETs). Each chromosome has a set of instructions, or genes, that are translated into expression through the use of a specialized language called Karva [44]. Genes encoding GEPs are all the same length, but the ETs they code vary in sizes and forms [46]. Figure 1 represents the typical ET in GEP. Connection functions keep these independent genes working together to build a chromosome. Typically, genes are split between a “head” (encoded function) and a “tail” (non-encoded function). Variables, functions, and constants all point to the encoded function, whereas only constants point to the non-encoded one [47]. In case where the terminal symbols in the gene’s head are insufficient to encode a function, the constants and variables in the gene’s tail are sometimes employed as supplemental terminal symbols [48]. The head may contain functions such as addition (+), subtraction (−), multiplication (×), and division (/), or any other mathematical operators (e.g., e, log, ln, etc.) to connect the data and form the mathematical expression.

The GEP initiates with the generation of random population of chromosomes complying with the Karva language. As described earlier, these chromosomes are represented by ETs of varying shapes and sizes. Next, the primary genetic operations such as mutation, crossover, transposition, and recombination are performed according to the ratios set by GEP modeler [48]. Figure 2a,b represents the process of crossover and mutation in GEP, respectively. The iterative process continues until the stopping criteria (maximum generations or optimum solution) are achieved. For more theoretical details regarding the GEP, readers may refer to the comprehensive studies available in literature [46,49].

3. Experimental Procedure and Numerical Analysis

For model development and probabilistic assessment, a database consisting of 475 numerical simulations conducted by Raja and Shukla [37] has been utilized in this study. The simulations were carried out using PLAXIS-3D by validating the Adams and Collin [7] large-scale footing load tests on geosynthetic-reinforced soil foundations. A brief detail of the experimental setup, numerical modeling, and data collection is presented in this section.

3.1. Experimental Setup

Adams and Collin [7] carried out extensive tests on model footings under static load to investigate the performance of isolated square footings placed on granular soil that was reinforced with geosynthetics. These tests took place at the Turner–Fairbank Highway Research Centre (T-FHRC) in the USA. The testing pit had dimensions of 7 m in length, 5.4 m in width, and 6 m in depth. They conducted a series of thirty-four experiments on square footings of varying sizes positioned on a bed of soil reinforced with three layers of geogrid. The application of force was achieved through a hydraulic jack, and a hand pump was employed to maintain steady pressure. The force applied was measured using a load cell and a strain indicator box, while the resulting settlements were accurately gauged and recorded by four linear variable displacement transducers (LVDTs), with a data acquisition system capturing the measurements.

3.2. Numerical Modelling

Numerical modeling was conducted by utilizing the PLAXIS-3D program. In order to verify the model accuracy, the Adams and Collin [7] square footing (0.91 m × 0.91 m) experimental results were compared with numerical model predictions for both unreinforced and reinforced soil. A brief explanation is provided herein, whereas a detailed description of the numerical simulations can be found elsewhere [37]. A rigid foundation was used in the modeling process. The hardening soil model (HSM) was utilized to depict the soil structure, enabling the simulation of the continuum elements and interfaces of the soil. The model represented the soil and its interfaces through ten-noded tetrahedral 3D finite elements. Additionally, the geogrids were modeled as linear elastic materials from the PLAXIS-3D material library. A conducted sensitivity analysis indicated that the coarse mesh with refinements sufficed for these simulations. As illustrated in Figure 3, there is a good correlation between the numerically derived stress–strain relationships and the experimental findings, with strain expressed as the ratio of the settlement to the footing width, s/B.

The reliability of the established numerical model for simulating and forecasting the settlement response of RSFs can be concurred from the results shown in Figure 3. After confirming the accuracy of the model, a comprehensive database was compiled, featuring results from 475 full-scale simulated model footing load tests. Key factors influencing RSF settlement, such as the footing’s width and shape (B and L/B), the magnitude of the applied load (q), the depth of the soil bed (d), soil strength parameters (cohesion (c) and angle of internal friction (ϕ)), characteristics of the reinforcement including the uniform reinforcement spacing ratio (z/B), the number of reinforcement layers (N), and the initial tensile modulus (J), were identified as significant. The previous studies showed that all these inputs are well recognized in investigating the behavior of RSFs [37,41].

4. Data Processing and Analysis

4.1. Data Processing

Data processing is an essential step in the ML process as it prepares the data for analysis and modeling. Data processing involves cleaning, transforming, and normalizing the data to ensure that it is in a format that can be easily understood and analyzed by the ML algorithm. Additionally, data processing can also include feature selection and engineering, which helps to identify and extract the most relevant and informative features from the data to improve the model performance. Overall, data processing is important for the success of any ML project, since it helps to guarantee that the data are of sufficient quality and are prepared for analysis.

Table 1. Descriptive information of the employed dataset.

Parameters	q	ϕ	c	K	z/B	B	L/B	N	d	s
	(kN/m2)	(°)	(kN/m2)	kN/m	(-)	(m)	(-)	(nos.)	(m)	(mm)
Minimum	50.00	30.00	0.00	500.00	0.10	1.00	1.00	1.00	0.80	0.01
Mean	297.37	39.67	2.41	1237.74	0.23	1.47	2.44	3.71	1.35	2.85
Median	300.00	40.00	1.00	1000.00	0.20	1.00	2.00	3.00	1.20	1.93
Mode	50.00	40.00	1.00	1000.00	0.30	1.00	2.00	3.00	1.20	0.78
Maximum	600.00	50.00	10.00	3000.00	0.40	3.00	10.00	7.00	2.40	33.90
Standard Error	8.92	0.21	0.13	24.23	0.00	0.03	0.08	0.06	0.01	0.16
Kurtosis	−1.34	0.24	1.70	1.20	−0.90	0.01	5.90	−0.05	1.71	20.40
Skewness	0.23	0.24	1.66	1.28	0.00	1.16	2.30	0.75	1.60	3.66

presents the descriptive statistics of the complete dataset used in this study. It is important to check for multicollinearity when building ML models because it ensures that the independent variables used in the model are not highly correlated with one another [50]. This can lead to unstable and unreliable model estimates, as well as difficulty in interpreting the importance of each variable in the model. The review of the literature suggests that, if the correlation coefficient (r) is greater than 0.7–0.8, the model is highly affected by the collinearity issue [51,52]. Figure 4 shows the correlation matrix which indicates that the relationship between any two dependent variables does not go beyond 0.7.

4.2. Computational Modeling for Probabilistic Analysis

Figure 5 represents the flowchart of the complete methodology used in the development of GEP models for a probabilistic analysis of the settlement of the RSF. The complete process is segregated into three distinct stages. In stage 1, several combinations of soil and geosynthetic parameters (input parameters) were used to determine the settlement of the RSF. Each combination represents one data point that was fed to PLAXIS 3D for the estimation of the settlement of the RSF [45]. In this way, the obtained 475 points were partitioned into the training and testing dataset. It is noteworthy that, using the holdout validation technique, several data splitting ratios, such as 70:30, 75:25, and 80:20, were applied. The most effective results were achieved by allocating 80% of the data for training purposes and the remaining 20% for testing the models. The training dataset was utilized to develop the GEP models as outlined in Section 2, while the testing dataset served to evaluate the models’ capability to generalize. It is important to note that all data were normalized based on the range of each variable (minimum–maximum) within the training dataset before being input into the GEP.

During the second stage, the performance of the GEP models was examined using a variety of statistical indicators. The ideal values for these indicators guided the assessment of each model’s accuracy. The model demonstrating the best performance on the testing data, in predicting the settlement of RSFs, was then chosen. Additionally, both a parametric study and independent validation were carried out to verify the GEP model’s rationale. For the parametric study, a one-at-a-time parametric sensitivity analysis was performed. This analysis aimed to assess the model’s performance in line with the established understanding in the literature regarding the behavior of RSFs. This process led to the selection of the most suitable model for predicting the RSF settlement.

In the third stage, the best GEP model was used in performing the RA of RSFs under service loading conditions. This was accomplished by generating various sets of input parameters with varying levels of the non-deterministic input parameters’ coefficient of variation (COV). Then, the new dataset is normalized using the actual dataset’s (training) minimum and maximum values. Thereafter, the best GEP model is used to predict the settlement values followed by the probabilistic analysis (calculation of $\beta$ and P_f) of the settlement estimation.

It is important to note that the holdout validation technique, parametric sensitivity analysis, and independent validation, as implemented in stages 1 and 2, respectively, also act as tools for controlling overfitting. These techniques have been commonly employed in numerous previous studies to ensure models generalize well to new, unseen data [53,54,55,56].

5. Results and Discussion

For developing GEP models, various combinations of chromosomes, genes, and head sizes were applied, and, for each combination, the best model was selected to predict the settlement of RSFs. Based on the previous research on GEP-based predictive modeling and the hit and trial approach, the quantity of chromosomes, genes, and head sizes were varied between 50–200, 3–5, and 8–10, respectively. For each combination, addition was used as the linking function, and [+, ×, −, +, /, log, ln, e, x²] as the function set. Each model was run for 10,000 generations, and the optimal model was selected based on the cost functions (RMSE and R²). The parametric framework of each GEP-based model employed to predict the settlement of the RSF is given in the Supplementary File (Table S1). It may be noted that ten configurations were applied (F1…, F10) and the respective models were abbreviated as GEP1…, GEP10. By utilizing the model configurations, the ETs were obtained from the Karva programming and translated into simple mathematical formulae to predict the settlement of geosynthetic-reinforced soil foundations (see Table S2). The best GEP model with the highest performance was obtained by Model 9, expressed as GEP-M9 with R² of 0.946 (average training and testing), as given in Equation (4).

$$\begin{array}{l}s=(((q)+(\varphi ))\times (\exp (((((\varphi )\times -(9.2282))-(\varphi ))\times \exp (-1.2522)){)}^2))+(\exp (((\exp ((q))+((N)\times (B)))+(((\mathrm{d})-1.9502)\\ -(2.733+(N)))))-(c))+\\ (\exp (((\exp ((q))+((\mathrm{d}{)}^2))-(((c)+6.543)+((\varphi )-(\mathrm{d})))))+(c)+\\ ((K)\times (((((K)-8.439)+((L/B)+(\varphi )))+(4.5804-(\mathrm{d})))\times \exp ((4.4445-8.4394))))+\\ (((((\mathrm{d})\times (q))/((z/B)+2.982))\times (((N)+(L/B))\times (z/B)))\times \\ (((\mathrm{d})\times (q))+((L/B)\times (L/B)))))\end{array}$$

(4)

5.1. Computational Modelling

Several performance metrics were calculated to evaluate the accuracy of the developed GEP models. These include the coefficient of determination (R²), adjusted (Adj.) R², mean absolute error (MAE), root-mean-square error (RMSE), Nash–Sutcliffe efficiency (NSC), Willmott’s index (WI), variance accounted for (VAF), and the performance index (PI). These metrics were used to assess the predictive power and reliability of each model. The equations for these metrics are detailed in

Table 2. Statistical indices used in this study and their mathematical form.

Indices	Mathematical Equations
R2	$\frac{{\sum _{i=1}^n(s_{{rsf}_m}-{\overline{s}}_{{rsf}_m})}^2-{\sum _{i=1}^n(s_{{rsf}_m}-s_{{rsf}_p})}^2}{{\sum _{i=1}^n(s_{{rsf}_m}-{\overline{s}}_{{rsf}_m})}^2}$
Adj.R2	$1-(1-R^2)\frac{n-1}{n-p-1}$
MAE	$\frac{1}{n}{\displaystyle \sum _{i=1}^n}\left|s_{{rsf}_p-}{s}_{{rsf}_m}\right|$
RMSE	$\sqrt{\frac{1}{n}{\displaystyle \sum _{i=1}^n}{\left((s_{{rsf}_m}-s_{{rsf}_p})\right)}^2}$
NSC	$1-\frac{{\sum _{i=1}^n(s_{{rsf}_m}-s_{{rsf}_p})}^2}{{\sum _{i=1}^n(s_{{rsf}_m}-{\overline{s}}_{{rsf}_p})}^2}$
WI	$1-\frac{{\sum _{i=1}^n(s_{{rsf}_m}-s_{{rsf}_p})}^2}{\sum _{i=1}^n{\left\{\left|s_{{rsf}_p-}{\overline{s}}_{{rsf}_p}\right|+\left|s_{{rsf}_m-}{\overline{s}}_{{rsf}_m}\right|\right\}}^2}$
VAF	$100\left(1-\frac{var\left(s_{{rsf}_m}-s_{{rsf}_p}\right)}{var\left(s_{{rsf}_m}\right)}\right)$
PI	$0.01VAF+adj.R^2-RMSE$
where n is the number of data points, p is the number of predictor variables, $s_{{rsf}_m}$ is the measured settlement of the RSF of the ith sample, $s_{{rsf}_p}$ is the predicted settlement of the RSF of the ith sample, ${\overline{s}}_{{rsf}_m}$ is the mean of the measured settlement values of the RSF, and ${\overline{s}}_{{rsf}_p}$ is the mean of the predicted settlement values of the RSF.

[57,58,59].

The ideal values of these parameters, namely, R², Adj.R², MAE, RMSE, NSC, WI, VAF, and PI are 1, 1, 0, 0, 1, 1, 100%, and 2, respectively [49,50,51]. These metrics effectively illustrate the accuracy, error, and variance in the comparison between actual and predicted RSF settlement values. The performance of the developed GEP models on both training and testing datasets is detailed in

Table 3. Performance indices for the training dataset.

Indices	GEP-M1	GEP-M2	GEP-M3	GEP-M4	GEP-M5	GEP-M6	GEP-M7	GEP-M8	GEP-M9	GEP-M10
R2	0.431	0.872	0.891	0.885	0.909	0.926	0.897	0.876	0.938	0.928
Adj.R2	0.417	0.869	0.889	0.882	0.907	0.924	0.894	0.873	0.936	0.926
MAE	0.089	0.035	0.042	0.03	0.031	0.027	0.038	0.029	0.026	0.038
RMSE	0.167	0.057	0.057	0.053	0.047	0.042	0.055	0.055	0.04	0.053
NSC	-0.152	0.867	0.865	0.883	0.909	0.926	0.877	0.875	0.936	0.884
WI	0.727	0.964	0.968	0.969	0.976	0.981	0.969	0.965	0.984	0.967
VAF	16.67	87.155	86.808	88.474	90.928	92.605	89.174	87.551	93.69	92.058
PI	0.417	1.683	1.699	1.714	1.77	1.808	1.731	1.693	1.833	1.794
Indices	Scoring and Ranking
R2	1	2	5	4	7	8	6	3	10	9
Adj.R2	1	2	5	4	7	8	6	3	10	9
MAE	1	2	5	7	6	9	4	8	10	3
RMSE	1	3	2	7	8	9	5	4	10	6
NSC	1	3	2	6	8	9	5	4	10	7
WI	1	2	5	7	8	9	6	3	10	4
VAF	1	3	2	5	7	9	6	4	10	8
PI	1	2	4	5	7	9	6	3	10	8
Total score	8	19	30	45	58	70	44	32	80	54
Rank	10	9	8	5	3	2	6	7	1	4

Note: Bold values indicate best obtained values.

and

Table 4. Performance indices for the testing dataset.

Indices	GEP-M1	GEP-M2	GEP-M3	GEP-M4	GEP-M5	GEP-M6	GEP-M7	GEP-M8	GEP-M9	GEP-M10
R2	0.727	0.887	0.929	0.906	0.879	0.909	0.95	0.902	0.961	0.946
Adj.R2	0.698	0.875	0.921	0.896	0.866	0.9	0.945	0.892	0.956	0.94
MAE	0.103	0.05	0.046	0.042	0.046	0.035	0.042	0.043	0.034	0.044
RMSE	0.17	0.098	0.072	0.092	0.1	0.082	0.057	0.097	0.049	0.071
NSC	0.529	0.844	0.915	0.863	0.835	0.89	0.946	0.845	0.961	0.917
WI	0.833	0.853	0.978	0.955	0.944	0.966	0.986	0.948	0.99	0.976
VAF	70.051	86.178	91.831	87.392	84.16	89.291	94.917	85.815	96.06	93.451
PI	1.229	1.639	1.767	1.678	1.608	1.711	1.836	1.653	1.868	1.804
Indices	Scoring and Ranking
R2	1	3	7	5	2	6	9	4	10	8
Adj.R2	1	3	7	5	2	6	9	4	10	8
MAE	1	2	4	8	3	9	7	6	10	5
RMSE	1	3	7	5	2	6	9	4	10	8
NSC	1	3	7	5	2	6	9	4	10	8
WI	1	2	8	5	3	6	9	4	10	7
VAF	1	4	7	5	2	6	9	3	10	8
PI	1	3	7	5	2	6	9	4	10	8
Total score	8	23	54	43	18	51	70	33	80	60
Rank	10	8	4	6	9	5	2	7	1	3

Note: Bold values indicate best obtained values.

. Along with the indices, the total scores and ranking were also provided. The model with the highest R², Adj.R², NSC, WI, VAF, and PI and lowest MAE and RMSE achieves the best ranking. It is noteworthy that each index is scored out of 10 and the total score is obtained by the summation of the partial scores. It can be observed that the GEP-M9 has achieved the highest accuracy (rank = 1) with an R² of 0.938 and an RMSE of 0.040 (total score = 80). Moreover, the GEP-based models (M5, M6, and M10) have attained more than 90% accuracy in estimating the settlement of RSFs with the total scores of 58, 70, and 54, respectively. Apart from the model (GEP-M1), which has obtained below 50% accuracy (total score = 8), other models have achieved 80% accuracy (R² = 0.800). Moreover, the VAF and PI values of 93.6% and 1.833, respectively, also substantiate that the GEP-M9 model is associated with the least bias and variance. In the testing phase, the GEP-M9 model continued its superior performance with a total score of 80, and standout R² and RMSE values of 0.961 and 0.049, respectively. In contrast, the GEP-M1 model displayed the weakest predictive ability for the RSF settlement, with R² and RMSE values of 0.727 and 0.170, respectively. Therefore, the GEP-M9 model is identified as having the highest predictive accuracy among the models for estimating the RSF settlement.

5.2. Sensitivity Analysis and Independent Studies Validation

Scholars and researchers have highlighted that evaluating the predictive capacity of an ML model should not rely solely on statistical metrics [37,58,60]. A robust model is one that not only predicts the output with reasonable accuracy, but also simulates it in a manner consistent with the overall underlying physical behavior of the system under investigation. Another vital measure of a ML model’s precision is its capability to accurately forecast data that are completely independent and not included in either the training or testing datasets. Therefore, the constructed GEP model is subjected to a one-at-a-time parametric sensitivity analysis. In this process, each input parameter is individually adjusted from its lowest to highest value within the data range that was used, while keeping all other parameters fixed at their average values [60,61]. The variation of input parameters is performed in terms of equal steps (20 steps), and the corresponding output (settlement) is recorded. Furthermore, after computing the settlement, a sensitivity index (SI) is calculated at each step according to the procedure described by Khosrojerdi [41]. Subsequently, the features are ranked based on the average sensitivity index according to their importance.

The strength and robustness of the GEP-based design formula is assessed by comparing the anticipated RSF settlement values with the actual physical behavior of RSF settlement predictions derived from well-established geotechnical data. Figure 6 presents the outcomes of the sensitivity analysis. The comparison reveals that the direction of the trends predicted using the synthetic data aligns with the expected trends informed by the physical principles underlying the settlement prediction.

For example, an increase in the q, z/B, B, L/B, and d logically and expectedly leads to higher settlement values. Similarly, extensive research on geosynthetic RSFs has demonstrated that enhancements in shear strength parameters (c and ϕ), K, and the N contribute to reduced settlement values under identical applied pressures (e.g., [34,37,41]). Consequently, it can be inferred that the developed GEP model accurately predicts settlement values in a rational manner, with the predicted trends aligning logically across the data range. Furthermore, the outcomes of the feature importance analysis, as depicted in Figure 7, highlight ϕ, q, d, and N as the most critical factors in determining the settlement of the RSF.

The developed model has predicted the training and testing data with reasonable accuracy and also demonstrated its ability to generalize well via a parametric sensitivity analysis. However, the developed GEP is subjected to further rigorous testing by evaluating its predictive veracity against entirely independent real-life experimental data. For this purpose, the result of the experiments of Adams and Collin [7] and Gabr and Hart [62] were compared with the predicted results (see Figure 8). It can be observed that the values predicted by the GEP model are in reasonable agreement with the experimental values of the independent studies. Hence, to this point, it can be concurred that the developed GEP model predicts the settlement of the RSF in an intelligent and reasonable way.

5.3. Probabilistic Assessment

In this section, the result of the probabilistic assessment in the form of the reliability index (β) and probability of failure (P_f) are presented. As stated earlier, the best GEP model was chosen to perform the probabilistic assessment of RSFs under service loading conditions. Multiple sets of input parameters with varying coefficients of variation (COVs) of non-deterministic input parameters were developed for this purpose. Assuming the constant mean for the non-deterministic soil variables (i.e., ϕ and c), the value of the standard deviation (SD) was calculated as follows: SD = mean × COV.

The review of the previous research shows that the range of COV for ϕ varies in the range of 3–13%, whereas, for c, the range is 10–70% [63,64,65]. For this study, five different combinations of COVs were considered, namely, C-1, C-2, C-3, C-4, and C-5, taken as 3%, 5%, 7%, 10%, and 13%, respectively, for ϕ, and as 5%, 10%, 15%, 20%, and 30%, respectively, for c. For each combination, various arrangements of deterministic parameters (i.e., K, z/B, B, L/B, N, and d) were developed by considering the mean value of ϕ as 40° and 35° with a c value of 1. The arrangements are summarized in

Table 5. Arrangements (configurations) of the deterministic parameters.

Arrangements	Deterministic Parameter Values
A-1	K = 1000; z/B = 0.3; B = 1; L/B = 2; N =3; d = 1.2
A-2	K = 1000; z/B = 0.3; B = 2; L/B = 4; N =3; d = 1.5
A-3	K = 1000; z/B = 0.3; B = 3; L/B = 6; N =3; d = 1.8
A-4	K = 1500; z/B = 0.2; B = 2; L/B = 2; N =1; d = 1.5
A-5	K = 1500; z/B = 0.2; B = 2; L/B = 2; N =2; d = 1.5
A-6	K = 500; z/B = 0.2; B = 3; L/B = 3; N =3; d = 2

. It may be noted that the configurations given in Table 5 are typical of RSF designs.

For each arrangement, the reliability index was calculated by considering μ_c values of 25, 30, 35, 40, 45, and 50. Accordingly, five reliability indices β25, β30, β35, β40, β45, and β50, along with the respective P_f values, were estimated for the service loads of 100, 200, 300, 400, 500, and 600 kPa. The steps involved in the FOSM analysis in this study can be summarized as follows: (i) the establishment of the mean values and COVs of non-deterministic parameters (ϕ and c); (ii) an estimate of the SD as the product of the mean and covariance; (iii) the generation of random values of the parameters according to the mean values and SD (100 instances for each combination in this study); (iv) the selection of deterministic input parameters (i.e., K, N, z/B, B, L/B, and d) as given in Table 5; (v) the generation of simulated datasets for the configurations given in Table 5 and its normalization with respect to the training data range; (vi) the prediction of the settlement of the geosynthetic RSFs of the simulated datasets using GEP-M9 under various service loads (100–600 kPa); and (vii) the estimation of the reliability indices (β25–β50) and probabilities of failure (P_f-25–P_f-50) for each loading step (100–600 kPa). As an illustration, an example of RA for A–I is represented in Figure 9. It can be noted that, with the increase in COV, the reliability indices decrease.

Terzaghi et al. [66] recommended a maximum allowable settlement of 25 mm for isolated shallow foundations, while Skempton and MacDonald [67] proposed a limit of 40 mm. Given a permissible settlement of 25 mm, the calculated reliability indices (β25) for the specified configuration were 3, 1.9, and 1.2 for service loads of 400 kPa, 500 kPa, and 600 kPa, respectively. For a permissible settlement of 40 mm, the reliability indices (β40) were found to be 5.5, 3.8, and 2.7 for the same service loads. The targeted range for β in geotechnical engineering systems is typically between 2 and 4, according to the literature [68].

Considering the upper-bound limit of the targeted reliability index, the maximum severity, that is, a decrease in the reliability can be noted for the higher loading conditions (400–600 kPa) with the increase in COV. The analysis of P_f was conducted at both levels, and the results are illustrated in Figure 10. If the maximum settlement was considered to be 25 mm, then the P_f showed a significant increase from 2.87 to 11.12 when C-4 and C-5 are compared, which is almost a four-fold increase in the probability of failure at 600 kPa. Similarly, for the same comparison of COV cases at 500 kPa, the P_f also increased from 0.20 to 2.55 (about 12 times). These results showed that the increase in the COV had a detrimental effect on the reliability index, which, in turn, increased the P_f. These trends are especially true for the higher loading conditions considered in the study (400 kPa to 600 kPa).

The results of the RA for A–I showed that the reliability indices decreased and the probability of failure increased with the increase in COV. For all the common arrangements presented in Table 5, the minimum values of reliability indices (obtained for C5) at all the loading levels were calculated and presented in Figure 11 and Figure 12, respectively. It can be noted that the mean values for the friction angle were 40° and 35° and the cohesion value was 1.

For the same pressure, if the footing size, B, is increased, the pressure bulb that has formed in the soil beneath the footing exerts compression over a considerably larger volume and depth of soil, resulting in more settling of the footing [69]. Hence, apart from the configurations summarized in Table 5, the effect of the width of the footing is also studied for all the COV cases under the service loading conditions. The results of the C-1 and C-5 for A–I are presented in Figure 13 and Figure 14, respectively. It is noteworthy that, in this case, the B values were raised from 1 m to 3 m and the maximum permissible settlement was considered to be 25 mm. Under service loads (300–600 kPa), the value of P_f increases with the increase in the B (Figure 13), but the values can be considered safe, with β significantly higher than the upper-bound limit (i.e., 4); hence, the RSF design may be considered reliable. However, the same is not true for the C-5 case as represented in Figure 14.

The most pronounced effects are observed under higher loading conditions (500 and 600 kPa) for case C-5 (Figure 14), where an increase in B results in P_f rising from 3.2% to 4.1% and from 12.7% to 15.1%, respectively, for 500 kPa and 600 kPa. These findings underscore the significant impact of a higher COV on the failure probability of RSFs, particularly under higher loading conditions. It can also be observed that, with the increase in the B value, β decreases and P_f increases at a higher COV level. This demonstrates that the B also influences the failure probability of RSFs. Moreover, the choice of COV values significantly affects the likelihood of RSF failure; therefore, it is necessary to identify the precise COV value throughout the reliability analysis process. Moreover, a detailed numerical example with calculations of the reliability index and the probability of failure is presented in the Supplementary File (see Tables S3 and S4), attached as a supplementary file. The input values of Gabr and Hart [62] were used to demonstrate the calculation procedure.

Another important advantage of this hybrid probabilistic modelling can be related to the amount of time that could be saved from expensive computational numerical modelling. For example, a single simulation run (loading stages 50–600 kPa) in PLAXIS 3D has been conducted using Intel Core (TM) i7-7700 CPU (2 × 3.60GHz) with 32 GB of installed memory, and each simulation approximately took 30–50 min on average. Considering only the arrangements and COV combinations presented in the paper and the lower bound of the run time (30 min), the simulation may take a minimum of 37.5 days to estimate the settlement values of RSFs. Contrarily, the GEP modelling along with the RA only took a few minutes to compute the settlement with reasonable accuracy.

6. Summary and Conclusions

For practitioners, the estimation of a geosynthetic RSF under service loads is of the utmost importance. The evaluation of the settlement of the RSF via numerical modeling is computationally expensive. In this study, the GEP modeling technique is utilized for the probabilistic analysis of RSFs. The database generated via validated large-scale PLAXIS-3D modeling in the past was used to train and test the ten GEP model frameworks. Based on the results obtained in the testing dataset (R² = 0.961 and RMSE = 0.049), GEP-M9 has shown the best performance among the ten developed GEP models. The predictive strength and rationality of the model (GEP-M9) were corroborated via a parametric sensitivity analysis and independent validation studies. Moreover, the suggested method provides an alternative to the numerical method and makes it possible to apply the GEP-based RA as a seamless addition to the deterministic analysis. The created model can be used by geotechnical engineers to predict the P_f of geosynthetic RSFs. For this study, the P_f is studied at various reliability levels (β25, β30, β35, β40, β45, and β50) in relation to different COV concentrations. Based on the results presented in this paper, the following general conclusions can be drawn:

• The best-performing GEP model (GEP-M9) for predicting the settlement of the RSF is obtained with 150, 10, and 5 chromosomes, head size, and genes. The proposed model predicts the settlement of the RSF with high accuracy for training (R² = 0.938, Adj.R² = 0.936, MAE = 0.026, RMSE = 0.040, NSC = 0.936, WI = 0.984, VAF = 93.69%, and PI = 1.833) and testing (R² = 0.961, Adj.R² = 0.956, MAE = 0.044, RMSE = 0.071, NSC = 0.917, WI = 0.990, VAF = 96.45%, and PI = 1.804) datasets. Moreover, the GEP-M9 has obtained the highest total score of 80, both in the training and testing dataset. The GEP-M1 showed the worst performance with a total score of 8, both in the training and testing dataset.
• The rationality and robustness of the GEP-M9 model were substantiated by the underlying physical behavior of the settlement of RSFs, in the parlance of geotechnical engineering. A feature importance analysis reveals that ϕ, q, d, and N are the most important features in estimating the settlement of RSFs.
• The suggested GEP-based RA technique provides rational and reasonable solutions to predict the P_f, taking into account the uncertainty of soil parameters. The results show that β and P_f are significantly affected by the COV of soil properties (ϕ and c). The increase in the COV results in a decrease in β and increase in P_f. Hence, it can be deduced that, by taking into account the uncertainty related to ϕ and c, the suggested probabilistic technique allows RSF settlement to be described as a cumulative probability distribution function, which estimates the settlement of the RSF corresponding to particular reliability levels.
• The choice of COV values significantly affects the likelihood of RSF failure. Moreover, RSF failure is affected by the width of the footing. The effect is more pronounced at higher loading levels (400–600 kPa) when compared to lower loading levels (100–300 kPa). Therefore, it is essential to classify the accurate COV value throughout the course of RA.

Although the developed GEP model predicts the settlement of a geosynthetic RSF with reasonable accuracy, however, the AI models are limited to the range of the data it is trained on and the GEP is no exemption to this limitation. Hence, care should be applied when predicting the settlement beyond the range of the data. Moreover, the model is only applied to estimating the immediate settlement. The creep effect is not taken into account. It may be noted that the effects of c and ϕ were not computed independently; rather, the combined influence of soil characteristics was explored across a range of COVs. The uncertainty associated with the GEP model is not considered in this study, and the only source of uncertainty is the soil parameters.