The Data-Driven Homogenization of Mohr–Coulomb Parameters Based on a Bayesian Optimized Back Propagation Artificial Neural Network (BP-ANN)

Abstract

Homogenization methods can characterize the mechanical properties of these materials based on appropriate constitutive models and data. They are also applied to the characterization of mechanical parameters under complex geotechnical conditions in geotechnical engineering because of the complexity and heterogeneous nature of geotechnical materials. Unfortunately, existing homogenization methods for geotechnical mechanical parameters often incur immense computational costs. Hence, a framework that utilizes finite element analysis for generating a dataset which is then trained using a Bayesian Optimized Back Propagation Artificial Neural Network (BP-ANN) to obtain the homogenized Mohr–Coulomb parameters of the soils is proposed. This is the first time that Bayesian optimization and a BP-ANN have been used in conjunction to predict the homogenized mechanical parameters of soils. The dataset used for training the data is generated using the commercial FEM software ABAQUS (6.10). The maximum difference between the top and bottom part of the tunnel of the heterogeneous model and homogeneous model of our test cases only varies by 5.3%, thereby verifying the excellence of the Bayesian Optimized BP-ANN.

1. Introduction

For centuries, scientists have discussed how to effectively carry out geotechnical engineering practices in the face of complex, non-uniform, and discontinuous geotechnical conditions. The discontinuity, inhomogeneity, and anisotropy of geotechnical materials pose great difficulties for engineering practices, such as collapses or landslides that may occur during underground excavation. In order to ensure the smooth progress of engineering practices, scientists and engineers have researched, discussed, and used many approaches to try to model and simulate complex geotechnical conditions. Many approaches in this area fall within the following categories: continuous or discontinuous numerical methods, such as the finite difference method (FDM), finite volume method (FVM), finite element method (FEM), discrete element method (DEM), and discrete fracture network (DFN), etc. In addition, the homogenization methods used to simplify complex heterogeneous materials into homogeneous materials with equivalent properties are also applied in this field.

Homogenization methods offer a valuable means to simplify and characterize the mechanical properties of complex heterogeneous materials. By utilizing appropriate constitutive models and relevant data, these methods enable the representation of complex materials as equivalent, homogeneous materials with simplified properties. This simplification allows for the computationally efficient calculations of material responses. So far, scientists have discussed how to achieve the characterization of these complex materials and widely apply them to the research and development of composite materials. In the field of geotechnical engineering, complex and discontinuous geotechnical conditions add many difficulties to engineering practices. Due to the performance of the homogenization method when facing complex heterogeneous materials, it is also applied in the characterization of geotechnical mechanical parameters. For example, Chang developed a meso-macro constitutive model for saturated frozen saline sand based on homogenization theory, predicting the deformation of frozen saline sand under different salt contents and confining pressures, and reproducing strain-hardening/softening, high-expansion, and pressure-melting characteristics [1]. Gasmi proposed a new numerical homogenization method called HLA-Dissim which simulates the discontinuity network of real rock masses based on the International Society for Rock Mechanics (ISRM) scan line field survey method [2]. Shi studied the impact of the dilatancy effect of layered rock with rough interfaces on its mechanical behavior, used a homogenization method to simplify the layered rock structure into an equivalent continuous medium, defined interfaces with different roughness to evaluate the dilatancy effect, and used a dilatancy model called the contact density model to estimate the stiffness characteristics of rough surfaces [3]. These studies suggest that homogenization methods exhibit excellent performance in the field of geomechanics and engineering. Overall, the utilization of homogenization methods in geomechanics and engineering demonstrates their exceptional performance in tackling the difficulties posed by complex geotechnical conditions. These methods present a powerful tool for researchers to gain insights, make accurate predictions, and effectively address challenges encountered in geotechnical engineering practices. However, the field of the homogenization of nonlinear materials, such as soil, is relatively underexplored compared to other areas of research in geotechnical engineering. While there have been significant advancements in the characterization of complex heterogeneous materials using homogenization methods, the application of these methods specifically to nonlinear materials like soil has received limited attention.

