A Machine Learning based uncertainty quantification for compressive strength of high-performance concrete

Abstract

High performance concrete (HPC) properties depend on both its constituent materials and their interaction. This study presents a machine learning framework to quantify the effects of constituents on HPC compressive strength. We first develop a stochastic constitutive model using experimental data and subsequently employ an uncertainty quantification method to identify key parameters in relation to the compressive strength of HPC. The resultant sensitivity indices indicate that fly ash content has the strongest influence on compressive strength, followed by concrete age at test and blast surface slag content.

Graphical abstract

Keywords

uncertainty quantification / machine learning / Artificial Neural Networks / compressive strength of concrete / dependent variables

1 Introduction

High-performance concrete (HPC) has been used widely in concrete construction industry for decades. Additional cementitious materials, i.e., blast furnace slag, fly ash and chemical admixtures (e.g., superplasticizers) are incorporated into the four fundamental ingredients of conventional concrete, namely fine and coarse aggregates, Portland cement and water to make the HPC material [1]. As a consequence, a high strength, high workability, dimensional stability, and durability material has been obtained.

Mechanical properties of the concrete are paramount importance in structural design. They have significant effect on the capacity of structures, the reliability of service-proven structural systems, etc. Therefore, a comprehensive understanding of concrete design methods is necessary to obtain the requited properties which ensure structural safety. Particularly, compressive strength is one of the crucial properties of the material and it has been studied for many years. Most concrete design methods are based on experiments or semi-experiments. Although experimental or semiexperimental results are accurate, they are quite expensive from economic point of view. Besides concrete are also designed based upon empirical methods. However, as the HPC is a mixture of many ingredients with different proportions and the mixing techniques are not indicated properly, a predictive strength from the available data is not completely reliable [2]. For instance, around ten different ingredients are mixed together to design the modern concrete and more and more ingredients tends to be added to enhance material characteristics. Consequently, the adjusted properties of the resulting concrete are increased so that the empirical methods no longer can be used to design the concrete.

The concrete material is stochastic in nature and its properties are affected by different factors (e.g., the above-mentioned constituent parameters). Many experimental studies indicated that concrete strength development is affected by different constituents [3]. Although the experimental data can be used to establish the material development as a function of the constituents, the relationship between the compressive strength and additional cementitious materials (i.e., fly ash, blast furnace, etc.) is still lacked. For the sake of better designing the concrete mixture, the influence of the concrete composition on the strength should be investigated.

According to Yeh [4], effects of some ingredient parameters on the compressive strength of the HPC are considered. A traditional approach has been performed using experimental data where one constituent parameter is varied while keeping the others fixed. By applying the same technique for all parameters, the effects of constituents on the compressive strength of low-and high-strength concrete can be estimated. However, the approach does not take interactions between constituent parameters into account leading to unsatisfied results. To construct a model of accounting for the HPC strength where the interactions between the constituent parameters are considered, design of experiments has to be employed.

Metamodeling has been used widely to replace complex computational models once experimental data are significantly large. Metamodels adopting experimental training data have been used to search for the optimum proportion mixture design. Statistical models are widely used for predicting the compressive strength tests. For instance, relationship between water-cement proportion and the compressive strengths of concrete was established by Ref. [5] using a regression technique. However, there are several factors influencing the compressive strength and the compressive strength predictions are not accurate if input variable is only the water-cement proportion. Machine Learning (ML) and Artificial Neural Networks (ANN) are regularly employed in constructing meta models. Several methods such as Neural Network (NN), XGBoost, Random Forest (RF), Support Vector Machine, k-Nearest Neighbors etc. have been performed for numerical predictions or regressions [6]. For instance, Mai et al. [7] constructed various ML models such as ANN, RF and Categorical Gradient Boosting (CatBoost) to search for the best predictive model for the compressive strength of fiber reinforced self-compacting concrete. Noureldin et al. [8] proposed a framework for supervised machine learning technique to analyze the seismic performance of frame structures where soil-structure interaction was taken into account. Then, sensitivity analysis was performed to detect the key parameters affecting the soil−structure interaction.

To date, ANNs have been employed extensively to predict the compressive strength of concrete. The ANNs consists of three layers, i.e., the input layer, the hidden layers and the output layer. Each layer includes a couple of neurons and at each of the neurons there is an associated activation function. Due to the presence of the activation function in each layer, ANN is able to learn the complex relationship between the inputs and output [9]. Although multi layers allow learning the nonlinear input-output relationship, it requests tuning several parameters and hyperparameters, e.g., batch size and the number of epochs [6]. Furthermore, it is burdensome to determine the proper network structure of ANN and to define the activation function, etc.

Ziolkowski and Niedostatkiewicz [10] utilized machine learning for concrete mix design where destructive tests were performed to obtain a large number of databases. Then, the optimal ANN architecture is determined to predict the compressive strength of the concrete resulting from different mix recipes. Lin and Wu [9] built an ANN architecture and the resultant formulae with the corresponding synaptic weights and the thresholds to predict the compressive strength of concrete in accordance with the proportion of concrete mix ingredients. Experimental data are then used to verify and validate their ANN model. Öztaş et al. [11] show applicability of NN in designing the compressive strength and slump of high strength concrete (HSC). Input variables such as the water to binder proportion, water content, fine aggregate proportion, fly ash content and superplasticizer etc. are fed into the NN model to predict the compressive strength and the slump of HSC. Other approaches, which are robust against noisy data such as Gaussian Process, polynomial, Penalized Spline regression models etc. are regularly used.

