Concrete Creep Prediction Based on Improved Machine Learning and Game Theory: Modeling and Analysis Methods

Abstract

Understanding the impact of creep on the long-term mechanical features of concrete is crucial, and constructing an accurate prediction model is the key to exploring the development of concrete creep under long-term loads. Therefore, in this study, three machine learning (ML) models, a Support Vector Machine (SVM), Random Forest (RF), and Extreme Gradient Boosting Machine (XGBoost), are constructed, and the Hybrid Snake Optimization Algorithm (HSOA) is proposed, which can reduce the risk of the ML model falling into the local optimum while improving its prediction performance. Simultaneously, the contributions of the input features are ranked, and the optimal model’s prediction outcomes are explained through SHapley Additive exPlanations (SHAP). The research results show that the optimized SVM, RF, and XGBoost models increase their accuracies on the test set by 9.927%, 9.58%, and 14.1%, respectively, and the XGBoost has the highest precision in forecasting the concrete creep. The verification results of four scenarios confirm that the optimized model can precisely capture the compliance changes in long-term creep, meeting the requirements for forecasting the nature of concrete creep.

1. Introduction

Creep significantly affects concrete’s mechanical properties. It provides a reference to evaluate the mechanical features of concrete structures under long-term loading, which is the basis of the mechanical response of concrete under this condition [1,2]. The generation of creep increases the deflection of the concrete and the loss of prestressing force, causing a break in the structure and reducing the performance of the construction [3,4,5,6]. Concrete creep arises from the expansion of cracks in the transition zone at the interface between mortar and aggregate within the concrete; cracks are first generated in the transition region between coarse aggregate and mortar. The cracks are not merely related to the acceptable deformation of the cementitious material itself, but also to the uneven shrinkage caused by the internal drying of the concrete; simultaneously, the internal restraint generated by creep can also slow down the development of microcracks [7,8,9,10,11].

Factors exerting an influence on creep have been extensively investigated, such as the Young’s modulus intrinsic to the cementitious material itself, the water-to-binder ratio, the curing temperature, the relative humidity, and the age of loading of the concrete [12,13,14,15,16,17]. Recently, scholars have attempted to decouple the effects of various variables from different perspectives to make more accurate predictions of creep under complex conditions. The interaction of these variables and the high sensitivity of creep with respect to these parameters pose an even greater challenge for large-scale measurements. In an effort to further surmount the impact that hydration reactions and the interaction of input variables exert on the results of a creep experiment during a loading procedure, Wyrzykowski et al. [18] partially substituted unhydrated cement with inert quartz powder, and then conducted uniaxial compression creep experiments to calibrate the mechanical properties of the non-aging system. Ladaoui et al. [19] carried out tests on the essential creep of four kinds of high-efficiency concrete under various temperatures and found that the basic creep coefficients at the temperature increases of 20 °C and 50 °C were 2 and 3. Gan et al. [20] investigated the impact of factors such as the water–cement ratio and stress level on creep measurements by means of the micro-cantilever bending method, indicating that microscopic variations affect the overall creep. However, they are confined to decoupling a finite number of variables, thereby restricting the precise prediction of creep in actual loading.

With the extensive development of creep research, a series of experiment-based analysis models, semi-empirical models considering the coupling effects of different variables, and classical models have been applied to creep studies. Bažant and Wu [21] proposed the utilization of the Maxwell chain model to represent the rate-dependent creep law of concrete. Zhang and Hamed [22] presented the rate-dependent creep law. Nikola Tošic’ et al. [23] analyzed the difference in creep flexibility coefficients between RAC and NAC under the same conditions using the MC₂₀₁₀ model based on the establishment of a database on the creep of recycled aggregate concrete (RAC) and natural aggregate concrete (NAC). Jian Cao et al. [24,25] investigated the consequences of different fly ash contents and loading ages on the creep recovery of fly ash concrete using the modified B₄ model. Cherif Belebchouche et al. [26] investigated the creep behavior of self-compacting concrete beams under different parametric mix proportions by means of the B₃, ACI₂₀₉, EC₀₂, and CEB MC ₉₀ models. The above models provide an effective alternative to experiments for the in-depth study of the creep mechanism under the influence of multiple parameters. However, they assume ideal conditions, and there are inevitable accuracy errors in actual engineering, while MC2010, MC09, EC02, and other models only consider a limited number of input characteristics, and the B3 and B4 models need to be modified under actual engineering conditions according to a large number of experimental data on the model coefficients when considering the creep of admixed concrete, self-consolidating concrete, etc.

ML models have been extensively utilized in fields such as rubberized concrete frost resistance [27], chloride diffusion in concrete [28], the compressive strength of silica fume concrete [29], lateral confinement factors of reinforced concrete columns [30], and the mechanical features of recycled aggregate concrete [31]. ML models provide new directions for creep prediction. The Support Vector Machine (SVM) [32,33,34] can utilize kernel functions instead of the inner product of mapping to a high-dimensional space, avoiding dimensionality catastrophes and providing good robustness. Random Forest (RF) [35,36] has characteristics including controllable generalization error, a strong ability to handle high-dimensional eigenvector space, resistance to overfitting, and fewer hyperparameters for optimization, but there are still some limitations. Extreme Gradient Boost (XGBoost) [37,38] introduces the second-order partial derivatives of the loss function for weighting based on the CART regression tree, characterized by high accuracy and speed. However, XGBoost has the disadvantages of being prone to falling into local optimal solutions and high requirements for learning samples. The prediction ability of the above model is influenced by hyperparameters. Intelligent optimization algorithms can search for optimal hyperparameters to upgrade the prediction precision and robustness of machine learning (ML) models. The Snake Optimization (SO) [39] algorithm aims to search for the optimal solution by simulating the foraging behavior of snake groups. However, the standard SO algorithm has problems with the slow convergence speed in the early stage and becoming trapped in local optimal solutions. Therefore, three improved strategies are proposed based on the snake optimizer to promote the algorithm’s generalization ability in optimizing various model parameters.