At the same time, with the development of computer science and technology, new methods of machine learning and neural networks are widely used in the field of computer vision (image recognition, object detection, etc.), natural language processing (speech recognition, sentiment analysis), and industries such as finance and transportation. A notable feature of these new methods is that they can learn the features of the given datasets well and make precise predictions, to a certain extent, and a more detailed introduction to neural networks, and the development of the whole methodology may be found in [4].

Scientists have keenly discovered the good prospects of these new technologies in the field of materials and have discussed them extensively. Logarzo used machine learning to homogenize the inelastic behavior of materials and trained intelligent constitutive laws using recurrent neural networks [5]. The strains generated by the materials created by the Terzaghi principle in a porous medium in the context of tunnel works have been resolved with the use of FEM implemented with Matlab [6]. Ismail used an ANN to identify the properties of uncertain materials in unidirectional layers and used plug-ins in ABAQUS to generate virtual data to verify that the ANN could clearly and accurately parameterize the relationship between fiber mechanical properties on a microscopic scale and fiber/matrix interface parameters and mechanical properties on a macroscopic scale [7]. Peng discusses the latest technologies, applications and opportunities of combining machine learning and multiscale modeling and raises some open questions and potential challenges and limitations in their article [8]. Eidgahee compared several data-driven, machine learning tools to predict the dynamic modulus of asphalt mixtures with high accuracy [9]. Frankel employed data-driven homogenization approaches to predict the mechanical response of oligocrystals composed of multiple crystals [10]. The application of neural networks and machine learning techniques in the field of material science can effectively handle the relationship between material properties and relevant parameters [11]. Rehman predicted the hydraulic conductivity of sandy soils using machine learning methods and gave a detailed introduction to AI optimization methods utilized in geotechnical engineering [12,13,14]. The existing research in the field of homogenization has two notable drawbacks. Firstly, both analytical and numerical homogenization methods, including the previously proposed asymptotic homogenization method, tend to primarily focus on solving elastic homogenization problems. In other words, there is a lack of sufficient research addressing inelastic homogenization problems. This lack of attention hampers the ability of engineers and researchers to effectively predict and analyze the behavior of complex geotechnical systems involving inelastic materials. Traditional analytical and numerical approaches require substantial computational resources and time-consuming iterations to achieve accurate results. This can be challenging when dealing with complex geotechnical problems that require efficient and practical solutions. However, one potential solution to address this issue is the utilization of neural network models. Neural networks offer a promising alternative due to their ability to provide rapid predictions through a single forward pass, resulting in negligible time costs. Once trained, a neural network can quickly generate homogenized results for a given set of inputs. This computational efficiency makes neural networks an attractive option for homogenization, potentially reducing the time and resources required for analysis.

In this work, a Bayesian Optimized BP-ANN for the homogenization of Mohr–Coulomb mechanical parameters, based on its excellent performance in predicting output values, is proposed [15]. The structure of a back propagation artificial neural network model (BP-ANN) is defined and then trained with a dataset generated using the Python scripting interface of the commercial FEM software ABAQUS. In the process of constructing the dataset, the mechanical parameters of heterogeneous soil layers are selected based on the distribution characteristics of the mechanical parameters of heterogeneous soil layers in actual situations to enhance the reliability and generalization performance of the dataset. At the same time, Bayesian optimization is also used to improve the performance of the BP-ANN model by optimizing the hyperparameters to minimize the loss function. By combining the power of the BP-ANN and Bayesian optimization, we aim to achieve fast, accurate, and reliable homogenization of the Mohr–Coulomb mechanical parameters of soil. Also, this work may promote the utilization of FEM as a means of generating virtual training data in large quantities.

2. Methodology

2.1. Back Propagation Artificial Neural Networks