Polynomial regressions and ANN were employed by Yeh [1,3,4] to predict the compressive strength of the HPC and to assess qualitatively the effects of concrete constituents on this response. Although, good performance of the ANN models was indicated, overfitting could be ignored as the author is only based on coefficients of determination (COD-R^2) or root-mean-square error (RMSE). It is shown that the proposed approach in this study does not quantify the effects of the mixture constituents on the compressive strength. Furthermore, the constituents are dependent on each other and their interaction effects of the compressive strength have not been tackled so far. Therefore, global sensitivity analysis (SA) can be an appropriate tool for assessment of the meta models.

Global SA based on the Analysis of variances (ANOVA) decomposition (i.e., the functional analysis of variance decomposition) was introduced by Ref. [12]. The total sensitivity measurements were then proposed by Ref. [13]. Among sensitivity analysis methods, variance-based methods have been examined by numerous studies and they have been considered as versatile and effective ones. The variance-based methods known as the computer experiment are replacement tools for the experimental design’s analysis. If Monte Carlo simulations are used the methods can thus take the whole factor distribution into account [14].

This study aims to provide a framework consisting of: 1) making random samples for dependent variables, 2) building effective mathematical metamodels that approximate the experimental data, and 3) performing SA in a global context. To this end, in Section 2 the experimental data for the HPC is approximated using metamodeling (i.e., penalized spline methodology). Then, Latin hypercube sampling (LHS) described briefly in Section 3 is used to generate random inputs. Transformation of independent input variables to dependent ones is then performed to prepare sample points. Section 4 presents different metamodels. An improvement of variance-based method presented in Section 5 is finally employed to estimate effects of the dependent inputs on the model output quantitatively followed by Section 6 that explain sequential steps of uncertainty quantification procedure. In Section 7 numerical results, i.e., the COD and the sensitivity indices are discussed before we make remarks in conclusion section.

2 Data description

Experimental data (around 1000 mixtures) collected from different sources by Yeh [1,15] were examined in this study. These HPC mixtures, made with ordinary Portland cement, fly ash, blast furnace slag, and superplasticizer were assembled into a data set. The HPC specimens were compacted on a vibrating table, with moist curing maintained for 24 h. After demolding, the concrete was cured in water at temperature of 20 °C (see description in Ref. [4]). According to Ref. [1], specimens with different sizes and shapes were used in different tests. However, based on guidelines for experimental design, all compressive strengths measurements were adjusted to their equivalent values for 150 mm cylinders.

As mentioned above, the HPC is composed of different constituents. Each constituent is represented by a random input variable in the predictive mechanical model, as the input values are of different magnitudes. In this study, the compressive strength of the HPC is considered as a function of eight variables, i.e., 1) cement (\mathrm{k}\mathrm{g}/{\mathrm{m}}^3); 2) fly ash (\mathrm{k}\mathrm{g}/{\mathrm{m}}^3); 3) blast furnace slag (\mathrm{k}\mathrm{g}/{\mathrm{m}}^3); 4) water (\mathrm{k}\mathrm{g}/{\mathrm{m}}^3); 5) superplasticizer (\mathrm{k}\mathrm{g}/{\mathrm{m}}^3); 6) coarse aggregate (\mathrm{k}\mathrm{g}/{\mathrm{m}}^3); 7) fine aggregate (\mathrm{k}\mathrm{g}/{\mathrm{m}}^3); 8) age at test (d). Tab.1 summarizes the histograms, corresponding fitted distributions, and ranges of the variables within a stochastic modeling framework.

3 Sampling method

LHS, due to its simplicity and effectiveness, is used to generate quasi-random variables for high-dimensional problems. As the training points provided by this method are scattered entirely throughout the design space, the computational effort is then significantly reduced [16]. The generated independent variables are then transformed into dependent variables using Gaussian copula and Cholesky factorization. More details of this method can be found in Ref. [17].

4 Metamodels

4.1 Penalized spline regression

Recalling the n-dimensional model function y=f\left(\mathit{x}\right) with \mathit{x}=\left[x_1,x_2,...,x_n\right] being the vector of n independent variables. Correspondingly, if we define the n-component smoothing vector \mathit{\lambda }=\left[{\lambda }_1,...,{\lambda }_n\right], at a design point \left({\mathit{x}}^{\left(j\right)},y^{\left(j\right)}\right) the additive model is rewritten as follows:

(1)y^{\left(j\right)}={\beta }_0+\sum _{i=1}^n{f}_i\left(x_i^{\left(j\right)}\right)+{ϵ}_j,

where x_i^{\left(j\right)},j=1,...,N, with N being the number of realizations, denotes the jth realization of the ith variable. Defining {\mathbf{\Xi }}_i knots for the ith variable and choosing the degree of the truncated power basis of p. Correspondingly, if we introduce vectors of coefficients \mathit{\beta } and \mathit{\zeta } with their respective ith rows {\mathit{\beta }}_i=\left[{\beta }_{1i},...,{\beta }_{pi}\right] and {\mathit{\zeta }}_i=\left[{\zeta }_{1i},...,{\zeta }_{{\mathrm{\Xi }}_{i}i}\right] the mixed mode representation can be written as:

(2)f\left(\mathit{x},\mathit{\beta },\mathit{\zeta }\right)={\beta }_0+\sum _{i=1}^n\left({\beta }_{1i}{x}_i+\cdots +{\beta }_{pi}{x}_i^p+\sum _{l=1}^{{\mathrm{\Xi }}_i}{\zeta }_{li}{\left(x_i-{\xi }_{li}\right)}_{+}^p\right).