In previous studies, feature importance, LIME, PDP, etc., have often been used for the interpretability analysis of models. However, these methods cannot stably explain the prediction analysis of complex models. The Shapley theory was first proposed in 1953 [40,41]. It has the characteristics of strong additivity and high stability and has been widely used in concrete strength prediction [42,43], concrete penetration resistance [44], and concrete early performance analyses [45].

In this study, the author proposes four improvement strategies for the SO algorithm to enhance its capability to leap out of local extremums. On this basis, three ML models (SVM, RF, and XGBoost) for creep prediction are developed according to the NU creep–shrinkage database. Finally, based on the SHAP method, the influence of single features and cross-features on prediction results is explored.

2. Data Processing Analysis

From 2011 to 2013,the NU database was expanded from the RILEM database by Professor Bažant’s team at Northwestern University with the support of the U.S. Department of Transportation. This expansion was based on collating more than 1400 creep test curves and 30,000 data points from many journal articles, conference records, and other documents [46,47]. Based on consolidation theory, this paper studies 4207 sets of experimental data in the NU database to construct the ML model. Thirteen parameters are chosen as model inputs: water–cement ratio (w/c), aggregate–cement ratio (a/c), cement content ((C) Kg/m³), cement type, 28-day compressive strength (fc₂₈ MPa), surface-to-volume ratio (V/S), sample height (mm), loading temperature (T °C), experimental humidity (RH_Test °C), loading stress ((Sigma) MPa), strength–stress ratio (sigma/f_ct0), time of loading or time exposed to the environment ((t′) (days)), and loading time ((t) days). Compared with other numerical parameters, in this study, the cement type is treated as a classification parameter of the input variable. In the classification of cement types, rapid-hardening cement (R), slow-hardening cement (SL), and rapid-hardening high-strength cement (RS) are, respectively, characterized as 1, 2, and 3 for classification. The measured creep compliance (Jcreep (μm/m/MPa)) in the experiment is taken as the output variable.

Table 1. The minimum, maximum, mean, and derivative values of input parameters.

Variable	Min	Max	Mean	SD	Type
w/c	0.236	0.8	0.4351	0.1145	Input
a/c	1.22	8.32	4.418	1.031	Input
C (kg/m3)	247	725	410.158	86.4	Input
fc28 (MPa)	118.9	10.8	53.07	23.13	Input
h mm)	76	1800	510.819	275.19	Input
V/S	13	129	31.58	10.97	Input
T (°C)	19	130	35.478	23.346	Input
RH_Test (%)	20	101	76.75	20.5	Input
Sigma (MPa)	0.69	46.3	16.56	9.93	Input
sigma/fct0	0.031	0.84	0.335	0.1167	Input
t (days)	1.421 × 10−14	6979	200.27	454.38	Input
Cem	1	3	/	/	Input
t′ (days)	0.66	90	30.62	18.78	Input
JCreep (μm/m/MPa)	−16.7	590.42	70.79	45.62	Output

describes the database employed for creep prediction, encompassing the minimum, maximum, mean, and derivative values (as depicted in Figure 1).

3. Methodology

The traditional grid search algorithm and random search algorithm for complex models with high dimensions and large amounts of data need to enumerate the as-summed values of each hyperparameter combination to optimize the model. Meanwhile, due to the shortcomings of the low convergence rate and local extremum of the snake optimizer, there is often a large error when enumerating hyperparameter combinations. Not only is the optimization efficiency low, but the optimization results are often unsatisfactory. In consideration of the limitations of the aforementioned optimization methods, this study puts forward the use of the HSOA optimization method to find the hyperparameters of models such as SVM, RF, and XGBoost during the cross-validation process to improve the model’s prediction efficiency. Finally, based on the optimized model, the Shapley value is calculated to characterize the final impact of different input variables on creep compliance.

3.1. Machine Learning Models

3.1.1. Support Vector Machine Models

The SVM was presented by Corinna Cortes and Vladimir Vapnik in 1995 as a generalized linear classifier [48,49]. A SVM uses interval maximization as its learning approach to transform linear problems into minimized regularized loss functions. Hence, the SVM’s learning method provides the optimum algorithm to solve convex quadratic programming. The parameters that need to be optimized under cross-validation are (1) the penalty coefficient (C), (2) the regularization parameter (penalty), (3) the kernel function parameter (sigmoid), and (4) the classification decision parameter (decision_function_shape).

3.1.2. Random Forest Model

A Random Forest (RF) is a tree-based ensemble model proposed by Breiman in 2001 [50]. The Random Forest (RF) adopts a random feature selection strategy [51]. In the sample selection stage, the Random Forest algorithm selectively and repeatedly samples k samples from the original dataset to construct weak learners, and the output result is the voting average of multiple weak learners [52,53]. The random selection of features implies that a limited number of features are utilized to split nodes during each round of data training, ensuring that all decision trees are uncorrelated with each other; the final prediction of the RF model is the mean of all predictions made by the decision trees.

The hyperparameters that need to be optimized for the Random Forest model are (1) the maximum depth of the decision tree (max_depth); (2) the minimum number of leaf nodes (min_leaf_nodes); (3) the number of decision trees (n_estimators) (the final forecasting of the RF model is the mean of all predictions made by the decision trees); and (4) the max feature subset parameter (max_features).

3.1.3. Extreme Gradient Boosted Tree Model (XGBoost)

Extreme Gradient Boosting creates a strong learner through training multiple weak learners, incorporating second-order derivatives to increase the loss function’s accuracy, regular terms to prevent overfitting, and Blick storage for parallel computing. The optimal model is built by minimizing the data’s loss function [54]. A shrinkage strategy is introduced to avoid overfitting the loss function, where the residual values fitted by each weak learner are multiplied by eta in the range [0–1]. The expression of the loss function is shown in Equation (1) [38]:

$$L=\underset{f\in F}{\min }{\displaystyle \frac{1}{N}}{\displaystyle \sum _{i=1}^{N}L(y_i},f(x_i))$$