Neural networks, especially the BP-ANN, have shown great potential in predicting and mapping input–output relationships in recent years. Due to their powerful prediction capabilities, they have been applied in structural damage identification, structural health monitoring, deep foundation pit deformation prediction, and concrete-related performance prediction [16,17]. Here, Multiple Layer Perceptron (MLP) will be used as the overall structure of the BP-ANN and the dropout method will be applied to further mitigate the possibility of overfitting [18].

Since BP-ANNs are widely used to find the mapping between the input and output data, which is sometimes not achievable using a simple regression algorithm because of the potential non-linearity between the input and output, activation functions can introduce non-linearity to the neural network by introducing an additional step at each layer during the forward propagation, resulting in a slightly improved computational cost [19].

Under the conditions of specifying the loss function and neural network architecture, it is important to select an appropriate optimizer to train the network in order to minimize the loss function. The ADAM (adaptive moment estimation) optimizer is a widely used optimizer, which is a fast and computationally efficient optimizer, and is very suitable for problems with many parameters [20].

2.2. Bayesian Optimization

Bayesian optimization originates from Bayesian theorem, which is a method of using Bayesian formulas to establish the probability distribution of the optimization process.

$$p\left(f\right|D_{1:t})=\frac{p\left(D_{1:t}\right|f\left)p\right(f)}{p\left(D_{1:t}\right)}$$

(1)

$$D_{1:t}=\{\left(x_0,y_0\right),\cdots ,\left(x_t,y_t\right)\}$$

(2)

Among them, $f\left(·\right)$ represents the objective function for evaluating performance. $x_t$ is the decision vector. $y_t=f\left(x_t\right)+{\epsilon }_t$ is the evaluation value of the objective function, ${\epsilon }_t$ is the observation noise. $D_{1:t}$ is the set of observed values. $p\left(D_{1:t}\right|f)$ is the likelihood distribution of the observed value $y$. $p\left(f\right)$ is the prior probability distribution of $f$, which assumes the state of the unknown objective function. $p\left(D_{1:t}\right)$ represents the marginal likelihood distribution. $p\left(f\right|D_{1:t})$ is the posterior probability distribution of $f$, which describes the confidence level of the unknown objective function after modifying the prior probability distribution through the observed dataset.

Bayesian optimization fully utilizes the information of the previous sampling point and finds the hyperparameter combination that globally optimizes the results by learning the shape of the objective function. The Bayesian optimization framework mainly consists of two core parts: the probability proxy model and the collection function. This article uses Gaussian processes (GP) as the proxy model and Expected Improvement (EI) as the collection function. The algorithm process is shown in Figure 1. Firstly, a probability proxy model is used to fit the objective function. then, the collection function is used to determine the next optimal sampling point based on the posterior probability of past observation points from areas that have not been sampled or may have optimal solutions. Finally, the corresponding parameter set is updated until the stop condition is triggered. The parameter optimization process is represented as follows:

$$\begin{gathered} {\theta }_{t+1}=argmaxg\left({\theta }_t\right) \\ \theta \in \mathsf{\Psi } \end{gathered}$$

(3)

In the equation, $g\left(·\right)$ is the collection function; $\theta$ is the optimization parameter combination at time t; and $\mathsf{\Psi }$ is a set of optimized parameter spaces.

2.3. Back Propagation Artificial Neural Network Based on Bayesian Optimization

The algorithm flow of the back propagation artificial neural network based on Bayesian optimization is shown in Figure 2.

(1)

Divide the dataset into training and testing sets.

(2)

Set the hyperparameters that need to be optimized and the range of values for each hyperparameter. In this work, we selected hyperparameters such as learning rate, dropout rate, the number of hidden layers, and the number of neurons in each hidden layer.

(3)

Set the initial values of hyperparameters, perform Bayesian optimization, and obtain the optimal hyperparameters.

(4)

Substitute the optimal hyperparameters into the back propagation artificial neural network to determine the final architecture of the back propagation artificial neural network.

(5)

Train models on the training set.

(6)

Test the model on the test set.

3. Generation of the Training Dataset