Let define \mathit{X} and \mathit{Z} as linear spline bases with their respective ith rows {\mathit{X}}_i=\left[1x_i,...,x_i^p\right] and {\mathit{Z}}_i=[{\left(x_i-{\xi }_{1i}\right)}_{+}^p{\left(x_i-{\xi }_{{\mathrm{\Xi }}_{i}1}\right)}_{+}^p]. Equation (2) is rewritten in terms of matrix form as:

(3)\mathit{y}=\mathit{X}\mathit{\beta }+\mathit{Z}\mathit{\zeta }+\mathit{\epsilon },\mathrm{C}\mathrm{o}\mathrm{v}\left(\mathit{\epsilon }\right)={\sigma }_{\epsilon }\mathit{I}.

Consequently, the approximated model is shown as:

(4)\hat{\mathit{y}}={\beta }_0+\sum _{i=1}^n\left({\beta }_{1i}{x}_i+...+{\beta }_{pi}{x}_i^p+\sum _{l=1}^{{\mathrm{\Xi }}_i}{\zeta }_{li}{\left(x_i-{\xi }_{li}\right)}_{+}^p\right),

(5)\hat{\mathit{y}}=\mathit{X}\hat{\mathit{\beta }}+\mathit{Z}\hat{\mathit{\gamma }}=\mathit{C}{\left({\mathit{C}}^{\mathrm{T}}\mathit{C}+\mathbf{\Lambda }\right)}^{-1}{\mathit{C}}^{\mathrm{T}}\mathit{y},

with

(6)\mathit{C}\equiv \left[\mathit{X}|\mathit{Z}\right]\mathrm{a}\mathrm{n}\mathrm{d}\mathbf{\Lambda }=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(0,0,0,{\lambda }_1^2{1}_{{\mathrm{\Xi }}_1\times 1},...,{\lambda }_n^2{1}_{{\mathrm{\Xi }}_n\times 1}\right).

The error between the mechanical model (experimental data) and the approximated regression spline model is then minimized subjected to the constraints on the coefficients \mathit{\beta } and \mathit{\zeta }, see Ref. [18] for details. Using Lagrange multiplier method with the penalty function \lambda (\cdot ) (i.e., smoothing parameter) for the ith predictor, the coefficients \hat{\mathit{\beta }} and \hat{\mathit{\zeta }} are defined as the minimizer of

(7)\sum _{j=1}^N\{y^{\left(j\right)}-f\left({\mathit{x}}^{\left(j\right)},\mathit{\beta },\mathit{\zeta }\right){\}}^2+\sum _{i=1}^n{\lambda }_i^2\left({\xi }_{l,i}\right){\zeta }_{l,i}^2.

Generalized cross-validation (GCV) is minimized to obtain the smoothing parameter {\hat{\lambda }}_{\mathrm{G}\mathrm{C}\mathrm{V}}, see Ref. [19] for details.

4.2 Artificial Neural Networks

ANNs are designed to model complex, nonlinear relationships in data. An ANN consists of layers of neurons that transform inputs using weights, biases, and activation functions. For a neuron in layer l, the output is computed as:

(8){\mathit{\text{z}}}_j^{\left(l\right)}=\sum _{i=1}^{n^{\left(l-1\right)}}{w}_{ij}^{\left(l\right)}{a}_i^{\left(l-1\right)}+b_j^{\left(l\right)},a_j^{\left(l\right)}=\varphi \left({\mathit{\text{z}}}_j^{\left(l\right)}\right),

where w_{ij}^{\left(l\right)} is the weight connecting neuron i in layer l-1 to neuron j in layer l, b_j^{\left(l\right)} is the bias term, \varphi \left(\mathit{\text{z}}\right) is the activation function (e.g., ReLU: \varphi \left(\mathit{\text{z}}\right)=\mathrm{m}\mathrm{a}\mathrm{x}\left(0,\mathit{\text{z}}\right)).

Learning occurs through backpropagation, an iterative optimization process that minimizes a loss function, such as Mean Squared Error (MSE) in our case:

(9)\mathcal{L}=\frac{1}{m}\sum _{i=1}^m{\left(y_i-\widehat{y_i}\right)}^2,

where m is the number of training samples, y_i the actual value, and \widehat{y_i} the predicted value. Backpropagation computes gradients of the loss with respect to weights and biases using the chain rule. The weight update is performed as:

(10)w_{ij}^{\left(l\right)}\leftarrow w_{ij}^{\left(l\right)}-\eta \frac{\mathrm{\partial }\mathcal{L}}{\mathrm{\partial }{w}_{ij}^{\left(l\right)}},b_j^{\left(l\right)}\leftarrow b_j^{\left(l\right)}-\eta \frac{\mathrm{\partial }\mathcal{L}}{\mathrm{\partial }{b}_j^{\left(l\right)}},

where \eta is the learning rate. Multiple iterations are required to progressively reduce the loss, as the optimization process starts far from the optimal solution.

To prevent overfitting, L2 regularization adds a penalty term \lambda \sum w^2 to the loss function, while dropout randomly deactivates neurons during training to enhance generalization. In this study, the ANN architecture included three hidden layers (400-400-300 neurons) with ReLU activations. The model was trained using MSE as the loss function, and early stopping was employed to halt training once validation performance stopped improving.

