Reliability assessment of compressive and splitting tensile strength prediction of roller compacted concrete pavement: introducing MARS-GOA-MCS

ABSTRACT It is required to determine the mechanical properties of Roller Compacted Concrete Pavement (RCCP) such as Compressive Strength (CS) and Splitting Tensile Strength (STS) for constructing and maintaining the road pavement because of the inevitable uncertainties in the experimental setup. In this study, a reliable hybrid intelligent model based on the integration of Grasshopper Optimisation Algorithm (GOA) and Multivariate Adaptive Regression Splines (MARS), called MARS-GOA, was developed to present new closed-form equations for CS and STS prediction of RCCP. Initially, GOA was developed to serve as a search engine of the proposed algorithm to optimise MARS control parameters by minimising the error of prediction. Statistical metrics indicate that proposed hybrid MARS-GOA outperformed ELM, M5p, and standard MARS for prediction of both the CS (CoD = 0.811, PMARE = 22.146% and U95 = 34.670) and STS (CoD = 0.816, PMARE = 16.192% and U95 = 2.725) of RCCP. The results of reliability analysis, namely Monte-Carlo Simulation (MCS), revealed that TCS, TTS, and uncertainty level of the used dataset play a major role in the determination of the reliability of the RCCP mix design.


Introduction
In the last decades, the highway sustainability, pavements maintenance and construction have been regarded as a main challenge in road and pavement engineering. Roller Compacted Concrete Pavement (RCCP) is an economical, ecoenvironmentally and short-run construction alternative for pavement challenges because of using the minimum cement volume and maximum aggregate content, which have been led to its increased consistency and durability compared to bituminous asphalt (Hashemi et al. 2018, Modarres et al. 2018, Lam et al. 2018, NoParast et al. 2021). In addition, RCCP is characterised by low water absorption, high-temperature stability, satisfactory mechanical properties and low deformation under loads. RCCP has been mainly utilised for pavements with heavy duty including long handling courtyards, local streets and airport terminals among others (Chhorn et al. 2019). Also, it is considered as an eco-friendly material due to its impermeability characteristic. It is noteworthy that using the Supplementary Cementitious Materials (SCMs) for ensuring the appropriate compaction in the blends with natural sand for producing RCCP has recently attracted much due to improved mechanical properties with the lower cement cost production (Adamu et al. 2020). Thus, this study aims to investigate the RCCP blends containing pozzolanic SCMs. SCMs when used in conjunction with Portland cement, Portland limestone cement or blended cements, contribute to the properties of hardened concrete through hydraulic and/or pozzolanic activity.
The parameters of the experimental programme including proportions of mixture, cast, curing condition, as well as the age of specimens (AS), may affect the behaviour of concrete strength (Ashrafian et al. 2018). During the mixture designing, it was sought to produce reliable, cost-effective concrete. Of course, this is a time-consuming process, which requires a complicated controlling of the conditions and assumptions underlying the concrete ingredients (Taheri Amiri et al. 2019, Ashrafian et al. 2019, Shahmansouri et al. 2020). Therefore, this study proposes a number of linear and non-linear regression (LR)-based predictive models relying on a variety of statistical methods for the prediction of the concrete behaviour.
In the last decade, Data-Driven Models (DDMs) have been developed and employed in concrete industry in order to predict its properties such as mechanical and structural ones as well as its durability, (Gholampour et al. 2017, Amlashi et al. 2019, Asteris et al. 2019a, Ly et al. 2019, Sun et al. 2019, Ashrafian et al. 2020b, Liu et al. 2020. For example, Deng et al. (2018) analzyed an Artificial Neural Network (ANN) for the estimation of the recycled aggregate concrete's strength capacity. Moreover, Ashrafian et al. (2020aAshrafian et al. ( , 2020bAshrafian et al. ( , 2020c presented a metaheuristic regression method for simulation of the strength characteristics in lightweight concrete with foam. Feng et al. (2020) used an adaptive boosting approach based on the machine learning to present a regression model for concrete's CS. Sun et al. (2019) developed a Support Vector Machine (SVM) to simulate the durability and hardness. Iqbal et al. (2020) used a Gene Expression Programming (GEP) to propose the empirical models for predicting the mechanical properties of environmentally friendly concrete incorporating waste foundry sand. Shahmansouri et al. (2019) developed a linear GEP model to deal with the role of the mixture designing of ecologically sound concrete with zeolite and CS. Han et al. (2019) investigated the Random Forest (RF) technique for estimation of the strength properties of concrete with high-performance. Asteris et al. (2019b) analyzed the M5prime model tree (M5p) and Multivariate Adaptive Regression Spline (MARS) to formulate the self-compacting characteristic of the concrete. Gholampour et al. (2018) employed a tree model to estimate the recycled concrete's CS. Ashrafian et al. (2020aAshrafian et al. ( , 2020bAshrafian et al. ( , 2020c proposed a classified method to reveal the RCCP's mechanical properties. Zhang et al. (2021) proposed an auto-estimated model for lightweight foamed concrete's mixture design based on the multi-objective firefly algorithm and SVM method. Ahmad et al. (2021) developed an evolutionary artificialintelligence-based approach for predicting the hemp biocomposites' performance via GEP. Shahmansouri et al. (2021) analyzed the geopolymer concrete's mechanical properties containing the natural zeolite as well as silica fume using an optimum design based on the response surface method. Ren et al. (2021) sought to predict the concrete's CS using the manufactured sand and ensemble classification as well as regression tree method.
The CS and Splitting Tensile Strength (STS) of RCCP are considered the non-linear characteristics of its ingredients because they are realised only by its discrete results. In this case, selecting the most suitable forms, the model's degree, and appropriate values of the model parameters considerably affect the efficiency and accuracy of the developed prediction model, thereby providing an effective representation of different RCCP mix parameters. Although classical DDMs provide estimation models that may have accuracy and robustness for the prediction of the mechanical characteristics of concrete, they are complicated and may have a high-cost training process. MARS is a non-parametric regression technique that can automatically model non-linearities and interactions between variables, even if no prior information on the form of the response function and/or no explicit relation between independent and dependent variables is available (Zhang and Goh 2016). Nevertheless, similar to ANNs (Momeni et al. 2015), MARS implementation may have some disadvantages, i.e. slow rate of training and getting trapped in local minima. Besides, effective tuning of the model input parameters such as the maximum number of basis functions (M max ), maximum number of interactions (I max ), and the penalty parameter (c) is also a challenging task that can affect the efficiency and reliability of the MARS-based developed models. To overcome these disadvantages, the use of powerful optimisation algorithms to enhance MARS performance is of great interest. Integrated models based on the combination of DDMs and global optimisation algorithms have an attractive viewpoint of limited control factors for the implementation of novel methods and reduction of the model's overfitting (Nieto et al. 2015).
In the present study, the Grasshopper Optimisation Algorithm (GOA) as an economical simple and high convergence stochastic search algorithm with flexible computational performance (Saremi et al. 2020) are firstly employed in this research for tuning of MARS control parameters. The GOA algorithm increases the average fitness of solutions, thereby effectively enhancing the initial random population. The target fitness is improved during the iterations, thereby leading to the more accurate estimation of the global optimum proportional to the number of iterations. A high number of studies have discussed the capability of GOA in comprising other optimisation algorithms (Mirjalili et al. 2018, Barman et al. 2018. The main contributions of the present study are as follows: (1) Development of the MARS-GOA-based hybrid intelligent model for simulation of the CS and STS of RCCP. Furthermore, the emerging closed-from predictive equations for CS and STS of RCCP are also presented using standard MARS and the proposed hybrid MARS-GOA.
(2) Comparison of the accuracy and efficiency of the proposed hybrid MARS-GOA optimised models with the other employed DDMs including the standard MARS, regression technique based classification called M5p, methods based on committee called Extreme Learning Machine (ELM) and ANN, as well as ensemble method named RF for CS and STS prediction of RCCP. (3) A reliable optimum mix design in RCCP is important because the trivial change in the content of materials can drastically change the freshness and hardness (Ashrafian et al. 2020a(Ashrafian et al. , 2020b(Ashrafian et al. , 2020c. In commonly developed models for predicting mechanical and structural characteristics of concrete mixtures, the regression-based predicted responses are calculated with regard to deterministic input parameters, while it is well known that uncertainties in structural design and construction processes are often inevitable. Due to the no consideration for uncertainties, the final deterministic designs may be unreliable and most of the time result in failure (Aggarwal et al. 2015, AzariJafari et al. 2019). Therefore, this study also aims to analyze the effect of uncertainties of input variables on the reliability level of the developed models for RCCP concrete mix design. To this end, by considering the influential input parameters as random variables, the effects of uncertain parameters on the reliability analysis results are investigated for different values of target CS and target STS by the Monte-Carlo Simulation (MCS) method.
This paper is organised as follows. The statistical parameters of the relevant experimental dataset, the employed DDMs including M5p, ELM, ANN, RF, MARS, and details of the proposed hybrid MARS-GOA model are presented in Section 2. The computational process and comparison of proposed developed models in conjunction with reliability analysis results for CS and STS prediction of RCCP are presented in Section 3. Finally, concluding remarks and limitations are presented in Section 4.

Theoretical framework and data description
The cost-effective design of the environmentally friendly high-quality concrete mixture is a challenging task (AzariJafari et al. 2019). Additionally, it is required to consider a number of experimental variables in constructing the economical RCCP mixtures with mechanical properties in a short time . RCCP has become very convent due to its simple and fast production while producing a great surface.
The strength of the pavement design and its application affect the selection of materials used in RCCP mixtures. The basic materials for RCCP production are cement, SCMs, water and aggregates. Generally, the cost of materials for RCCP is almost equal to the cost of materials used in ordinary Portland concrete. Nevertheless, using a number of economical ways is possible due to the lower cement contents usually required in RCCP mixtures to attain strengths equivalent to those of ordinary concrete (ACI 325-10R-95 2001).
In this study, for the development of DDMs-based reliable and adaptable estimation concepts, data samples at different ages of specimens were collected using the peer-reviewed experimental programme document (Atiş et al. 2004, Tangtermsirikul et al. 2004, Mardani-Aghabaglou and Ramyar 2013, Pavan and Rao 2014, Rao et al. 2015a, 2015b, 2016a, 2016b, Hesami et al. 2016, Aghaeipour and Madhkhan 2017, Ghahari et al. 2017, Shamsaei et al. 2017, Vahidi et al. 2017, Hashemi et al. 2018, Debbarma et al. 2019. The relevant data contained the main influential parameters about mixture proportions of RCCP in various constituents. Based on this collected data, modelling of the CS and STS of RCCP was carried out using 621 and 326 experimental data samples, respectively. The final dataset was divided into two subsets to develop the DDMs. The train (calibration) data was used for constructing the models for the CS and STS prediction. To evaluate the performance of recruited models, the test (validation) data were also employed. In this regard, about 75%, 466 for CS and 245 for STS, and 25%, 155 for CS and 81 for STS, were considered as the calibration and validation groups, respectively.
The statistical analysis of the most effective variables for designing the concrete with proper mechanical properties such as the coarse aggregate (CA), fine aggregate (FA), SCMs, water (W), binder (B) and AS, are shown in Table 1. Additionally, the r-Pearson correlation analysis of the relevant input parameters for developing the proposed DDMs has been illustrated in the heatmap as displayed in Figure 1, where there are no significant correlations between the selected matrixes versus CS and STS. In 1992, the M5p model was formed on the basis of decision trees and LR approaches. The primary terminal node and some more lead nodes comprise a binary decision tree. This creates a connection between the inputs and target instances, i.e. independent or dependent parameters (Quinlan 1992). Note that, although decision trees are mostly used for categorical type data, the type of data suitable for them is the quantitative type (Rezaie-Balf et al. 2019a). Figure 2 shows the M5p model. It has two major steps as follows: (1) Generate a decision tree by separating the input data, which can be done while resolving the standard deviation of every subset for reaching a reasonable primary node or parent node. Upon the splitting step, the child node would have a smaller std compared with the parent node. (2) Solving the error by having every node in the decision tree tested. Standard deviation for quantifying the error values can be measured using the following equation: where T j corresponds to the pattern subset with the jth output of the potential set and T stands for a set of instances with the ability to reach the primary node.
The standard deviation process of M5 captures a node. All subspace attributes are evaluated to derive an std; the computed std then is compared to the initial std and the expected error reduction will be computed, and then, the attribute with the maximum SDR is selected. When the primary value percentage of std becomes more than SDR, the splitting process stops. After that, LR relationship map is built. So, the most potential error-reducing node is selected through different splitting processes. A pruning technique, based on LR functions, is also used to eliminate sub-trees to avoid over-fittings. This model learns better and solves high-complex problems better than the M5rule model (Rezaie-Balf et al. 2019). Huang et al. (2006) proposed a new type of the training feedforward neural network (FFNN), named ELM, which is like ANN and selects the weights of hidden neurons based on FFNN. A Single-Layer FN (SLFN) with N hidden nodes, h(x) as an activation function, and M samples can be written as follows:

Extreme learning machine
where W i and b i denote the input weights and biases in L hidden nodes, respectively. ELM determines the input weights (W i ) and biases (b i ) using the least-square (LS) technique to calculate the output weights β. In addition, activation function [h(x)] can be expressed as follows: where H is the hidden layer output matrix which is defined as follows: The ELM mainly aims to minimise ||b|| that is equal to the 2 ||b|| maximisation (Rezaie-Balf et al. 2019a). ELM with a simple algorithm can reduce the training time of the model construction. In this model, the hidden nodes' parameters are randomly selected and the output weights are calculated by the LS method (Rezaie-Balf and Kisi 2018). The one-layer architecture ELM network is shown in Figure 3. More details on ELM are available in Rezaie-Balf and Kisi (2018) and Miche et al. (2008).

Artificial neural network
The ANN similar to real nervous systems is considered a high number of processing elements or neurons and the layers   which are connected by weights, parallelly (McCulloch and Pitts 1943). A three-layered ANN with layers i, j and k as well as weights W ij and W jk is shown in Figure 4.
The input x is determined by the weighted output sum of the first layer and then allocated to each neuron in the second and third layers. For instance, the y in the second layer j is computed as follows (Haykin 1999): where u j , O pi and W ij denote on the bias for neuron j, the ith denotes on the output of the first layer and the weights between first and second layers, respectively. Logistic function as a popular activation function of ANN can be expressed in the form of following equation:

Random forest
The proposed RF by Breiman is a non-parametric and classification-based regression model (Breiman 2001). Random Forests incorporate a high number of easy-to-interpret decision trees rather than the parametric models. By integrating the results of the decision tree models, a more comprehensive prediction model is obtained. This study mainly aims to estimate the RCCP mechanical properties through the regression models' capability. The random forests regression model as a non-parametric regression model consists of a set of K trees {T 1(X) , T 2(X) , … , T K(X) } where X = {x 1 , x 2 , … , x β } is the input vector of β-dimension which forms a forest. The ensemble produces P outputs corresponding to each tree Y p (p = 1, 2, … , P).
The final result is obtained by averaging all the tree predictions.

Proposed hybrid intelligent MARS-GOA algorithm
In this section, at first the MARS and GOA are briefly described. Then, details of the proposed hybrid MARS-GOA are presented in the next subsection.  2.2.5.1. Multivariate adaptive regression splines (MARS). One of the non-parametric and regression-intelligence models that have been introduced by Friedman (1991) is MARS according which the independent and dependent variables of a knowledge system can be modelled using splines. In MARS, the Basis Function (BF) based on the slope of the regression line derived in the final model (one spline to more) has been altered. MARS is considered flexible in terms of evaluating the non-linear procedures in predictor and target variables, compared to other widely used methods in regression-based data intelligent approaches (Hastie et al. 2008, Jekabsons 2010). An overview of the MARS paradigm is shown in Figure 5.
After acquiring the MARS model, the relationship is obtained by combining all BFs such as intercept BF and pairwise BF with interaction.
Two types of BFs could be used for mapping using input parameters (X p ) to response Y.
where t presents the threshold value. MARS could search probable situations in all degrees and reveal a solution for variable interactions because of the extraction of the complicated data structures from multi-dimensional datasets. A further parameter is investigated on the forward step by defining a maximum allowable degree of interaction (I max ). Typically, the first and second orders of interaction are just allowed, however, higher degrees are used if the data warrants it. The general equation for the MARS function is expressed as follows (Rezaie-Balf et al. 2019a, 2019b: Here f (x) is the predicted target, which corresponds to the predictor variable. Also, x, b 0 and b m are the predicted constant coefficients so as to obtain the most appropriate data fit, and M denotes the number of BF. The constant coefficients b and l in Equation (10) can be calculated using the LS method. In MARS, the locations of the knots are determined by an adaptive regression algorithm and BFs are obtained by a stepwise search (Hastie et al. 2008). The forward and backward stages are two steps for constructing MARS model. According to Equation (11), MARS adopts the Generalised Cross-Validation (GCV) algorithm to remove generated BFs in forward stage, thereby preventing over-fitting in the backward stage (Rezaie-Balf et al. 2019a, 2019b, Pham et al. 2020. The GCV formula is as follows: where N, M, and c denote the number of training data records, number of BFs, and the penalty for each BFs, respectively. 2.2.5.2. Grasshopper optimisation algorithm. The GOA which was introduced by Saremi et al. (2020) simulates grasshopper swarms' behaviour in nature and is used in different engineering fields. The swarming behaviour of grasshoppers is found in both nymphs (with no wing) and adults (with wings). The adults are considered for searching the entire search space globally (exploration) and finding more appropriate food source regions whereas the nymphs are employed for exploiting a specific region or neighbourhood of a particular position (exploitation) (Mirjalili et al. 2018). Here, the algorithm is clarified as follows (Barman et al. 2018): Step 1: At first, a population of size S c for GOA is prepared by the following equation: Here S c denotes population size and N is the dimension of the problem. u j and l j are the upper and lower bound of the jth variable.
Step 2: The best place in this step is calculated based on the fitness as follows: where the f(x j ) represents the objective function or the Mean Square Error (MSE) in the present study.
Step 3: For balancing the exploration and exploitation, GOA is equipped with 'c' which is varied based on the iteration where M I denotes the maximum cycle number.
Step 4: Then, a new position is determined using the following equation: where s(r) = f .e (−r/l) − e −r (16) T d = best-found solution so far, f and l denote the attraction intensity and the length of attraction. The new position is evaluated by Equation (15) and the better position is preserved.
Step 5: Finally, an iteration counter (iter) increases and the procedure is repeated from Step 4.
Step 6: This process continues until the maximum cycle number (M I ) is obtained for iter.
For more details on GOA structure and its optimisation, see Saremi et al. (2020).
2.2.5.3. Details of the proposed hybrid MARC-GOA prediction model. In the DDMs, the selection of the hyperparameters and optimal platforms of the employed prediction model is a challenging task and a crucial process that can affect the accuracy and efficiency of the prediction model. This issue is manifested based on the complexity of the concept of the modelling and database comprehensiveness. One of the solutions for the detection of the optimum DDM parameters is reviewing the relevant literature by checking the optimistic solutions based on the restricted heuristic search. However, this kind of trial and error is inaccurate and might fail and return an unsatisfied solution due to a restricted heuristic search.
MARS technique was developed based on the key factors such as the maximum basis functions (M max ), the penalty factor (c), and the order of interaction (I max ). Searching for optimum key factors is a very difficult challenge, but detection of the proper parameters may considerably improve the MARS performance in modelling the phenomena. Therefore, an intelligent non-parametric model based on the hybrid GOA -MARS is proposed to help scholars to overcome this challenge ( Figure 5). There exists a higher convergence probability at a local minimum by original MARS, while GOA could find a global optimum. Hence, the GOA-based MARS model proposed in the present study benefits from the capabilities of both GOA and MARS in which GOA looks for the optimal key factors (c, M max and I max ) in the search space; in addition, MARS uses the optimised parameters to find the satisfactory results. Figure 6 shows the general flowchart of the proposed hybrid MARS-GOA.
As shown in Figure 6, the training procedure in the GOAbased MARS model begins with initialisation of a grasshopper population. In this step, the grasshopper positions that reveal the key factors of MARS are randomly assigned in the search space according to Equation (10). Subsequently, the GOAbased MARS model is trained using the initial position of grasshoppers (c, M max , and I max ). Then, errors, which occur between predicted and actual values are calculated. These errors at each iteration are reduced by updating the grasshopper position. This procedure is repeated until the termination criteria are met.
During the optimisation, MARS is performed to construct a new model for each set of key factors which are prepared by each candidate solution of GOA. The capability of models due to fitness function is compared using the greedy selector in GOA. When the calibrating (train) subset is carried out, MARS is proposed for the validation (test) dataset. For achieving the highest performance, minimisation of prediction error is considered as an objective function (f ): where E train and E test denote the error of 75-25% data dividing feature, which is considered for calibrating (train) and validating (test) phases, respectively.
In order to minimise MSE as the error measurement, candidate solutions move within the search space in each iteration; subsequently, the MARS model is retrained and the new error is calculated using the new key factors till satisfactory error is attained. It is noteworthy that the maximum generation number is considered as the termination condition of the GOA search process. Finally, the obtained optimum parameters are used for the proposed prediction model.

Statistical performance for evaluation
The following common performance indices are used to compare the hybrid MARS-GOA with standalone ELM, ANN, MARS, M5p and RF models in prediction of CS and STS properties of RCCP.
where t exp and t pre denote the experimental and predicted target variable values (e.g. CS and STS), respectively. t exp and t pre are the average values of the experimental and estimated target, respectively. Also, N indicates the total amount of data. The R 2 metric which is in the range of (0-1) (with R 2 = 1 as the best value) determines the adequacy of selected input variables in the estimation of the target. RMSE with the range of zero to positive infinity and an optimum value of zero is utilised to evaluate the performance. OI with the range of (0, +1) and optimistic value equal to zero is utilised for evaluation of the sufficiency and goodness of fit of the developed models. PMARE (0, +1) on the basis of relative error is measured to evaluate the descriptive and logical behaviour of the developed model, and U95 is also evaluated as a 95% uncertainty band for verification of model validity.

Development of the standalone models for CS and STS of RCCP
In order to implement M5p model for CS and STS of RCCP, a series of tuning parameters should be set for the purpose of classification-based DDMs. For extracting LR relationships of CS and STS, pruning parameter and smoothing condition were considered for evaluation of M5p prediction capacity. With 6 independent variables and 2 dependent variables, 12 and 18 linear models were extracted from the M5p model for the prediction of target variables. The optimistic number of trees in M5p models has been set when this parameter can provide minimum error prediction in the training and constructing model according to Figure 7. Moreover, as in the M5p rule, CS and STS predictive DDMs are proposed for AS both lower and higher than 10.5 and 21, respectively. Coefficients of M5p models for linear models CS and STS characteristics are listed in Tables A1 and A2 in the appendix. It is noteworthy that all of the input variables are considered in predicting the strength characteristics of RCCP; they are implemented by M5p's linear models significantly. In the case of the ELM predictive model, a neural platform containing one-layer network was carried out. The most widely used transfer function consisting of the sigmoid formula was applied by reviewing the relevant literature for feature extraction (Rezaie-Balf et al. 2019a, 2019b, Pham et al. 2020. As shown in Table 2, for identifying the optimum  ELM architecture, 10-120 hidden neuron nodes were analyzed and each architecture was evaluated using a testing group (25% of the dataset). For initialisation of the random values of the hidden layer in ELM, the model was run 1000 times and the minimum MSE was selected for the best structure.
The default of the Bagger algorithm for the implementation of RF was used in a bag size percent tuned at 150, delta criterion set was 0.1003 and the leaf number was 8. It should be stated that no general mathematical formulation is used to find the optimum tree number. The ANN method can use the model construction of a neural network for perceiving the complex non-linear mathematical model, notably when there is no mathematical expression between independent variables and target. In this study, to estimate CS and STS of RCCP, the models were proposed including six input layers with six neurons. The optimal neurons selection using this structure for ANN topology trained by the Levenberg Marquardt (LM) algorithm was presented in Table 3.
In the present study, Jekabsons code was used to construct a non-parametric regression-based model for developing the MARS model. The best MARS model had 26 and 35 basis functions for CS and STS characteristics of RCCP, whereas the I max level was determined in second and third-order interaction, respectively. Also, smooth parameters were investigated using trial and error in the range of 1-8 according to Jekabsons (2010), and the best value was obtained as c best = 1 for CS and c best = 2 for STS parameters based on the minimum MSE values. The total number of effective parameters was 27 and GCV = 58.620 for CS, and 71 and GCV = 0.4689 for STS predictions. Table B1 in the appendix section represents details of BFs of the MARS model. The interpretable model based on the non-parametric formula was defined as follows: where C i and BF i were referred as coefficient and basis function of MARS procedure.

Proposed hybrid MARC-GOA models for CS and STS of RCCP
The parameters of GOA (e. g., C min , C max , and f ) should be tuned up before optimisation of MARS control parameters. In this study, the values of C min = 0.1, C max = 0.9 and f = 0.5 are considered for regularising the capability of GOA. Figure  8 shows the MSE values as an objective function versus the number of iterations for the proposed GOA-based MARS model. The high convergence speed of the algorithm for the prediction of both CS and STS of RCCP is quite obvious. As can be seen from Figure 8, the employed algorithm provides an efficient and reliable estimation for both CS and STS of RCCP after approximately 60 iterations. It should be mentioned that the volumetric-weighted forms of predictor variables of RCCP mixtures were utilised in the proposed MARS-GOA as it was found that this substantially enhances the model's training performance. The predictive formulation analysis using MARS-GOA adopted 47 and 41 BFs of piecewise-linear spline functions with second and third-order interaction for CS and STS, respectively. Also, the best values for the MARS penalty factor were obtained by the GOA as c best = 0 and c best = 1 for CS and STS, respectively.
The BFs-based equations of MARS-GOA for the formulation of CS and STS are presented in Table B2 in the appendix. It can be concluded from Table B2 that interaction term has occurred between BFs (34 of the 46 BFs for CS and 30 of the 39 BFs for STS). The existence of this condition identified that the developed MARS-GOA model is not simply additive and that interactions play a crucial performance in the construction of a reliable and accurate predictive model for CS and STS prediction of RCCP. This characteristic defined that the non-linear MARS-GOA is capable to capture independent variables and the target relationships without making any particular assumption about the underlying functional contribution. Furthermore, it can be shown that there are six knots for CA, at values of 1012, 1209, 1194, 1095, 1097.9 and 962.5 kg/m 3 for CS while there are three knots for binder, at values of 295, 372.45 and 400 kg/m 3 for STS. These knot values clarified researchers to have an insight into the perspective of where significant changes in the experimental dataset may occur. Finally, the interpretable MARS-GOA predictive equations for CS and STS of RCCP are developed as

Results comparison of proposed hybrid and standalone models
In this section, the accuracy and applicability of the employed DDMs called MARS, M5p, ELM and the proposed hybrid MARS-GOA for the prediction of the STS and CS of RCCP are compared. The statistical measures for the assessment of the proposed DDMs for predicting the CS and STS are reported in Table 4 and Table 5, respectively. As shown in Table 3, the proposed hybrid MARS-GOA model (R 2 =  The performance metrics of presented DDMs for STS prediction in train and test phases revealed that MARS-GOA proposed accurate model in terms of R 2 = 0.869, OI=1.028 MPa for training and R 2 = 0.816, OI=1.046 MPa for testing compared to other developed methods. For the testing step, MARS-GOA has the minimum RMSE (7.077 and 0.555 MPa) and U95 (2.662 and 2.725) for compressive and tensile strength, respectively; it improved the accuracy of the testing subset for prediction of CS and STS in terms of PMARE for the ELM (5.219% and 2.075%), M5p (3.021% and 2.469%) and MARS (5.219% and 2.075%), respectively. Figures 9-12 showed the plots for comparing the experimental predictions using DDMs. The results of the proposed models for CS prediction of RCCP over the training and testing data sets are reported in Figures 9 and 10, respectively. Figures 11  and 12 show the results for STS of RCCP over the training and testing data sets, respectively.
The closer unity (red and dotted line X = Y ) in Figures 9-12 reveals that the proposed hybrid MARS-GOA had a better visual agreement of measured and predicted values of CS and STS characteristics. It is noteworthy that the experimental values of CS and STS parameters illustrated robust coherence with the experimental data points and the most accurate results of CS and STS are provided by proposed MARS-GOA for both the training and testing subsets. Furthermore, most deviations from ±30% black lines (lowest accuracy) corresponded to M5p and ELM for the CS and STS prediction of RCCP, respectively.
In order to present more statistical details on the performance of the developed intelligent GOA-based hybrid and independent predictive models, the Taylor diagrams for CS and STS prediction of RCCP are depicted in Figure 13. According to the predicted results, all of the implemented DDMs may be displayed in terms of their close relationship and in terms of prediction of the experimental STS and CS.
As shown in Figure 13, the statistical performance measurements including the standard deviations and correlation coefficients of the presented DDMs are used to determine the comprising between the experimental observations and predicted values. According to graphical visualisation, the proposed hybrid MARS-GOA intelligent model used quite the most accurate performances compared to the ANN, ELM, RF, M5p, and the standard MARS model. Eventually, it is evident that MARS models integrating GOA could increase the accuracy of the prediction values for both properties of RCCP compared to other DDMs.

Reliability analysis using MCS
In this section, reliability analysis for CS and STS prediction of RCCP is performed using standard MARS and the proposed hybrid MARS-GOA models to deal with uncertainties involved in the concrete mix design process (BIS 1959). For this purpose, all input variables (CA, FA, SCM, W, B and AS) of the RCCP concrete mix are considered as random variables, X = [CA, FA, SCM, W, B, AS]. The RCCP mix design prediction reliability problem considering probabilistic concrete mix input parameters showing random characteristics is then formulated as: where X is a vector of random variables, f X denotes the joint probability distribution function (PDF) of the random vector, y denotes the calculated response for the CS or STS of RCCP obtained by the employed prediction model, and y t is the desired value for the target compressive strength (TCS) or target tensile strength (TTS) of RCCP concrete mix design. The MCS, as the most accurate reliability analysis method (Cardoso et al. 2008, Safaeian Hamzehkolaei et al. 2016, Safaeian Hamzehkolaei et al. 2018, Pham et al. 2020, is employed for evaluating the reliability of RCCP concrete mix design. The key step in MCS is generating random samples for uncertain variables based on the respective PDFs and, then, simulating the performance of the system for the generated samples. Using MCS, the probability integral in Equation (27) can be calculated as  where E f is the expectation operator, {x i } N i=1 is the i-th generated sample, and I[ ] is the indicator function defined as In order to reliability analysis of RCCP concrete mix design, a normal PDF is considered for all random variables. The mean values of random variables (m X ) are fixed at their respective mean values of the dataset listed in Table 1 ( . In this case, simulated data are generated around the average values of variables are decreased by increasing the standard deviation. For investigating the effect of uncertainty levels of the input data, standard deviations of random variables (s X ) are gradually increased within the interval of [0.1-1] times of the measured standard deviations of the used dataset (s d ).   Equations (25) and (26) for the standard MARS and the hybrid MARS-GOA models, respectively.
As can be seen from Figures 14 and 15, for all cases under consideration, the reliability of CS and STS prediction for both the standard MARS and the proposed hybrid MARS-GOA models are greater than 50%. The results show that the target CS and target TS have a significant influence on the reliability levels of the developed models. It is clear that higher reliability levels for RCCP mix design can be obtained when the TCS or TTS is decreased. Besides, for all uncertainty levels (a) of the random variables, the reliability levels of the hybrid MARS-GOA for both the CS and STS prediction of RCCP are much greater than that of standard MARS models. It is worth mentioning that the reliability of prediction for both the CS and STS models of RCCP could be hugely enhanced by decreasing the standard deviation scale factor. As can be seen in Figures  14 and 15, the reliability levels of the proposed MARS-GOA models are higher than 85% for a , 0.5. Therefore, it can be concluded that, more robust and more reliable prediction models can also be obtained when the uncertainty level of the used dataset is decreased. Results of this section verified that the proposed MARS-GOA enhanced both the accuracy and reliability levels of the developed models compared to standard MARS.

Study limitation and future research
The findings of this study identified the superiority of the MARS-GOA in constructing robust and efficient regressionbased DDMs for CS and STS formulation of RCCP. It was shown that the proposed hybrid intelligent algorithm with the optimised MARS control parameters not only enhanced the accuracy but also increased the robustness and reliability of the CS and STS prediction of RCCP compared to the standard MARS model. Therefore, it can be employed to develop reliable relationships that can serve as good substitutes for empirical models to design material properties of RCCP. However, the time-consuming process of the hyperparameters tuning of the MARS may be the main challenge of the developed algorithm. Coupling efficient DDMs and data mining approaches such as the feature selection of inputs, data preprocessing, etc., in conjunction with advanced global optimisation algorithms with high convergence speed could be addressed in a future study to cover modelling limitation.   Figure 15. Reliability sensitivity results for STS prediction of RCCP.

Conclusion
RCCP is an eco-friendly and cost-effective type of concrete for road and highway applications. This study focused on modelling and formulation of the CS and STS for construction and maintenance operation of RCCP using a novel MARS integrated GOA. In the proposed hybrid models, MARS control parameters were tuned using GOA by minimising the error of the prediction. The efficiency and capability of the developed models with the optimised input control parameters were compared with different benchmark DDMs including standard MARS, M5p, ANN, RF and ELM, by utilising several evaluation measures to assess the predictability of the predictive DDMs. Based on the modelling process, error assessment, and visual interpretation such as scatter plots, radar plots, and Taylor diagrams, MARS-GOA models improved the precision of the CS and STS prediction through hybridising of GOA.
The reliability sensitivity analysis of the developed models for RCCP concrete mix design was also performed using MCS by considering the most influential input parameters as random variables with normal PDFs. Reliability analysis results for a wide range of uncertainty levels verified that the target CS and target STS have a significant influence on the reliability of the RCCP concrete mix design. For all values of the target CS and target STS, the proposed hybrid MARS-GOA models outperformed standard MARS models. Besides, the reliability of the developed modes was highly increased by decreasing the uncertainty level of the used dataset.   (