The generation of our training dataset can be broken down into three stages, which is further illustrated in Figure 3. In the first one, the finite element method is employed using the commercial software ABAQUS on 1296 heterogeneous models, each with a different set of soil parameters, and the maximum difference in displacement is then exported. In the second stage, FEM is employed again, but this time on homogeneous models with various soil parameters, and the maximum difference in displacement is exported. In the third stage, the acquired maximum difference in displacement of the heterogeneous models and homogeneous models is compared. The two groups of parameters would be considered equivalent only if the difference between the two values is below a certain threshold, which would be $5\times {10}^{-5}$ in our work.

3.1. Problem Formulation

Before we discuss the details of the numerical simulations, we must first specify the parameters of the physical system. Due to the length of the tunnel, we have simplified the problem into a planar strain problem.

A schematic diagram of our heterogeneous geometric model is shown in Figure 4. Our geometric model, from top to bottom, consists of 4 laminae of soil with a thickness of 2 m. The tunnel, encompassed by the four laminae of soil, is of a hollow rectangular shape, with a width of 10 m and a height of 4 m. The thickness of the tunnel is 50 cm. Below the tunnel, a thick layer of rock is added to simulate the propagation process of the seismic wave. The approximate size of an element is 0.25 m, which is then validated through mesh sensitivity analysis.

Although there are numerous constitutive laws available to describe the behavior of soil under dynamic loads [21,22,23], the Mohr–Coulomb model is selected for the sake of simplicity. The constitutive law of the four different types of soils may be specified using the Mohr–Coulomb model, a constitutive model widely used in modeling the stress–strain behavior of rocks [24,25]. These soil layers, from top to bottom, consist of silty clay, fine sand, pebble, and silty clay.

The conventional Mohr–Coulomb model uses four parameters to define the stress–strain behavior of the material: elastic modulus E, Poisson’s ratio υ, cohesion force c, and friction angle ϕ. Additionally, we characterize the behavior of the tunnel material as C20 concrete using linear elasticity.

To prevent seismic waves from reflecting, an infinite boundary condition is imposed on the left and right side of the planar model [26,27]. CINPS4 elements have been used for the infinite boundaries, while the remainder of the model utilized CPE4R elements. The upper surface of the model was assumed to have a free boundary.

The loading condition for this model was implemented using seismic amplitude data from the Tangshan earthquake, whose acceleration magnitude is illustrated in Figure 5. The frequency domain analysis of the input wave is depicted in Figure 6. The loading condition is then defined by adding acceleration at the bottom nodes of the model.

3.2. Finite Element Simulation of the Heterogeneous Model

By utilizing the Python scripting interface of ABAQUS, the dynamic physical problems are solved numerically, and the training dataset is finally acquired. Specifically, we focused on the difference in horizontal displacement between the top and bottom of the tunnel. This difference value was recorded as a function of time, and the maximum value among all time points is selected. All the maximum displacement difference–time curves are collected, and one of these curves is illustrated in the Figure 7 and Figure 8.

The horizontal displacement contour on the deformed configuration of the selected set of parameters may also be obtained.

The maximum difference in horizontal displacement is measured at 1296 different sets of material parameters based on the randomness and distribution characteristics of them.

As mentioned above, mesh sensitivity analysis is then employed to ensure the size of our elements is small enough to obtain an accurate solution. The computed displacement is now plotted against the element size in Figure 9.

Simulations with several different element sizes are conducted and the maximum displacement difference values are obtained. The displacement difference converges at around 0.014 m, and a mesh size of 0.25 is chosen since it may accurately represent the features of the physical problem without incurring a huge computational overhead.

3.3. Finite Element Simulation of the Homogenized Model

In this section, grid search is employed on the homogenized parameters to find the equivalent values of the heterogeneous soil parameters. We start from an average of the heterogeneous parameters. Again, the Python scripting interface of ABAQUS is invoked to generate an adequate maximum difference in horizontal displacement at any given time.

3.4. Matching the Homogenized Parameters to the Heterogeneous Parameters