4.3 XGBoost

XGBoost is a machine learning algorithm based on gradient-boosted decision trees. It builds models iteratively, where each tree corrects the residual errors of the previous model. The prediction at iteration k is given by:

(11)\widehat{y_i^{\left(k\right)}}=\widehat{y_i^{\left(k-1\right)}}+f_k\left(x_i\right),

where f_k\left(x_i\right) is the kth tree, trained to minimize the residual error. The overall objective function is:

(12)\mathcal{L}=\sum _{i=1}^{m}l\left(y_i,\widehat{y_i}\right)+\sum _{k=1}^K\mathrm{\Omega }\left(f_k\right),

where l\left(y_i,\widehat{y_i}\right) is the loss function (e.g., squared error), and \mathrm{\Omega }\left(f_k\right)=\gamma T+{\displaystyle \frac{1}{2}}\lambda |w{|}^2 is the regularization term, penalizes large trees with many leaf nodes (T) or high weights |w{|}^2 representing increased model complexity. Here, T represents the number of leaf nodes, and |w{|}^2 is the L2 regularization term for leaf weights. XGBoost incorporates additional optimizations such as parallel processing, sparsity awareness, and tree pruning to enhance computational speed and accuracy. In this study, Bayesian optimization was used to fine-tune hyperparameters like learning rate, tree depth, and regularization terms. This approach ensured efficient modeling of nonlinear interactions among HPC constituents while maintaining strong generalization to test data.

ANNs are adept at capturing high-dimensional, nonlinear relationships through dense neuron connections, while XGBoost handles feature interactions and residual errors efficiently with interpretable tree-based models. Due to these differing strengths, these approaches have proven effective for accurate prediction of HPC compressive strength.

5 Sensitivity analysis

5.1 The ANOVA decomposition

Let consider a mathematical model y=f\left(\mathit{x}\right)\in {\mathbb{L}}^2 with \mathit{x}={[x}_1,x_2,...,x_n] being in a unit hypercube domain. We can then expand f\left(\mathit{x}\right) as follows:

(13)f(\mathit{x})=f_0+\sum _{i=1}^n{f}_i\left(x_i\right)+\sum _{j>i}^n{f}_{ij}\left(x_i,x_j\right)+f_{\mathrm{1...}n}\left(x_1,...,x_n\right),