(1)

where $X_i$ and $y_i$ describe the input and output values, respectively; $f\left(x_i\right)$ describes the ith weak learner; and N represents the total number of weak learners. To reduce the model’s complexity, Equation (1) can be rewritten after introducing $J\left(f\right)$ [38]:

$$L=\underset{f\in F}{\min }{\displaystyle \frac{1}{N}}{\displaystyle \sum _{i=1}^{N}L(y_i},f(x_i))+\lambda J(f)$$

(2)

Among them, λ is the structural risk coefficient parameter. After this is introduced, the objective function is as follows [38]:

$$obj(\theta )={\displaystyle \sum _i^{n}L(y_i},y_i)+{\displaystyle \sum _{k=1}^K\mathsf{\Omega }(f_k})$$

(3)

The objective function $Obj\left(\cdot \right)$ is divided into two terms, where the right and left terms indicate the model’s training error and regularization. The XGBoost model contains multiple CARTs, defining the complexity of each tree as follows [38]:

$$\mathsf{\Omega }(f)=\gamma T+{\displaystyle \frac{1}{2}}\lambda \parallel \omega {\parallel }^2$$

(4)

where T denotes the total number of leaf nodes, $‖\omega ‖$ denotes the leaf node vector modulus, and $\gamma$ is a model hyperparameter. A second-order Taylor expansion of the objective function can be obtained [38,55]:

$$f(x+△x)=l(y_i,{y_i}^{k-1})+g_i{f}_k(x_i)+{\displaystyle \frac{1}{2}}{h}_i{f}_k{(x_i)}^2$$

(5)

Assume that the output of the k-integrated tree models is the following [38,55]:

$${\widehat{y}}_i=\varphi ({\mathrm{x}}_i)={\displaystyle \sum _{k=1}^K{f}_k(x_i)}$$

(6)

where ${\widehat{y}}_i$ represents the prediction results of the first i decision trees for the sample. The objective function is described as follows [56]:

$$obj={\displaystyle \sum _{i=1}^{n}l(y_{i,}{{\widehat{y}}_i}^{(t-1)}}+f_t(x_i))+\mathsf{\Omega }(f_t)$$

(7)

After performing a second-order Taylor expansion and omitting the k − 1 rounds of pseudo-residuals, the objective function becomes the following [56]:

$$obj={\displaystyle \sum _{i=1}^n(L({\widehat{y}}^{-1}+f_t}{(x_i)}_i,y_i)+g_i{f}_i(x_i)+{\displaystyle \frac{1}{2}}{h}_i{f_t}^2(x_i))+\mathsf{\Omega }(f_t)$$

(8)

$$g_i={\displaystyle \frac{𝜕L({{\widehat{y}}_i}^{t-1})}{𝜕L({\widehat{y}}^{t-1})}}$$

(9)

$$h_i={\displaystyle \frac{{𝜕}^{2}L({{\widehat{y}}_i}^{t-1},y_i)}{𝜕L{({\widehat{y}}^{t-1})}^2}}$$

(10)

where $g_i$ and $h_i$ denote the first- and second-order derivatives.

To minimize the objective function, four hyperparameters of the XGBoost model need to be iteratively optimized: (1) max_depth; (2) the learning rate of the model (learning_rate); (3) the minimum sum of weights for node splitting (min_child_weight); and (4) the structural risk coefficient parameter. At the same time, the node splitting coefficient (gamma) of the objective function is optimized.

3.2. Snake Optimization Algorithm

The snake optimizer [39] is a recently proposed nonlinear, multi-objective population optimization algorithm inspired by snakes’ foraging, fighting, and mating behavior. The SO algorithm is divided into a food search phase, a fighting phase, and a mating phase. However, the snake optimizer has problems with slow convergence in the early stage and quickly falling into the local optimal solution, so four optimization methods are used: chaotic mapping replaces the population initialization stochastic distribution function, optimizing the food search phase of the algorithm with an average subtraction optimization strategy; a reverse learning strategy is utilized to expand the algorithm’s search range to avoid falling into the local optima; and the adaptive τ-distribution perturbation approach is used to promote the population’s iterative diversity and ensure the local search efficiency of the algorithm.

3.2.1. Food Search

Food (Q) and temperature (T) influence the behavioral patterns of snakes. Female and male snakes are uniformly distributed in a particular spatial area, and when the temperature is low and the total amount of food (Q) < 0.25, females and males randomly search for food and update their position information according to each other’s location. When the temperature (T) is less than the 0.6% threshold, females and males enter fighting or mating modes; otherwise, they focus on feeding.

At this stage, the population-position update equation becomes the following [39]:

$$X_{i,j}(t+1)=X_{food}\pm C_3\times Temp\times rand\times (X_{food}-X_{ij}(t))$$

(11)

where $X_{ij}(.)$ and $X_{food}$ indicate the position of an individual (male or female) and the best individuals, respectively, and $C_{3 }$ is a constant.

3.2.2. Combat Phase

When Temp < 0.6% threshold, males fight each other for the location of females; the fighting ability of the snake depends on the amount of food searched for in the food search phase. At this stage, due to the insufficient search range of the algorithm in the local space, there is significant differentiation in the combat capability of the male snakes, causing the algorithm to fall into the local optimum. The following is the updated formula for the optimal position of males [39]:

$$X_{i,m}(t+1)=X_{i,m}(t)+C_3\times FM\times rand\times (Q\times X_{best,f}-X_{i,m}(t))$$

(12)

where $X_{i,m}(.)$ and $X_{best.f}$ represent the position of the ith male and the best female individuals in the population, respectively, and $FM$ indicates the male snake’s fighting capability.

3.2.3. Mating Patterns

During the mating phase, the female snakes assume a dominant position. They preferentially search for locations with a cold climate and an abundance of food as the optimal breeding spots. As the snakes’ search ability is affected by both temperature and food, the snake optimizer can easily fall into the local optimal solution due to the overly small search range. The equation for the update of the female position is the following [39]:

$$X_{i,f}(t+1)=X_{i,f}(t)+C_3\times Mm\times rand\times (Q\times X_{i,f}(t)-X_{i,m}(t))$$

(13)

where $X_{i,f}(.)$ indicates the position of the ith female in the population, and $M_m$ indicates the male and female mating capabilities, respectively.

3.3. Improved Algorithms

3.3.1. Chaotic Mapping

Chaotic mapping is a complex nonlinear state with irregularity, randomness, and ergodicity, often used to generate chaotic sequences. Bernoulli chaotic mapping has the characteristic of a more uniform distribution; in this study, this approach is used instead of the population random distribution formula of the algorithm so that the population search is more efficient. Expressions such as [57]:

$$Z_{n+1}=\int \begin{array}{ll}\frac{Z_n}{1-\beta }& 0\le Z_n\le 1-\beta \\ \frac{Z_n-(1-\beta )}{\beta }& 1-\beta \le Z_n\le 1\end{array}$$

(14)

3.3.2. Average Subtraction Optimization Strategy

The SABO algorithm presents a novel computational notion called subtraction of search agents B and A, defined as follows [58]:

$$A-\nu B=\mathrm{sign}(F(A)-F(B))(A-{\nu }^{-}B)$$

(15)

${\nu }^{-}$ is an m-dimensional vector, including random numbers between [1, 2], and F(A) and F(B) describe the objective functions corresponding to search agents A and B, respectively. This strategy calculates the displacement of an individual search agent ${ X}_i$ in the search space by taking the arithmetic average of ${\nu }^{-}$ subtraction of the position information of female snakes ($X_{i,f}$) and male snakes ($X_{i,m}$). Reducing the total amount of food (Q) and the influence of temperature (T) on the population’s random search promotes the algorithm’s global search ability. So, Equation (11) becomes:

$$X_{best}(i,j)=X_i(t)+((rand(1:\dim ).*(X_i(t)-{\nu }^{-}{X}_i(X_{if},X_{i,m})))./(Nm)$$

(16)

where $X_{best}(.)$ and $X_{i,}\left(t\right)$ represent the positions of the best and ith individuals in the population at the ith iteration, respectively.

3.3.3. Reverse Learning Strategy

The reverse learning strategy is a strategy for particles to find an optimal solution at the current location for other regions [59], which can enlarge the SO algorithm’s search scope to prevent it from falling into the local optimal solution. In a given space, an individual (p) of height h is projected onto the x-axis to obtain the globally optimum position $X_{best},$ and a convex lens with focal length f is located at the origin; imaging through a convex lens gives a $P^{*}$ with a height of $h^{*}$ in the X-axis by imaging the resulting reverse optimal position ${X^{*}}_{best}$. Thus, the equation can be defined as follows [39,59]:

$${X^{*}}_{best}(t)={\displaystyle \frac{(a_j+b_j)}{2}}+{\displaystyle \frac{(a_j+b_j)}{2n}}-{\displaystyle \frac{X_{best}(t)}{n}}$$

(17)

$a_j$ and $b_j$ denote the current solution’s upper and lower bounds in the j-dimensional dimension. Thus, Equations (12) and (13) become the following [39,59]:

$$X_{best}(i,j)=Xm\left(i,j\right)+C_3*FM*rand*\left(Q*{X^{*}}_{best}\_f\left(t\right)-Xm\left(i,j\right)\right){\displaystyle \frac{1}{2}}$$

(18)

$$X_{best}\left(i,j\right)=Xm\left(i,j\right)+C_3*Mf*rand*\left(Q*{X^{*}}_{best}\_f\left(t\right)-Xf\left(i,j\right)\right){\displaystyle \frac{1}{2}}$$

(19)

3.3.4. Adaptive $\mathsf{\tau }\text{-}$Distribution Perturbation Approach

The population diversity reduces in the latter iterations of the SO. This leads to the algorithm having a tendency to become trapped in the local optima. Thus, an adaptive t-distribution dimension-by-dimension mutation approach is proposed for perturbing the individuals with the optimum fitness in the food search, combat, and mating stages, respectively, improving the SO algorithm’s capability to jump out of local optima. When $t\to \left(n\to \infty \right)\to N\left(0,1\right), t=1=C(O,1)$, where N (0, 1) and C (0, 1) describe the Gaussian and Cauchy distributions, respectively. The two boundaries of the t-distribution are Gaussian distribution and Cauchy distribution; these possesses strong local exploitation ability based on maintaining population diversity, which guarantees the population’s convergence rate in the later step. The formula transforms as follows [60]:

$${X^i}_{best}(\mathrm{i},\mathrm{j})=X_{\mathrm{best}}(i,j)+\tau (C\_iter)\times X_{best}(i,j)$$

(20)

where iter describes the current number of iterations; t(iter) describes the t-distribution with the degree of freedom t. Figure 2 presents the flowchart of the HSOA.

In this paper, four assessment indices, including the coefficient of determination (R²), mean absolute error (MAE), root mean square error (RMSE), and mean absolute percentage error (MAPE), are used to evaluate the efficiency of different models. MAE expresses the mean value of the residuals of creep flexibility and creep flexibility at the sample points. RMSE denotes the average degree of residuals at the sample points. Although R² does not maximize the residuals with the database variance, it more accurately expresses the differences in the efficacy of each model. The expressions for the four metrics are as follows [61,62]:

$$R^2={\displaystyle \frac{{\displaystyle \sum _{k=1}^N(}{q}_{o,k}-{\overline{q}}_0)(q_{t,k}-{\overline{q}}_t)}{\sqrt{{\displaystyle \sum _{K=1}^N{(q_{0,k}-{\overline{q}}_0)}^2{\displaystyle \sum _{K=1}^N(q_{,t,k}-{\displaystyle \sum {\overline{q}}_t{)}^2}}}}}}$$

(21)

$$MAE={\displaystyle \frac{1}{N}}({\displaystyle \sum _{K=1}^N|\frac{q_{0,k}-q_{t,k}}{q_{0,k}}}|)$$

(22)

$$RMSE=\sqrt{{\displaystyle \frac{1}{N}}{\displaystyle \sum _{K=1}^N(q_{0,k}-q_{t,k}}{)}^2}$$

(23)

$$MAPE={\displaystyle \frac{100\%}{N}}{\displaystyle \sum _{K=1}^N|\frac{q_t-q_0}{q_0}|}$$

(24)

where N is the total sample size, $q_{0,k}$ indicates the real value corresponding to k samples, $\overline{q_0}$ denotes the mean value, $q_{t,k}$ denotes the sample output value, $\overline{q_t}$ denotes the sample output mean value, and k = 1/N.

4. Results and Discussion

The training platform was hosted using a 12th Gen Intel(R) Core(TM) i7-12700H 2.30 GHz(TM), 16.0 GBRAM, and an NVIDIA GeForce RTX 3060 with 16G GPU. The software configuration is as follows: MATLAB 2024b and Python 3.9.

The purpose of the computational analysis lies in validating the predictive capabilities of the three models, namely HSOA-SVM, HSOA-RF, and HSOA-XGBoost. HSOA will optimize the hyperparameters of each ML model using five-fold cross-validation. The calculation process also uses cross-validation operation on 70% of the training data and 30% of the test data, and the outcome is the mean of the five-fold cross-validation.

4.1. Model Prediction Results

4.1.1. Optimal Hyperparameter Settings for ML Models

Table 2. HSOA-SVM, HSOA-RF, and HSOA-XGBoost hyperparameter settings.

HSOA-SVM		HSOA-RF		HSOA-XGBoost
Number of iterations	60	Number of iterations	60	Number of iterations	60
Population size	40	Population size	40	Population size	40
c	1.0	N_estimators	18	Max_depth	7
Kernel function coefficient	1/K	Min_leaf_nodes	8	Learning rate	0.47
Decision_function_shape	0.9	Max_depth	3.2	Min_child_weight	0.8
Penalty	0.4	Max_features	7	$\gamma$	0.85
/	/	/	/	$\lambda$	0.9

lists the hyperparameters optimized for each ML model, and the value of each parameter is the optimal value based on multiple iterations. The population quantity and the number of iterations is determined using the experimental approach.

4.1.2. Six Model Prediction Results

After determining the three models’ optimum hyperparameters using the HSOA method, 4207 data separated in advance from the database are utilized to evaluate the ML model’s predictive efficiency and generalization capability. The results of five-fold cross-validation are shown in

Table 3. Comparison of ML model performance.

ML	R2		MAE		MAPE		RMSE
	Train	Test	Train	Test	Train	Test	Train	Test
SVM	0.812	0.826	5.76	5.08	7.43	6.78	9.87	8.98
RF	0.825	0.845	5.67	5.13	7.31	6.42	9.65	8.55
XGBoost	0.877	0.849	5.35	5.36	7.23	6.25	9.45	8.54
HSOA-SVM	0.901	0.908	1.78	1.96	2.23	3.57	5.09	5.16
HSOA-RF	0.924	0.926	1.43	1.66	2.05	2.79	4.17	4.08
HSOA-XGBoost	0.945	0.968	1.26	1.45	2.24	2.33	4.01	3.88

. From the validation results, it can be seen that all six ML models attained relatively high precision, and the accuracies of the three models, namely SVM, RF, and XGBoost, on the test set witnessed increases of 9.927%, 9.58%, and 14.1%, respectively. Simultaneously, the accuracies of the optimized model on the training set and the test set exhibit only a subtle difference, indicating that the HSOA effectively reduces the overfitting of the model. Figure 3 presents the metric radar plots for the six models’ training and testing sets. The radar chart intuitively reflects the differences in four evaluation indicators between the original model and the optimized model in training set a and testing set b.

Figure 4 describes the prediction outcomes of the original models SVM (a), RF (b), and XGBoost (c), as well as the optimized models HSOA-SVM (d), HSOA-RF (e), and HSOA-XGBoost (f) on the test set. Analyzing and comparing the six models’ regression prediction outcomes in Figure 4, it is seen that the predicted values of the improved models have a smaller overall difference from the actual values, and the prediction accuracy is higher, demonstrating that the optimized model has considerable anti-overfitting ability and most of the data points have sufficient accuracy, even though there are some large discrete values in all three improved models.

4.1.3. Improved Model Performance Analysis

Residual analysis is an essential index for appraising the efficiency of diverse models. The residuals of different nonlinear regression models are closely related to the optimization algorithms and the mean and dispersion of the models themselves. The computational residuals of the three optimization models are shown in Figure 5. The particles in the residual plot of HSOA-XGBoost are more evenly distributed on both sides of the central axis, and the number of abnormal particles is significantly lower than that of HSOA-RF and HSOA-SVM. This is because, compared with the SVM and RF models, XGBoost effectively enhances the base learner due to the optimization of the $\lambda$ coefficient. Meanwhile, Figure 6 shows the three ML models’ residual interval distributions. The HSOA-XGBoost model’s residual interval distribution is more reasonable, indicating higher prediction stability.

Figure 7 Depicts the curves of the absolute errors of the three models versus the experimental values, and this curve presents an overview of the three models’ overall errors. The minimum and maximum values of the absolute errors of the three models are all between 1 and 4, respectively, indicating that all three improved models have reached a high-precision prediction level. This is also largely consistent with the residual ranges of the three models. At the same time, HSOA-XGBoost performs optimally in particle regression, residual analysis, and absolute error.

Although the HSOA-XGBoost model shows high prediction performance, the attention given to the black box mechanism behind the ML model prediction has been insufficient. Therefore, it is crucial to further explain the logic of HSOA-XGBoost.

4.2. Feature Importance Analysis Using SHAP

4.2.1. Characteristic Importance Analysis

To interpret the ML model, its output is manifested through a linear summation of its input features multiplied by the corresponding SHAP values. The exclusive solution for the SHAP value is derived by weigthing the average over all possibilities and summing all possible combinations of the eigenvalues on this basis [63].

$${\varphi }_{\mathrm{i}}={\displaystyle \sum _{S\subseteq N\backslash \left\{i\right\}}\frac{\left|S\right|!(|N|-|S|-1)!}{\left|N\right|!}}(\nu (S\cup \{i\})-\nu (S))).$$