In this research, the homogenization of the four laminae of soils by matching the displacement values is achieved. Two sets of parameters are considered equivalent if their difference is lower than a certain threshold, and, in this case, lower than 0.00005 m. Ultimately, only 247 entries met this criterion and were added to the final output database. Then, the dataset is further divided into the training dataset and the testing dataset. The statistics of the input and output parameters of our dataset are shown in

Table 1. Statistics of the input and output parameters.

Parameter	Min	Max	Avg	STD
$E_1$	$9.60\times {10}^7$	$8.64\times {10}^8$	$4.99\times {10}^8$	$3.40\times {10}^8$
$C_1$	$2.00\times {10}^4$	$2.50\times {10}^4$	$2.31\times {10}^4$	$2.42\times {10}^3$
${\varphi }_1$	15.00	20.00	17.63	2.50
$E_2$	$2.40\times {10}^8$	$2.16\times {10}^9$	$1.11\times {10}^9$	$8.42\times {10}^8$
$E_3$	$3.97\times {10}^8$	$3.58\times {10}^9$	$2.03\times {10}^9$	$1.40\times {10}^9$
$E_4$	$1.33\times {10}^8$	$1.20\times {10}^9$	$6.35\times {10}^8$	$4.47\times {10}^8$
$C_4$	$2.00\times {10}^4$	$2.50\times {10}^4$	$2.22\times {10}^4$	$2.48\times {10}^3$
${\varphi }_4$	15.00	20.00	17.65	2.50
$E_{pred}$	$3.50\times {10}^8$	$2.95\times {10}^9$	$1.79\times {10}^9$	$5.43\times {10}^8$
$C_{pred}$	$1.67\times {10}^3$	$1.67\times {10}^4$	$7.86\times {10}^3$	$3.34\times {10}^3$
${\varphi }_{pred}$	8.00	48.00	34.50	11.86

Units: E (pa), C (kpa), Phi (°).

.

4. Network Architecture and Parameter Optimization

4.1. Data Pre-Processing

The input parameters of the neural network have varying dimensions, so it makes sense to normalize the input parameters in order to increase the predicting precision of the neural network [28,29]. The input parameters are mapped to a range from 0 to 1 using the normalization formula below:

$$x_{norm}=\frac{x-x_{min}}{x_{max}-x_{min}}$$

(4)

where $x_{norm}$ stands for the normalized input, $x$ stands for the current original input, $x_{min}$ stands for the minimum value for all inputs, and $x_{max}$ stands for the maximum value.

4.2. Dataset Partitioning

In this study, we randomly divided the dataset into training sets and test sets, of which 80% (197 groups) were training sets and 20% (50 groups) were test sets. In the specific process of training the model, a five-fold cross-validation method was utilized to optimize each training process; that is, the training sets were divided into five parts, four of which were used for training in turn, and the remaining part was used for verification.

4.3. The Bayesian Optimized Neural Network

This study used the Pytorch library to implement the BP-ANN as well as the Bayesian optimization process. For this study, the hyperparameters of the training process of the neural network is optimized using Bayesian optimization. The results of the Bayesian optimization yields are shown in

Table 2. The results of the Bayesian optimization yields.

Neuron	Dropout_Rate	Learning_Rate	Hidden Layers
25	0.6	0.006	7

.

For the fully connected neural network (FC-NN) implemented in this paper, there are eight input variables, three output variables, and seven hidden layers. The tanh () function is selected as the activation function, the mean square error (MSE) is selected as the loss function, and the ADAM optimizer is used to optimize the model parameters. The mean square error (MSE), root mean square error (RMSE), and mean absolute error (MAE) are used to evaluate the predictive performance of the model. The relevant calculation formulas are:

$$MSE=\frac{1}{n}\sum _{i=1}^n{\left(y_{\mathit{tr}ue}-y_{pred}\right)}^2$$

(5)

$$RMSE=\sqrt{\frac{1}{n}\sum _{i=1}^n{\left(y_{\mathit{tr}ue}-y_{pred}\right)}^2}$$