with {\mathit{x}}_{\mathit{v}}=\left\{x_1,\dots ,x_m\right\}, E\left[f_{i_1,\dots ,i_m}\left({\mathit{x}}_{\mathbf{v}}\right)|{\mathit{x}}_{\sim i_{\mathrm{q}}}\right]=0, i_{\mathrm{q}}=i_1,\dots ,i_m and {\mathit{x}}_{\sim i_{\mathrm{q}}}={\mathit{x}}_{\mathbf{v}}\left\{x_{i_{\mathrm{q}}}\right\}. Note that the functions are orthogonal yielding E\left[f_i\times f_j\right]=V_i{\delta }_{ij} with {\delta }_{ij}=\{\begin{array}{c}1,\mathrm{i}\mathrm{f}i=j,\\ 0,\mathrm{i}\mathrm{f}i\ne j,\end{array} see Refs. [14,20].

(14)\begin{array}{rl}& f_0=E\left[y\right],\\ & f_i=E_{\mathit{x}\sim i}\left[y|x_i\right]-E\left[y\right],\\ & f_{ij}=E_{\mathit{x}\sim ij}\left[y|x_i,x_j\right]-f_i-f_j-E\left[y\right].\end{array}

The higher order terms are written in the same manner.

(15)\begin{array}{rl}{V}_i=& V\left[f_i\left(x_i\right)\right]=V_{x_i}\left[E_{\mathit{x}\sim i}\left[y|x_i\right]-E\left[y\right]\right]\\ V_{ij}=& V\left[f_{ij}\left(x_i,x_j\right)\right]=V_{x_i,x_j}\left[E_{\mathit{x}\sim ij}\left[y|x_i,x_j\right]\right]\\ & -V_{x_i}\left[E_{\mathit{x}\sim i}\left[y|x_i\right]\right]-V_{x_j}\left[E_{\mathit{x}\sim j}\left[y|x_j\right]\right]\end{array}.

As the inputs are independent, Eq. (13) is unique. The total variance of f\left(\mathit{x}\right) is then derived by

(16)V_y=\sum _{i=1}^n{V}_i+\sum _{j>i}^n{V}_{ij}+\cdots +V_{1,...,n},

where V_i denotes the partial variance of fdue to x_i alone and V_{i_1,...,i_m} known as the partial variance of due to the interactions between x_i and x_j. If we divide both sides of Eq. (16) by V_y, that yields

(17)\sum _{i=1}^n{S}_i+\sum _{j>i}^n{S}_{ij}+\cdots +S_{1,...,n}=1.

The Sobol’s sensitivity index is obtained consequently

(18)S_{i_1...i_m}=\frac{V_{i_1,...,i_m}}{V_y},

where S_i=V_i/V_y estimates the first-order effect of x_i on y and S_{i_1,...,i_m} estimates the high-order effect of joint contribution of x_1,...,x_m on y. The total effect of all input variables {\mathit{x}}_{\sim i} on y reads

(19)S{T}_i=S_i+\sum _{j\ne i}^n{S}_{ij}+\cdots +S_{1,...,n}.

It is indicated by Ref. [14], the unicity of the ANOVA decomposition is ensured only if the input variables are independent. In other words, the Sobol’s sensitivity index can be affected if the input variables depend on one another. Extended formulas of the variance-based sensitivity indices proposed by Ref. [21] will be presented in the following sections where the main and total effects of the dependent input variables on the model response are considered.

5.2 Formulas for the first-order and the total indices

If we split the input variables \mathit{x}=\left(x_1,...,x_n\right) into 2 subsets arbitrarily, i.e., the subset \mathit{w}=\left[x_{i_1},...,x_{i_l}\right] and its complementary one \mathit{z}=\left[x_{i_{l+1}},...,x_{i_n}\right], the total variance of f\left(\mathit{x}\right) in Eq. (16) can be rewritten as:

(20)V=V_{\mathrm{w}}\left[E_{\mathit{\text{z}}}f\left(w,\widehat{\mathit{\text{z}}}\right)\right]+E_{\mathrm{w}}\left[V_{\mathrm{w}}\left(f\left(w,\widehat{\mathit{\text{z}}}\right)\right)\right],

where the subset \mathit{\text{z}} is sampled using the joint probability density function (PDF) p\left(x,y\right) while \widehat{\mathit{\text{z}}} is sampled using the conditional PDF p\left(y,\widehat{\mathit{\text{z}}}|y\right). The first-order and the total-effect sensitivity indices for the subsets are described by Ref. [21] as follows:

(21)\begin{array}{rl}& S_{\mathrm{w}}=\frac{V_{\mathrm{w}}\left[E_{\mathit{\text{z}}}\left(f\left(w,\widehat{\mathit{\text{z}}}\right)\right)\right]}{V},\\ & S{T}_{\mathit{\text{z}}}=\frac{E_{\mathit{\text{z}}}\left[V_{\mathrm{y}}\left(f\left(\hat{w},\mathit{\text{z}}\right)\right)\right]}{V},\\ & S{T}_{\mathrm{w}}=\frac{E_{\mathrm{w}}\left[V_{\mathit{\text{z}}}\left(f\left(\hat{w},\mathit{\text{z}}\right)\right)\right]}{V}.\end{array}

The first-order index is estimated numerically by:

(22)S_{\mathrm{w}}=\frac{1}{V}\left[{\int }_{{\mathbb{R}}^l}^{}p\left(w\right)\mathrm{d}w{\left[{\int }_{{\mathbb{R}}^{n-1}}^{}f\left(w,\widehat{\mathit{\text{z}}}\right)p\left(w,\widehat{\mathit{\text{z}}}|w\right)\mathrm{d}\widehat{\mathit{\text{z}}}\right]}^2-f_0^2\right],

with f_0={\displaystyle \frac{1}{N}}{\displaystyle \sum _{i=1}^{N}f\left({\mathit{x}}_i\right)} and p\left(w\right) denoting a probability density function (PDF). It can be expanded further as follows

(23)\begin{array}{rl}{S}_{\mathrm{w}}=& \frac{1}{V}[{\int }_{{\mathbb{R}}^l}p\left(w\right)\mathrm{d}w[{\int }_{{\mathbb{R}}^{n-1}}f\left(w,\widehat{\mathit{\text{z}}}\right)p\left(w,\widehat{\mathit{\text{z}}}|w\right)\mathrm{d}\widehat{\mathit{\text{z}}}\\ & \cdot {\int }_{{\mathbb{R}}^{n-1}}f\left(w,{\widehat{\mathit{\text{z}}}}^{\prime }\right)p\left(w,{\widehat{\mathit{\text{z}}}}^{\prime }|y\right)\mathrm{d}{\widehat{\mathit{\text{z}}}}^{\prime }]-f_0^2],\end{array}

where the combined vectors \left(w,\mathit{\text{z}}\right), \left(w^{\prime },{\mathit{\text{z}}}^{\prime }\right), and \left(w,\mathit{\text{z}}\right), \left(w,{\check{\mathit{\text{z}}}}^{\prime }\right) are generated using the joint PDF p\left(w,\mathit{\text{z}}\right) along with the conditional PDF p\left(w,{\check{\mathit{\text{z}}}}^{\prime }|w\right). The first-order indices shown by Eq. (23) can be then obtained by

(24)\begin{array}{rl}{S}_{\mathrm{w}}=& \frac{1}{V}[{\int }_{{\mathbb{R}}^l}f\left(w,\mathit{\text{z}}\right)p\left(w,\mathit{\text{z}}\right)\mathrm{d}y\mathrm{d}\mathit{\text{z}}[{\int }_{{\mathbb{R}}^{n-l}}f\left(w,{\check{\mathit{\text{z}}}}^{\prime }\right)p\left(w,\check{\mathit{\text{z}}}|w\right)\mathrm{d}\check{\mathit{\text{z}}}\\ & \cdot {\int }_{R^{n-1}}^{}f\left(w^{\prime },{\mathit{\text{z}}}^{\prime }\right)p\left(w,{\mathit{\text{z}}}^{\prime }\right)\mathrm{d}w\text{'}\mathrm{d}{\mathit{\text{z}}}^{\prime }]].\end{array}

From both Eqs. (20), (21), and (24), {ST}_{\mathrm{y}} can be expressed by

(25)\begin{array}{rl}{ST}_{\mathrm{w}}=& \frac{1}{2V}{\int }_{{\mathbb{R}}^{n-1}}{\left[f\left(w,\mathit{\text{z}}\right)-f(\check{w},\mathit{\text{z}})\right]}^2\\ & \cdot fp\left(w,\mathit{\text{z}}\right)p\left(\check{w},\mathit{\text{z}}|\mathit{\text{z}}\right)\mathrm{d}w\mathrm{d}\check{w}\mathrm{d}\mathit{\text{z}}.\end{array}

For numerical computation, realizations of the variables \mathit{\text{z}}, {\mathit{\text{z}}}^{\prime }, and \widehat{\mathit{\text{z}}} have to be generated and the details can be found in Ref. [17].

5.3 Monte Carlo estimates

Monte Carlo (MC) algorithm is used in practice to estimate the first-order and the total-effect indices numerically. Introducing the subsets {\mathit{\text{z}}}^{\prime }, {\widehat{\mathit{\text{z}}}}^{\prime }, and {\hat{w}}^{\prime } generated from the conditional PDFs p\left(w,{\mathit{\text{z}}}^{\prime }|y\right), p\left(w,{\widehat{\mathit{\text{z}}}}^{\prime }|w\right) and p\left({\hat{w}}^{\prime },\mathit{\text{z}}|\mathit{\text{z}}\right), respectively, the first-order index reads

(26)S_{\mathrm{w}}=\frac{{\displaystyle \frac{1}{N}}{\displaystyle \sum _{j=1}^N\left(f\left(w_i,{\mathit{\text{z}}}_i\right)f\left(w_i,{\widehat{\mathit{\text{z}}}}_i^{\prime }\right)\right)-{\left({\displaystyle \frac{1}{N}}{\displaystyle \sum _{j=1}^{N}f\left(w_i,{\mathit{\text{z}}}_i\right)}\right)}^2}}{V},

where the total variance V is approximated numerically by

(27)V=\frac{1}{N_{w\mathit{\text{z}}}}\sum _{l=1}^{N_{y\mathit{\text{z}}}}{\left(f\left(w_l,{\mathit{\text{z}}}_l\right)\right)}^2-f_0^2,

with

(28)f_0=\frac{1}{N_{w\mathit{\text{z}}}}\sum _{l=1}^{N_{y\mathit{\text{z}}}}f\left(w_l,{\mathit{\text{z}}}_l\right),

being evaluated at N_{w\mathit{\text{z}}} data points of the combined vector (w,\mathit{\text{z}}). An estimation of the total-effect index yields

(29)S{T}_{\mathrm{w}}=\frac{{\displaystyle \frac{1}{N}}{\displaystyle \sum _{j=1}^N{\left(f\left(w_j,{\mathit{\text{z}}}_j\right)-f\left({\hat{w}}_j^{\prime },{\mathit{\text{z}}}_j\right)\right)}^2}}{2V}.

Numerically, the MC estimation for f\left(w,\mathit{\text{z}}\right), f\left(w,{\check{\mathit{\text{z}}}}^{\prime }\right), f\left(\check{w},{\mathit{\text{z}}}^{\prime }\right), f\left(w^{\prime },{\mathit{\text{z}}}^{\prime }\right), and f\left(w,\widehat{\mathit{\text{z}}}\right) at each data point is estimated using the above-mentioned meta model. It should be noted that n-dimensional multivariate normal distribution and multivariate conditional distribution are thus used to generate dependent variables through MC sampling algorithm. More details of the sampling method and the MC estimates of the sensitivity indices can be found in Ref. [17]. A schematic flow chart of the sensitivity analysis procedure shows the sequential steps linking experimental data preparation, surrogate model construction using ANN/ML, and sensitivity estimation.

6 Flow chart

The pipeline shown in Fig.1 for predicting HPC compressive strength involves the following steps: 1) data preprocessing, 2) implementation of machine learning models (ANN/XGBoost or Penalized Spline regression), hyperparameter tuning, 3) training, evaluation, and visualization, 4) estimating the sensitivity indices.