(25)

where $\left|S\right|$ is the number of non-zero indices, ${\displaystyle \frac{\left|S\right|!\left(\left|N\right|-\left|S\right|-1\right)!}{\left|N\right|!}}$ represents the weight of an index, $\left(\nu \left(S\cup \left\{i\right\}-\nu \left(S\right)\right)\right)$ represents the marginal contribution of I to $S$, and $S\subseteq N\setminus \left\{i\right\}$ represents summing over all possible indices.

The SHAP value for feature importance is regarded as the average predicted value of the sample for that feature minus the average predicted value of the other specimens. The SHAP value for each feature is shown in Figure 8, and the features are ranked according to their average SHAP value. In the compliance flexibility experiments, the loading time is the most crucial parameter for forecasting the creep compliance of concrete, the water–cement ratio is the second most significant variable, the cement content ranks third, the compressive strength follows closely behind, and input characteristics such as loading stress and strength–stress ratio also play a vital role. Although input variables such as stress level, temperature, relative humidity, and the weight ratio of aggregate to cement are placed behind other characteristic elements, this does not mean that these characteristic elements cannot exert significant functions in the creep process.

Figure 9 shows the boundary of the feature values separating the negative and positive contributions for the five influential features. In this paper, the impact of the top five features explained by the global SHAP value on the final forecasting of creep compliance is compatible with the ranking given by feature importance, and is also consistent with the characteristic factors of creep in the curing theory. Meanwhile, the linear incorporation of SHAP values can contribute to all features. For a fixed base value, the creep flexibility’s final prediction result is the offset outcome of the SHAP values with positive and negative contributions, and for the following text force plot is used to visualize this offset process.

4.2.2. Single Factor Analysis

To gain a further understanding of how HSOA-XGBoost can deliver accurate predictions, it is essential to further analyze the SHAP values of the five influential input variables. It is evident from Figure 10 that the distribution of SHAP values exhibits a specific development tendency along with the variation in parameter values. The influence rules of the six parameters on the prediction of HSOA-XGBoost are as follows: The prolongation of the creep age and the increase in cement content and water–cement ratio can significantly increase the creep prediction value. Meanwhile, at the low stress level, the internal of the concrete specimen will cause the surface aggregation of aggregate particles due to compressive stress, which will result in the reduction in the creep effect. At a high stress level, the increase in stress will significantly increase the creep prediction result. Concurrently, high-level stress causes the local materials to present nonlinear viscoelastic behaviors, exposing concrete to the interaction between creep and fatigue, eventually leading to creep damage. A growth in compressive strength will reduce the creep prediction result. The above phenomenon is in line with the creep hydration law in consolidation theory.

4.2.3. Explanation of Local Features

The overall influence of the five most influential factors on creep prediction was addressed in the preceding text. However, as the temperature increases, the cementitious material’s Young’s modulus, ultimate yield strength, and linear expansion coefficient change. This will cause a change in the capillary water content to affect the generation and development of microcracks. The variation in the water–cement ratio will affect the viscosity and exudation amount of water in concrete. Meanwhile, the sliding of cement gel sheets between the concrete gel water-absorbing layers and draining the water layer within the cement gel sheet can have a significant effect on the essential creep deformation. Thus, for a fixed base value of 67.23, the final forecasting of creep flexibility is the outcome of the interaction of the characteristic variables. Consequently, there exists a compelling need to delve deeper into the impact exerted by diverse input variables on creep prediction.

As shown in Scenario 1 in Figure 11a, the predicted result is greater than the base value, which may be due to the higher temperature, the increase in the activation energy of water between the water molecules of the gel material, the increase in the infusibility of water, the reduction in the water content inside the gel material, and the reduction in the instantaneous elastic stiffness and delayed decreasing in the stiffness of the concrete. Under load, the ability to resist stress relaxation is reduced at 65% relative humidity, resulting in a significant escalation of the creep deformation magnitude.

In contrast to Scenario 1, the predicted outcomes presented for Scenario 2 in Figure 11b fall short of the benchmark value. The meager compressive strength along with the elevated relative humidity constitute significant negative correlative factors in the prediction of creep compliance. Under the humidity condition of RH_Test = 99%, the increase in sliding water molecules in the gel makes it more likely for interlayer sliding to occur, and simultaneously, the gel pore spacing is more prone to being transformed into the space between adjacent layers. This results in a relatively increased specific surface area of the concrete, lowering the difficulty of interlayer sliding of the matrix and diminishing the creep modulus of the concrete.

In contrast to Scenario 1 and Scenario 2, the high compressive strength in Scenario 3 (as shown in Figure 11c) offers a positive correlation for creep prediction. In Scenario 3, the water–cement ratio exhibits the highest negative correlation, and the water–cement ratio and stress are key factors in the creep of early-age concrete. Under a low initial water–cement ratio, the creep modulus increases rapidly and nonlinearly as the degree of hydration grows. This implies that concrete, which contains a relatively high water–cement ratio during the hydration process, possesses a higher creep compliance.

4.3. Exemplification of General Scenarios

To illustrate the efficiency of the improved model in predicting the concrete’s creep nature, the author selected four typical cases from the test set and predicted the creep compliance within the range of 1 to 2190 days.

Table 4. Four classical scenarios used as examples.

	w/c	a/c	C (kg/m3)	fc28 (MPa)	H (mm)	V/S	T (°C)	RH_Test (%)	Sigma (MPa)	sigma/fct0	cem	t′ (Days)