(6)

$$M=\frac{1}{n}\sum _{i=1}^n\left(y_{\mathit{tr}ue}-y_{pred}\right)$$

(7)

In the formulas, $y_{\mathit{tr}ue}$ represents the variable value of the test set, and $y_{pred}$ represents the corresponding predicted value.

5. Results

5.1. Training Loss of the BP-ANN Model

After Bayesian optimization, the best hyperparameter combination was applied to the BP-ANN to obtain the loss curve under each epoch. The logarithm of the MSE loss of the training process is depicted in Figure 10 below.

5.2. Verification of Our Model on an Independent Dataset

The test set was used to evaluate the trained model by using MSE, MAE, and RMSE indicators, and the results are shown in

Table 3. Results.

	MSE	RMSE	MAE
E	0.06	0.24	0.18
$\varphi$	0.036	0.19	0.14
C	0.10	0.32	0.27

.

From 1296 sets of different Mohr–Coulomb parameters, after excluding the 247 sets of data that were successfully matched and used for training models, 70 sets of material parameters were randomly selected and the trained BP-ANN used them as inputs to obtain homogenized parameters, which were then substituted into ABAQUS to obtain the corresponding 70 sets of displacement difference. The mean absolute percentage error (MAPE) between the displacement difference obtained by these 70 sets of heterogeneous parameters and the displacement difference of the corresponding 70 sets of averaging parameters were calculated. The calculated MAPE was only 5.3%.

Thirty sets of Mohr–Coulomb parameters, which were randomly generated and not present in the original 1296 datasets, were used. These 30 sets of data were then substituted into the BP-ANN model to obtain homogenized parameters. The resulting homogenized parameters, along with the corresponding 30 sets of data, were calculated in ABAQUS. Subsequently, the 30 sets of displacement differences were obtained. Remarkably, the MAPE values of these displacement differences for all 30 sets were found to be only 2.6%. All the 100 sampled data points’ predicted displacement value and true displacement values are plotted in Figure 11. This signifies that the proposed method in this study demonstrates excellent stability and generalization capabilities.

5.3. Discussion

The methodology proposed in this study has the following advantages.

(1)

The proposed method exhibits generalizability and can be applied in various fields of geotechnical mechanics.

(2)

The established model accurately predicts the displacement difference. The mean absolute percentage error (MAPE) between the predicted and actual displacement difference is only 5.3% for the test set, and 2.6% for the independent validation set. This again demonstrates the model’s strong generalization ability for the geotechnical problems examined in this study.

(3)

Utilizing the trained model for practical geotechnical problems significantly reduces the need for computational resources and time costs when compared to traditional approaches.

6. Conclusions

We used a dataset consisting of 247 sets of heterogeneous–homogeneous Mohr–Coulomb parameters. Subsequently, we implemented the Bayesian Optimized BP-ANN using the Pytorch package. We proposed a neural network that utilized eight input parameters to evaluate the output Mohr–Coulomb parameters of the homogeneous soil.

Our work led to the following conclusions:

• Bayesian optimization, in our work, managed to optimize our network architecture as well as the hyperparameters.
• Applying the BP-ANN to predict the yields of the corresponding sets of predicted parameters and the good agreement between the difference in displacement computed using the predicted parameters and the simulation data of the heterogeneous model verify the superiority of our model.

While the Bayesian Optimized BP-ANN demonstrated high precision in predicting homogenized Mohr–Coulomb parameters, further enhancements are possible. Exploring alternative numerical or computational homogenization methods could lead to more accurate results. Additionally, constructing and utilizing a larger database would improve both the accuracy and generalization ability of the neural network model.

The fast and accurate predictions provided by the BP-ANN model for Mohr–Coulomb model parameters in homogenized soil layers make it a valuable tool in engineering applications. Also, the incorporation of simulation-generated datasets into engineering applications is a novel and valuable idea. This methodology offers the opportunity to generate large amounts of data through simulations, providing a controlled and efficient approach.