7 Numerical results

In this section, the proposed machine learning framework is applied to the experimental data of the HPC described in Section 2 to measure the SA indices. The input variables are dependent and the correlation matrix \mathit{C} between the input variables is shown by:

\mathit{C}=\left[\begin{array}{cccccccc}1& -0.28& -0.40& -0.08& 0.09& -0.11& -0.22& 0.08\\ -0.28& 1& -0.32& 0.11& 0.04& -0.28& -0.28& -0.04\\ -0.40& -0.32& 1& -0.26& 0.38& -0.01& 0.08& -0.15\\ -0.08& 0.11& -0.26& 1& -0.66& -0.18& -0.45& 0.28\\ 0.09& 0.04& 0.38& -0.66& 1& -0.27& 0.22& 0.19\\ -0.11& -0.28& -0.01& -0.18& -0.27& 1& -0.18& 0.00\\ -0.22& -0.28& 0.08& -0.45& 0.22& -0.18& 1& -0.16\\ 0.08& -0.04& -0.15& 0.28& -0.19& 0.00& -0.16& 1\end{array}\right].

Note that the components of the correlation matrix are computed using Pearson correlation formula, see Ref. [22] or details. The correlation among the input variables is illustrated in Fig.2. From the correlation matrix and the plot, we realize that a pair of inputs x_4 and x_5 (C_{45}=-0.66) and another one x_4 and x_7 (C_{47}=-0.45) are strongly correlated.

Scatter plots of the compressive strength with respect to (w.r.t.) each individual input variable is plotted in Fig.3. At a glance, we realize that the fly ash (x_2) and the age at test (x_8) may affect the compressive strength of the HPC the most. The effect of the input variables on the compressive strength will be quantified in the sequel.