S1	0.31	4.44	405	88	600	24	21	101	17.68	0.20	SL	28
S2	0.55	5.39	336	33	800	44	20	65	7.36	0.402	RS	7
S3	0.41	5.59	332	41	1400	47	20	65	9.4	0.32	R	14
S4	0.48	5.86	325	54	600	33	54	40	5.52	0.117	SL	7

presents the detailed circumstances of the four distinct schemes (Scenario S1–Scenario S4). Scenario S1 and Scenario S2 are creep experiments commencing from the concrete ages of 28 (t₀) and 7 days (t₀), with the creep durations being 200 and 553 days, respectively. Scenario S3 and Scenario S4 are creep experiments starting at the concrete ages of 14 and 7 days, respectively, and the creep periods are 1315 and 2190 days, respectively. In four typical scenarios, the author uses the optimized XGBoost model and the B₄ model for creep prediction.

Table 5. Creep flexibility coefficient calculation formula B4.

B4 Models [25]:

$J(\widehat{\mathrm{t}},{\widehat{\mathrm{t}}}^{\prime })={\mathrm{q}}_1+R_T{C}_0(\widehat{t},{\widehat{t}}^{\prime })+C_{\mathrm{d}}(\widehat{t},{\widehat{t}}^{\prime },{\tilde{t}}_0)$ $R_T=\exp [{\displaystyle \frac{U^{\prime }}{R}}({\displaystyle \frac{1}{239}}-{\displaystyle \frac{1}{T+273}})]$ $C_0(\widehat{t},{\widehat{t}}^{\prime })=q_{2}Q(\widehat{t},{\widehat{t}}^{\prime })+q_3\ln [1+{({\displaystyle \frac{\widehat{t}-{\widehat{t}}^{\prime }}{1\mathrm{days}}})}^{0.1}]+q_4\ln ({\displaystyle \frac{\widehat{t}}{{\widehat{t}}^{\prime }}})$ $C_d(\widehat{t},{\widehat{t}}^{\prime },{\tilde{t}}_0)=q5\langle \exp [-p_{5H}H(\widehat{t},{\tilde{t}}_0)]-\exp {[-P_{5H}{H}_c({{\widehat{t}}^{\prime }}_0,{\tilde{t}}_0)]\rangle }^{0.5}$ $H(\widehat{t},{\widehat{t}}_0)=1-(1-h)\tanh \sqrt{{\displaystyle \frac{\widehat{t}-{\widehat{t}}_0}{{\tau }_{sh}}}}$ $H_c(t,t)=1-(1-h)\tanh \sqrt{{\displaystyle \frac{{{\widehat{t}}^{\prime }}_0-{\tilde{t}}_0}{{\tau }_{sh}}}}$ ${\mathrm{q}}_2={\displaystyle \frac{p_2}{1\mathrm{Gpa}}}{({\displaystyle \frac{w/c}{0.38}})}^{p_{2w}}$ $q_3=p_3{q}_2{({\displaystyle \frac{a/c}{6}})}^{p_{3a}}{({\displaystyle \frac{w/c}{0.38}})}^{p_{3w}}$ $q_5={\displaystyle \frac{p_5}{1\mathrm{Gpa}}}{({\displaystyle \frac{a/c}{6}})}^{p_{5a}}{({\displaystyle \frac{w/c}{0.38}})}^{p_{5w}}{\left|{\mathrm{k}}_h{\epsilon }_{sh\infty }({\tilde{t}}_0)\right|}^{p_{5\epsilon }}$ ${\tau }_0={\tau }_{cem}{({\displaystyle \frac{a/c}{6}})}^{p_{\tau a}}{({\displaystyle \frac{w/c}{0.38}})}^{p_{\tau w}}{({\displaystyle \frac{6.5c}{p}})}^{p_{\tau c}}$

$Q(\widehat{t},{\widehat{t}}^{\prime })=Q_f({\widehat{t}}^{\prime })[1+({\displaystyle \frac{Q_f({\widehat{t}}^{\prime })}{Z(t,t^{\prime })}}{)}^{r({\widehat{t}}^{\prime })}{]}^{-{\displaystyle \frac{1}{r({\widehat{t}}^{\prime })}}}$ $Q_f({\widehat{t}}^{\prime })=[0.086({\displaystyle \frac{{\widehat{t}}^{\prime }}{1\mathrm{days}}}{)}^{{\displaystyle \frac{2}{9}}}+1.21({\displaystyle \frac{{\widehat{t}}^{\prime }}{1\mathrm{days}}}{)}^{{\displaystyle \frac{4}{9}}}{]}^{-1}$ $Z(t,t)=({\displaystyle \frac{{\widehat{t}}^{\prime }}{1\mathrm{days}}}{)}^{-0.5}\ln [1+({\displaystyle \frac{\widehat{t}-{\widehat{t}}^{\prime }}{1\mathrm{days}}}{)}^{0.1}]$ $r({\widehat{t}}^{\prime })=1.7{({\displaystyle \frac{{\widehat{t}}^{\prime }}{1\mathrm{days}}})}^{0.12}+8$ $q_1={\displaystyle \frac{p_1}{E_{28}}}$ $E_{28}=4734\mathrm{Mpa}\sqrt{{\displaystyle \frac{{\overline{f}}_c}{\mathrm{Mpa}}}}$ $q_4={\displaystyle \frac{p_4}{1\mathrm{Gpa}}}{({\displaystyle \frac{a/c}{6}})}^{p_{4a}}{({\displaystyle \frac{w/c}{0.38}})}^{p_{4w}}$ ${\tau }_{sh}={\tau }_0{K}_{\tau a}(k_s{\displaystyle \frac{D}{1\mathrm{mm}}}{)}^2$

presents the formula of the B₄ model. This model considers the influences of relevant factors such as loading age, cement type, water–cement ratio, aggregate–cement ratio, 28-day compressive intensity, and environmental temperature when forecasting the creep nature of the concrete. Meanwhile, it has received theoretical mathematical support from nanostructure micro-stress relaxation theory, water molecule diffusion theory, and microcrack damage model. Therefore, the B₄ model can effectively predict the creep compliance changes in concrete under different time spans.