In addition to the regression models discussed in the manuscript, we implemented two additional machine learning approaches, an ANN and XGBoost, to enhance prediction accuracy of HPC compressive strength. The ANN consists of three hidden layers with 400, 400, and 300 nodes, respectively, using ReLU activation functions. We included 30% dropout between layers and L2 regularization to prevent overfitting. The XGBoost model’s hyperparameters were optimized using Bayesian optimization, resulting in a learning rate of 0.159, maximum tree depth of 4 and 202 estimators, and regularization parameters α = 0.59 and λ = 0.046, see Tab.2 for details. Both models were evaluated by splitting the data into training and validation sets (0.3) to ensure generalizability on out-of-distribution data. XGBoost outperformed the ANN on the validation data set, achieving a lower RMSE and higher R^2 score, demonstrating better generalization capability for predicting concrete strength. These results suggest both models, particularly XGBoost, are effective at accurately mapping how multiple concrete ingredients—particularly the interactions between cement, fly ash, and blast furnace slag—influence the final compressive strength. For more theoretical and implementation details of each model, see Refs. [23,24].

Both models were evaluated on training and testing sets using metrics computed with scikit-learn. Metrics Used: mean_squared_error for MSE and RMSE; r2_score for R²; mean_absolute_error for MAE. Tools and Software Used: Data Preprocessing: Pandas (1.5.3), NumPy (1.24.4), and scikit-learn (1.2.2); ANN Implementation: TensorFlow (2.13.0) and Keras (2.13.1); XGBoost Implementation: xgboost (1.7.6); Hyperparameter Tuning: bayes_opt (1.4.2) for Bayesian optimization (XGBoost) and manual tuning with best practices (ANN); Visualization: Matplotlib (3.7.1) and Seaborn (0.12.2) for plotting parity plots, learning curves, and feature importance. History convergences of the loss functions over the number of epoch iterations are shown in Fig.4.The obtained results are shown in Tab.3 as follows.

Tab.3 and Fig.5(a) and Fig.5(b) show that the results obtained by the ANN, particularly R_{\mathrm{A}\mathrm{N}\mathrm{N}}^2=0.95\left(0.88\right) and RMSE_ANN = 3.72 (5.66) are in good agreement with the previous results shown in Ref. [4]: R_{\mathrm{R}\mathrm{E}\mathrm{F}}^2=0.94\left(0.93\right) and RMSE_REF = 3.01 (4.32) for training and testing data, respectively. Slightly discrepancies between the results are due to different hyperparameters used in our and their models.

To avoid carrying out a large number of experiments requested by the SA, metamodeling is employed in this study. As reported by Ref. [1], a quadratic regression model fails to approximate nonlinear relationships in the data (i.e., ). Thus, both quadratic [25,26] and penalized spline regression models, shown in Section 4, are constructed to approximate the experimental data on which the SA will be performed. The approximation ability of the two methods, indicated by either the COD (R^2) or the RMSE [27], is shown in Tab.4. The quadratic regression model, with , clearly fails to approximate the data, while the penalized spline regression model, with R_{\mathrm{s}\mathrm{p}\mathrm{l}}^2=0.92, fits the data almost perfectly. Therefore, the penalized spline regression model is employed in this study. The training points and the approximated regression model are plotted w.r.t. the two variables x_2, x_3 and the ones x_2, x_8 in Fig.6(a) and Fig.6(b), respectively.

The SA is then performed to measure the sensitivity indices of the input variables on the compressive strength. First, a large number of samples of dependent input variables are generated using the technique presented in Section 3. To this end, random variables following n-dimensional multivariate normal distribution and conditional normal distribution are generated using MC algorithms. Next, the predicted output values at the samples are obtained using the penalized spline regression model. Finally, the MC estimates are used to evaluate the sensitivity indices. The first-order and the total-effect sensitivity indices for the cases of independent inputs and dependent inputs are shown in Fig.7. The results demonstrate that the sensitivity indices differ when dependencies among the inputs are considered versus when they are ignored. In other words, interactions among the input variables have to be taken into account when measuring their effects on the compressive strength. Note that since the first-order and the total-effect indices here are approximated numerically, they are reduced by the R^2 to account for how well the metamodel explains the experimental data, that is

\begin{array}{rl}& {\hat{S}}_i=R^2{S}_i,\\ & {\widehat{ST}}_i=R^{2}S{T}_i.\end{array}

When measuring the sensitivity index of each individual variable on the compressive strength without considering the correlations (interactions) among them, the most sensitive variables are the coarse aggregate, followed by the age at test and the cement. However, when considering correlations among the input variables, the key influencing variables are fly ash and age at test, as shown in Tab.5. The effect of cement on the compressive strength is reduced relatively when considering the dependencies (interactions) among the concrete constituents. This can be explained from practical point of view: when the fly ash is added into the concrete, the amount of the cement in concrete can be reduced [28]. Furthermore, fly ash particles reacting with the cement can produce additional cementitious material, resulting in an increase in long-term compressive strength of the HPC. This suggests that the fly ash plays a more important role for the compressive strength when it interacts the cement matrix.

8 Conclusions