Where $\widehat{\mathrm{t}}$ and ${\widehat{\mathrm{t}}}^{\prime }$ are the current age adjusted for temperature and age at loading adjusted for temperature; $J(\widehat{\mathrm{t}},{\widehat{\mathrm{t}}}^{\prime })$ represents total creep;$q_1$ represents instantaneous compliance; $C_0(\widehat{t},{\widehat{t}}^{\prime })$ indicates the essential compliance function; and $C_{\mathrm{d}}(\widehat{t},{\widehat{t}}^{\prime },{\tilde{t}}_0)$ indicates the additional creep compliance caused by drying. $Q(\widehat{t},{\widehat{t}}^{\prime })$ indicates a binomial integral. $Q_f({\widehat{t}}^{\prime })$, $Z(t,t^{\prime })$, and $r({\widehat{t}}^{\prime })$ describe the time-dependent parameters obtained through calculation. $q_2~q_{5 }$ are parameters related to the water–cement ratio (w/c), aggregate–cement ratio (a/c), and cement type (C). $R_T$ is a variable related to creep activation energy and temperature. ${\tau }_0$ is the shrinkage factor of concrete, ${\epsilon }_{sh\infty }({\tilde{t}}_0)$ is ultimate shrinkage strain, and ${\tau }_{sh}$ is shrinkage strain. $P_1$-$P_5$, $P_{2w}$-$P_{5w}$, $P_{3a}$-$P_{5a}$, $P_{\tau a}$, $P_{\tau w}$, $P_{\tau c}$, $p_{5\epsilon }$, ${\tau }_{cem}$, and $P_{5H}$ are the coefficients related to the type of cement.

Figure 12 shows the forecast outcomes of the HSOA-XGBoost and B4 models for the experimental values of creep compliance in four scenarios. The fitting curves of HSOA-XGBoost in the four scenarios are consistent with the experimental results. In four typical scenarios, although the B4 model shows high accuracy in the early hydration process, it still fails to capture the creep compliance under long-term creep. At the same time, the creep prediction value of the B4 model in Scenario S1–Scenario S4 is generally greater than the actual value observed in the experiments. This changing tendency of the B₄ model is even more pronounced in Scenario S3 and Scenario S4 within a time span where the experimental data have relatively extensive coverage. The reason for this phenomenon is that the B4 model is inadequate in depicting the rapid changes in the physical state of the material under the influence of the hydration effect at the stage of creep, thus resulting in a larger calculated result for the parameters.

4.4. Limitations and Discussion

As presented in Figure 12, the predicted creep curve based on HSOA-XGBoost is non-smooth. This is a result of the nature of the input database. Even though each calculation for each set of experimental data in the figure is accomplished under the identical mixture parameters and experimental characteristic conditions, creep experiments persist for several months or years, so errors in the experimental outcomes are unavoidably present. Concurrently, ML models are highly dependent (R² = 0.708) on the integrity of the input data. Incomplete and scattered data points will influence the fitting effect of ML models for compliance creep test values. Although the HSOA optimization strategy has been proposed to improve the robustness and generalization ability of models such as SVM, RF, and XGBoost to reduce the influence of outliers, the data in the NU database are scattered points combined under various experimental conditions. Therefore, before training the data, interpolation processing of missing values needs to be performed using logarithmic functions to improve the integrity of the database. However, unstructured databases may pose certain challenges to prediction. Therefore, two possible measures need to be adopted to improve the applicability of ML models in creep prediction:

• Improving the quality of the database can be carried out by collecting more creep data under different experimental conditions, by considering more factors such as the shape of the test block (prism or cylinder) and the mechanical state during the experiment (three-point bending test or axial compression test) that will have an influence on creep.
• Prior to the training of an ML model, the data within the database should be preprocessed to eliminate the impacts of outliers and noise on the model training process. One prevalent approach is to substitute the missing values with the predicted values derived from a specific functional model. However, this method entails certain risks.

5. Conclusions

In this study, three machine learning models were established based on the introduction of the HSOA, and the creep behavior of concrete was predicted. Firstly, the NU database was preprocessed and divided into the training and test datasets. The five-fold cross-validation method was adopted to enhance the model’s prediction precision, and the ML prediction model’s generalization ability was verified based on 1262 sets of creep data. Subsequently, the predictions of the HSOA-XGBoost model were explicated considering the SHAP theory. The effects of diverse input variables on creep prediction were analyzed and accounted for by means of common creep theories. In the end, the efficiency of the B4 and HSOA-XGBoost models was evaluated through four scenarios. The subsequent conclusions can be stated:

• Through the optimization of the HSOA and cross-validation, all three optimized models achieved a relatively high level of accuracy. In the test dataset, HSOA-XGBoost demonstrated higher precision, with R² reaching 0.908, 0.926, and 0.968, respectively. HSOA-XGBoost exhibited a more robust performance. The fitting ability for the creep experimental values is far superior to that of the widely used B4 model.
• The SHAP theory offered a rational explanation of the ML model and provided five input features exerting a considerable influence on the prediction of concrete creep: (1) creep age; (2) loading stress; (3) cement type; (4) water–cement ratio; and (5) compressive strength. The five most impactful input features revealed by the SHAP theory were fundamentally in line with the influencing factors in the creep theory.
• As the ML model considers the balance between bias and variance, the prediction results of the HSOA-XGBoost model in long-term creep are close to the experimental observed values. Meanwhile, the creep law captured by the HSOA-XGBoost model is consistent with the general creep law in the consolidation theory, which verifies the rationality of the ML model. Given the scattered and imperfect data in the NU database, the creep compliance curve is non-smooth. Further supplementation of the database may contribute to more complete predictions. Although the HSOA optimizes the hyperparameters of models such as SVM, RF, and XGBoost in multiple iterations, due to the limitations of the prediction principle of ML models themselves and the limitations of the search principle of the SO algorithm, 100% prediction accuracy does not exist.