A machine learning methodology is presented in this study to quantify the effects of concrete constituents on the compressive strength of the HPC. The dependencies (interactions) between input variables, represented by the correlation matrix, are taken into account. The numerical example indicates that the effect of one input on the compressive strength is contributed not only by itself but also by its interactions with other variables. In addition to sensitivity measurement, a stochastic constitutive law representing the relationship between the compressive strength and the eight input variables has been constructed using a penalized spline regression model within a stochastic modeling framework. Consequently, the stochastic nature of the input variables has been included to account for the uncertain mechanisms in both virtual simulations and practical experiments.

References

[1] Yeh I C. Design of high-performance concrete mixture using neural networks and nonlinear programming. Journal of Computing in Civil Engineering, 1999, 13(1): 36–42

[2] Kasperkiewicz J, Racz J, Dubrawski A. HPC strength prediction using artificial neural network. Journal of Computing in Civil Engineering, 1995, 9(4): 279–284

[3] Yeh I C. Modeling of strength of high-performance concrete using artificial neural networks. Cement and Concrete Research, 1998, 28(12): 1797–1808

[4] Yeh I C. Analysis of strength of concrete using design of experiments and neural networks. Journal of Materials in Civil Engineering, 2006, 18(4): 597–604

[5] AbramsD. Design of Concrete Mixtures. Chicago, IL: Structural Materials Laboratory, Lewis Institute, 1918

[6] Urbas U, Zorko D, Vukašinović N. Machine learning based nominal root stress calculation model for gears with a progressive curved path of contact. Mechanism and Machine Theory, 2021, 165: 104430

[7] Mai H V T, Nguyen M H, Trinh S H, Ly H B. Optimization of machine learning models for predicting the compressive strength of fiber-reinforced self-compacting concrete. Frontiers of Structural and Civil Engineering, 2023, 17(2): 284–305

[8] Noureldin M, Ali T, Kim J. Machine learning-based seismic assessment of framed structures with soil–structure interaction. Frontiers of Structural and Civil Engineering, 2023, 17(2): 205–223

[9] Lin C J, Wu N J. An ANN model for predicting the compressive strength of concrete. Applied Sciences, 2021, 11(9): 3798

[10] Ziolkowski P, Niedostatkiewicz M. Machine learning techniques in concrete mix design. Materials, 2019, 12(8): 1256

[11] Öztaş A, Pala M, Özbay E, Kanca E, Çagˇlar N, Bhatti M A. Predicting the compressive strength and slump of high strength concrete using neural network. Construction and Building Materials, 2006, 20(9): 769–775

[12] I.M.Sobol’. Sensitivity analysis for non-linear mathematical models, mathematical modelling and computational experiment. Mathematical Modeling & Computational Experiment, 1993, 4: 407–411

[13] Saltelli A, Homma T. Importance measures in global sensitivity analysis of model output. Reliability Engineering and System Safety, 1996, 552: 1–17

[14] Saltelli A, Annoni P, Azzini I, Campolongo F, Ratto M, Tarantola S. Variance based sensitivity analysis of model output. Design and estimator for the total sensitivity index. Computer Physics Communications, 2010, 181(2): 259–270

[15] YehI C. Concrete Compressive Strength. 2007 (available at the website of UC Irvine Machine Learning Repository)

[16] Stein M. Large sample properties of simulations using latin hypercube sampling. Technometrics, 1987, 29(2): 143–151

[17] Vu-Bac N, Lahmer T, Zhuang X, Nguyen-Thoi T, Rabczuk T. A software framework for probabilistic sensitivity analysis for computationally expensive models. Advances in Engineering Software, 2016, 100: 19–31

[18] RuppertDWandM PCarrollR J. Semiparametric Regression (No. 12). Cambridge: Cambridge University Press, 2003

[19] Vu-Bac N, Silani M, Lahmer T, Zhuang X, Rabczuk T. A unified framework for stochastic predictions of mechanical properties of polymeric nanocomposites. Computational Materials Science, 2015, 96: 520–535

[20] Papaioannou I, Straub D. Variance-based reliability sensitivity analysis and the FORMα. Reliability Engineering and System Safety, 2021, 210: 107496

[21] Kucherenko S, Tarantola S, Annoni P. Estimation of global sensitivity indices for models with dependent variables. Computer Physics Communications, 2012, 183(4): 937–946

[22] Vu-Bac N, Rafiee R, Zhuang X, Lahmer T, Rabczuk T. Uncertainty quantification for multiscale modeling of polymer nanocomposites with correlated parameters. Composites. Part B, Engineering, 2015, 68: 446–464

[23] ChenTGuestrinC. XGBoost: A scalable tree boosting system. In Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining. New York, NY: Association for Computing Machinery, 2016, 785–794

[24] CholletF. Deep Learning with Python. Shelter Island, NY: Manning Publications, 2018

[25] Vu-Bac N, Lahmer T, Keitel H, Zhao J, Zhuang X, Rabczuk T. Stochastic predictions of bulk properties of amorphous polyethylene based on molecular dynamics simulations. Mechanics of Materials, 2014, 68: 70–84

[26] Vu-Bac N, Lahmer T, Zhang Y, Zhuang X, Rabczuk T. Stochastic predictions of interfacial characteristic of polymeric nanocomposites (PNCs). Composites. Part B, Engineering, 2014, 59: 80–95

[27] ForresterASobesterAKeaneA. Engineering Design via Surrogate Modelling: A Practical Guide. Hoboken, NJ: John Wiley & Sons Ltd, 2008

[28] Americancoal ash association. Fly Ash Facts for Highway Engineers. Technical report FHWA-IF-03-019. 2003

RIGHTS & PERMISSIONS

Higher Education Press