Application of Machine Learning Techniques for Predicting Compressive, Splitting Tensile, and Flexural Strengths of Concrete with Metakaolin

The mechanical properties of concrete are the important parameters in a design code. The amount of laboratory trial batches and experiments required to produce useful design data can be decreased by using robust prediction models for the mechanical properties of concrete, which can save time and money. Portland cement is frequently substituted with metakaolin (MK) because of its technical and environmental advantages. In this study, three mechanical properties of concrete with MK, i.e., compressive strength (f′c), splitting tensile strength (fst), and flexural strength (FS) were modelled by using four machine learning (ML) techniques: gene expression programming (GEP), artificial neural network (ANN), M5P model tree algorithm, and random forest (RF). For this purpose, a comprehensive database containing detail of concrete mixture proportions and values of f′c, fst, and FS at different ages was gathered from peer-reviewed published documents. Various statistical metrics were used to compare the predictive and generalization capability of the ML techniques. The comparative study of ML techniques revealed that RF has better predictive and generalization capability as compared with GEP, ANN, and M5P model tree algorithm. Moreover, the sensitivity and parametric analysis (PA) was carried out. The PA showed that the most suitable proportions of MK as partial cement replacement were 10% for FS and 15% for both f′c and fst.


Introduction
Concrete is the second most-consumed substance in the world after water. Portland cement clinker is the base of the majority of cementitious binders that are being used in concrete. The production of cement consumes a huge amount of energy and releases about 7% of CO 2 into the atmosphere [1]. During cement production, the calcination process, in which CaO is formed by extracting CO 2 from CaCO 3 , is responsible for about 50% of CO 2 production while the remaining 50% CO 2 is produced by energy use [2]. The input component. Age improves hydration products. By increasing time, the hydration of anhydrous cement particles occurs, which increases the amount of C-S-H gel and increases strength [32]. This paper has two aims: (1) to model and compare the f c , f st , and FS of concrete incorporating MK by an evolutionary algorithm (GEP), ANN, and decision trees (M5P and RF). (2) To explore the influence of different input parameters on the mechanical properties of concrete with MK by using parametric analysis. Moreover, sensitivity analysis was also carried out in order to find out the relative contribution of input parameters on mechanical properties. It is important to mention that the compressive strength of concrete with MK was not modelled by using M5P because of the large database (around 982 data points) which generates a significant number of linear models (greater than 40). More suitable ML techniques for the prediction of mechanical properties of concrete with MK may be found.

Data Collection
The database was collected from the literature and outliers were deleted. The remaining data comprised 982 data points for f c , 204 data points for f st , and 63 data points for FS. Our aim was to collect a large database and, therefore, all test results of mechanical properties of concrete with MK conducted based on international standards were collected for use in this study. This database was collected from 45 published studies [4,[15][16][17][18] and is shown in the supplementary document. In the supplementary document, Tables S1-S3 contain the information of input parameters and output for f c , f st , and FS database, respectively. The collected database includes C, MK, w/b ratio, FA, CA, SP, and days as input parameters and f c , f st , and FS as output parameters. The descriptive statistics of input and output variables used in the training set are shown in Table 1. These values help to give insights into the range and distribution-related properties of independent and dependent variables used in the training set.
All the data points used for modelling of f c of concrete with MK were for cube 150 mm. The f c of cylinder ∅100 mm × 200 mm was converted to f c of cube 150 mm by multiplying with factor 1.1 for normal-strength concrete (NSC) and 0.98 for HSC (according to ACI 363R, HSC is concrete that has specific f c for the design of 55 MPa or greater). In addition, f c of cube 100 mm was converted to cube 150 mm by multiplying with 0.96 for NSC and 0.9 for HSC [74]. All the results of f st were for cylinder ∅150 mm × 300 mm. The f st results obtained from cylinders ∅150 mm × 150 mm and ∅100 mm × 200 mm were converted into the equivalent of f st obtained from cylinder ∅150 mm × 300 mm by multiplying with factors 0.93 and 0.91, respectively [75,76]. The database for FS was for the prism size 100 × 100 × 500 mm.

Gene Expression Programming
GEP is a branch of genetic programming (GP) and it was originated by Ferreira [77]. The GP is a method for solving problems that are not domain specific. It uses Darwinian reproduction and the survival of the fittest principles to solve problems. In order to obtain a solution whose length can change throughout a run, GP employs a parse tree structure. Function set, terminal set, fitness function, control parameters, and terminal condition are the five distinct elements of GEP. The first three components control the algorithm's search space, while the latter two components control the search's speed and quality [78]. A character string of fixed length is used to get a solution in the GEP algorithm. The parse trees of different sizes and shapes are then used to present the solution and these trees are called expression trees (ETs). The complex and nonlinear programs can be generated with the help of the multi-genic nature of GEP. Arithmetic operations are represented as a function set in each gene of the GEP (for example, +, −, ×, /, etc.), and constants and fixed-length variables are represented as a terminal set (for example, 1, 2, a, b). In general, the number of chromosomes controls how long the program will run. An error may be reduced but running time is lengthened by increasing the number of chromosomes. Figure 1 represents the ET with one gene, three head sizes, and a function set of +, −, and * . In order to obtain the mathematical formula, ET has to be read from left to right and top to bottom. The mathematical formula of ET in Figure 1 is ((C * MK) + (0.5 − SP)).  Figure 2 shows the flowchart of the GEP algorithm. Random generation of a c some with a fixed length for each individual is the starting point of the GEP alg The individuals are evaluated and chosen based on their fitness for reproductio process continues with a new individual for a few generations and stops when a s is found. Conversion in population is performed on the selected program by using operators, such as mutation, rotation, and crossover [77].

Artificial Neural Network
Artificial neural network (ANN) is a widely used artificial intelligence (AI) m There are many types of ANN such as radial basis function network, feedforward network (FNN), spiking neural network, etc. Among them, the widely used me FNN [80]. Single-Layer Perceptron (SLP) and Multi-Layer Perceptron (MLP) are tw  Figure 2 shows the flowchart of the GEP algorithm. Random generation of a chromosome with a fixed length for each individual is the starting point of the GEP algorithm. The individuals are evaluated and chosen based on their fitness for reproduction. This process continues with a new individual for a few generations and stops when a solution is found. Conversion in population is performed on the selected program by using genetic operators, such as mutation, rotation, and crossover [77].  Figure 2 shows the flowchart of the GEP algorithm. Random generatio some with a fixed length for each individual is the starting point of the G The individuals are evaluated and chosen based on their fitness for repro process continues with a new individual for a few generations and stops w is found. Conversion in population is performed on the selected program by operators, such as mutation, rotation, and crossover [77].

Artificial Neural Network
Artificial neural network (ANN) is a widely used artificial intelligence There are many types of ANN such as radial basis function network, feedf network (FNN), spiking neural network, etc. Among them, the widely u

Artificial Neural Network
Artificial neural network (ANN) is a widely used artificial intelligence (AI) method. There are many types of ANN such as radial basis function network, feedforward neural network (FNN), spiking neural network, etc. Among them, the widely used method is  [80]. Single-Layer Perceptron (SLP) and Multi-Layer Perceptron (MLP) are two types of FNN. Because of a single perceptron, SLP cannot execute non-linear problems; therefore, MLP is often used for non-linear problems [80].
An input layer, hidden layer(s), units (neurons), weights, an activation function, and an output layer are the typical components of an MLP. The input layer receives information from the outside environment. Without completing any calculations, the input layer sends these data to neurons in the hidden layer. The majority of a network's internal processing takes place in hidden layers, which are situated in between the input and output layers. The calculations are displayed to the environment outside by the output layer. The weights are used to connect adjoining layers and the function of the activation function is to decide how neurons will generate output value for the next layer [80]. Figure 3 shows a three-layer MLP with two inputs, one hidden layer with four hidden neurons, and two outputs.
Materials 2022, 15, x FOR PEER REVIEW processing takes place in hidden layers, which are situated in between the input put layers. The calculations are displayed to the environment outside by the outp The weights are used to connect adjoining layers and the function of the activati tion is to decide how neurons will generate output value for the next layer [80]. shows a three-layer MLP with two inputs, one hidden layer with four hidden n and two outputs. Three common forms of activation functions are linear transfer (purelin), hy tangent (tanh), and sigmoid (logsig). For problems with function fitting, linear functions work effectively. The ranges of the outputs of the tanh and logsig fu respectively, are −1 to 1 and 0 and 1. The fact that tanh can simulate input val contain negative, neutral, and positive numbers is one of its advantages over log There are a few advantages of ANN, such as: (1) it can model the relatively process and it does take outliers into account which makes its scope broad; (2) it c from examples and can build a relationship between dependent and independen bles [82]. More information about the ANN can be found in [81].

M5P Model Tree Algorithm
M5 algorithm was originally discovered by Quinlan [83] and the M5P algori is its expanded form. In order to handle enumerated attributes and missing va attributes, the M5P method was modified from M5. The M5P algorithm converts merated properties into binary variables prior to tree creation [85].
The illustration of the M5 algorithm is presented in Figure 4. The input dat vided into a number of sub-spaces, each of which contains data with shared charac (Figure 4a). To lessen a variation in the data inside a specific sub-space, linear re models are applied. The splitting process is then carried out on a number of no pending on data gathered from the preceding stage, and each node is separated ac to a certain attribute (Figure 4b). This stage enables the creation of a structure tha Three common forms of activation functions are linear transfer (purelin), hyperbolic tangent (tanh), and sigmoid (logsig). For problems with function fitting, linear transfer functions work effectively. The ranges of the outputs of the tanh and logsig functions, respectively, are −1 to 1 and 0 and 1. The fact that tanh can simulate input values that contain negative, neutral, and positive numbers is one of its advantages over logsig [81].
There are a few advantages of ANN, such as: (1) it can model the relatively complex process and it does take outliers into account which makes its scope broad; (2) it can learn from examples and can build a relationship between dependent and independent variables [82]. More information about the ANN can be found in [81].

M5P Model Tree Algorithm
M5 algorithm was originally discovered by Quinlan [83] and the M5P algorithm [84] is its expanded form. In order to handle enumerated attributes and missing values for attributes, the M5P method was modified from M5. The M5P algorithm converts all enumerated properties into binary variables prior to tree creation [85].
The illustration of the M5 algorithm is presented in Figure 4. The input data are divided into a number of sub-spaces, each of which contains data with shared characteristics (Figure 4a). To lessen a variation in the data inside a specific sub-space, linear regression models are applied. The splitting process is then carried out on a number of nodes depending on data gathered from the preceding stage, and each node is separated according to a certain attribute (Figure 4b). This stage enables the creation of a structure that resembles an inverted tree. When it comes to fresh data, they begin at the top of the tree's root, travel The input space is initially partitioned into a tree-like structure of different subspaces. At the node, the variability is calculated using the standard deviation of the values. Standard deviation reduction (SDR) is used to reduce the expected error at the node and helps to build a tree as follows: where is the standard deviation, are the sets produced by splitting node in accordance with a given attribute, and is the set of data that comes to the node [84]. Using a pruning strategy, the over-training issue is managed. However, the trimming procedure might result in jarring breaks between adjacent linear models. The final phase is the smoothing procedure to solve this issue. The final model of the leaf is created during the smoothing phase by combining all models from the leaf to the root. This filters the The input space is initially partitioned into a tree-like structure of different sub-spaces. At the node, the variability is calculated using the standard deviation of the values. Standard deviation reduction (SDR) is used to reduce the expected error at the node and helps to build a tree as follows: where sd is the standard deviation, S i are the sets produced by splitting node in accordance with a given attribute, and S is the set of data that comes to the node [84]. Using a pruning strategy, the over-training issue is managed. However, the trimming procedure might result in jarring breaks between adjacent linear models. The final phase is the smoothing procedure to solve this issue. The final model of the leaf is created during the smoothing phase by combining all models from the leaf to the root. This filters the estimated value of the leaf [23].

Random Forest
RF is a supervised ML method that comprises an ensemble of tree structures. ML techniques like bagging and random feature selection are used in RF [86]. In bagging, the bootstrap sample is generated using training data, and each tree is individually formed based on this sample. The estimate process then makes use of the average of the tree outputs [87]. A modified variation of bagging is RF. In RF, instead of selecting all features for a tree, a random subset of features is chosen. Due to its randomness, RF is resistant to overfitting and outperforms other ML methods like ANN and support vector machines [86].
RF has good generalization capability [86] and it provides a flexible framework with room for selecting objective functions (task-specific), various classes of splitting functions, or posterior models. Tree depth and the number of trees are the main hyperparameters in RF. The depth of a tree directly impacts the generalization ability of each tree, so its maximal allowable limit should be optimized [88]. An increasing number of trees helps to decrease prediction error by average out of noisy predictions. The schematic of the RF is shown in Figure 5. A detailed description of RF methodology can be found in [86].
R PEER REVIEW 9 of 37 outputs [87]. A modified variation of bagging is RF. In RF, instead of selecting all features for a tree, a random subset of features is chosen. Due to its randomness, RF is resistant to overfitting and outperforms other ML methods like ANN and support vector machines [86]. RF has good generalization capability [86] and it provides a flexible framework with room for selecting objective functions (task-specific), various classes of splitting functions, or posterior models. Tree depth and the number of trees are the main hyperparameters in RF. The depth of a tree directly impacts the generalization ability of each tree, so its maximal allowable limit should be optimized [88]. An increasing number of trees helps to decrease prediction error by average out of noisy predictions. The schematic of the RF is shown in Figure 5. A detailed description of RF methodology can be found in [86].

Model Development and Evaluation Criteria
The collected database was randomly divided into two sets: 70% and 30% for the training and testing sets, respectively. The ′ , , and FS of concrete incorporating MK were considered to be a function of the following input parameters while developing models: where ′ , , and FS are in MPa, while C, MK, FA, CA, and SP are in kg/m 3 . For GEP modelling, three models were developed for the mechanical properties of

Model Development and Evaluation Criteria
The collected database was randomly divided into two sets: 70% and 30% for the training and testing sets, respectively. The f c , f st , and FS of concrete incorporating MK were considered to be a function of the following input parameters while developing models: f c , f st , and FS = f (C, MK, w/c, FA, CA, SP, days) where f c , f st , and FS are in MPa, while C, MK, FA, CA, and SP are in kg/m 3 . For GEP modelling, three models were developed for the mechanical properties of concrete incorporating MK named: GEP I for f c , GEP II for f st , and GEP III for FS. The parameters used in the GEP algorithm for three GEP models are shown in Table 2. The Sqrt, Exp, Ln, Log, Inv, X2, X3, X4, X5, 3Rt, 4Rt, 5Rt denote square root, exponential, natural logarithm, inverse, X to the power of 2, X to the power of 3, X to the power of 4, X to the power of 5, cube root, quartic root, and quintic root, respectively. The input, output, and hidden layer specifications are the initial stage in creating an ANN model(s). All of the ANN models used in this study have one hidden layer, one input layer, and one output layer. All of the proposed ANN models have a 7-n-1 architecture. Trainlm, which changes bias and weight values in line with Levenberg-Marquardt (LM) optimization, was utilized for the training function [81]. Additionally, for the performance function and adaptation learning function, respectively, learngdm and mean squared error were used. The log-sigmoid was chosen as a transfer function in all created ANN models.
Three models were developed for the mechanical properties of concrete with MK named: ANN I for f c , ANN II for f st , and ANN III for FS. In this study, the number of epochs, max_fail, µ, and min_grad values were kept as 1000, 35, 0.001, and 1 × 10 −7 , respectively for all developed ANN models.
The M5P algorithm generates linear regression mathematical equations after making different classes of data. The general form of the M5P algorithm can be written as follow: Two M5P-based models were developed for concrete with MK inclusion named M5P II for f st and M5P III for FS. Similar to other GEP and ANN techniques used in this study, three models were developed for estimating the mechanical properties of concrete with MK by using RF named RF I for f c , RF II for f st , and RF III for FS.
For comparison purposes among models developed by GEP, ANN, M5P, and RF, graphical presentations of absolute error (AE) were drawn for both training and testing data sets. A horizontal line on an absolute error of 10 MPa was drawn for f c results and the percentages of data below it were mentioned. For f st and FS results, this operation was performed on an absolute error of 0.75 MPa.
For the four ML techniques, several trials were run in order to obtain a higher value of the R 2 and R, and lower values of relative squared error (RSE), mean absolute error (MAE), root mean squared error (RMSE) for both training and testing sets. Moreover, a performance index (ρ) was used to assess the model performance as a function of both R and relative root mean squared error (RRMSE).
where e is the average experimental value, n represents the total number of samples, and e i and p i are the experimental and predicted values, respectively. Poor performance is indicated for models with R 2 < 0.7 [89] while a model with R > 0.8 indicates a significant positive correlation between estimated and experimental outcomes [90]. The RSME, MAE, and RSE illustrate how accurate the proposed model is; a high value demonstrates how far the estimated results differ from the experimental results, whilst a low value demonstrates an acceptable level of accuracy in the estimated outcomes.
In addition to the aforementioned statistical indicators, the discrepancy ratio (DR) was also utilized to evaluate the performance of developed models. DR is expressed as: where all the terms are as described previously. A precise match between the estimated and real values is indicated by a DR of zero. Between actual and predicted values, a negative DR denotes underestimation and a positive DR denotes overestimation. In this study, in models developed for f c of concrete with MK, the accuracy is defined as a percentage of DR values fall in the range of −0.1 to 0.1; this range was also used by Benhood et al. [23]. However, for f st and FS database, the accuracy is defined as percentage of DR values that fall in the range of −0.05 to 0.05.

GEP I Model
To develop all GEP models in this study, GeneXproTools 5.0 software was used. Different researchers have used different parameters of the GEP algorithm in order to obtain a model with high accuracy and generalization capability as shown in Table 3. In this study, several runs were tried in order to obtain: (1) a relatively simple model by trying to minimize the number of genes and head size as a number of genes increases sub-ET size and head size increases the complexity within each gene. (2) A model that gives a high value of R 2 and R for training and testing sets and a low value of MAE, RMSE, RSE, and ρ. After several trials, the parameters of the optimum GEP model (GEP I) for f c of concrete with MK are given in Table 2 and generated ETs are shown in Figure 6. The mathematical formula was obtained from Figure 6 by following the procedure as mentioned in Section 3.1 and given as follows: where Table 3. Different parameters of the GEP algorithm used by researchers to obtain a reliable and robust model. The comparison of actual and estimated results for training and testing data sets by GEP I for f c of concrete with MK and their absolute error is shown in Figure 7. The R 2 value of 0.81 for the training set shows that GEP predicted values correlate well nonlinearly with actual results, while the R 2 value of 0.81 for the testing set indicates that the GEP algorithm can predict the output well by using unseen data as input, indicating its high generalization capability.
AE was also plotted in order to assess the performance of GEP prediction. It is shown in Figure 7c that a relatively large portion of training data (76.3%) are below the AE of 10 MPa with an average error or MAE of 7.1 MPa. For testing data, 74.2% of data lie below 10 MPa with an average error of 7.3 MPa. Training data that are below 10 MPa AE are slightly higher (2.8%) as compared with testing data.

ANN I Model
The ANN I model and other ANN models (ANN II and ANN III) developed in this study were trained by using MATLAB R2019a neural network toolbox. After several trials, the best ANN model was achieved with the architecture of 7-12-1 named ANN I. Figure 8a,b show that the value predicted by the ANN I model matches excellently with actual experimental results, having a value of R 2 equal to 0.94 for both training and testing sets. This value of R 2 for both data sets is 16% higher as compared with the value obtained from the GEP I model. The high R 2 value for the testing set indicates that ANN has a high generalization capacity and can predict the output based on unseen data once it is adequately trained on given input parameters and output. Moreover, the slope of regression for both data sets is close to the ideal fit (1 for the ideal case) indicating a slight difference between actual and estimated results.

ANN I Model
The ANN I model and other ANN models (ANN II and ANN III) developed in this study were trained by using MATLAB R2019a neural network toolbox. After several trials, the best ANN model was achieved with the architecture of 7-12-1 named ANN I. Figure 8a,b show that the value predicted by the ANN I model matches excellently with actual experimental results, having a value of R 2 equal to 0.94 for both training and testing sets. This value of R 2 for both data sets is 16% higher as compared with the value obtained from the GEP I model. The high R 2 value for the testing set indicates that ANN has a high generalization capacity and can predict the output based on unseen data once it is adequately trained on given input parameters and output. Moreover, the slope of regression for both data sets is close to the ideal fit (1 for the ideal case) indicating a slight difference between actual and estimated results. A very large portion of the training (92%) and testing (93%) data sets is below the AE of 10 MPa as shown in Figure 8c,d. The percentage of data below AE of 10 MPa obtained by ANN I is 20.6% and 25% higher for training and testing data sets as compared with GEP I. This large quantity of data below AE of 10 MPa shows that difference between actual and predicted results is small.

RF I
The RF I model and other RF models in this study (RF II and RF III) were developed by using WEKA version 3.9.5 (developed by University of Waikato, Hamilton, New Zealand). All the settings of parameters were kept as default for all the RF models developed in this study. Figure 9a,b present the comparison between predicted values obtained from RF I and experimental results for training and testing sets. An excellent value of R 2 (i.e., 0.99), which is approaching ideal condition 1, for both data sets indicates that RF is an excellent tool for the prediction of f c of concrete with MK and has a high capability to forecast output based on un-seen data set. The high prediction capability of RF I is also obvious from the slope of the regression line which is 0.95 and 0.96 for training and testing data sets, respectively. For both training and testing sets, the R 2 values given by RF I are 22.2% and 5.3% higher than that of GEP I and ANN I, respectively. A very large portion of the training (92%) and testing (93%) data sets is below the AE of 10 MPa as shown in Figure 8c,d. The percentage of data below AE of 10 MPa obtained by ANN I is 20.6% and 25% higher for training and testing data sets as compared with GEP I. This large quantity of data below AE of 10 MPa shows that difference between actual and predicted results is small.

RF I
The RF I model and other RF models in this study (RF II and RF III) were developed by using WEKA version 3.9.5 (developed by University of Waikato, Hamilton, New Zealand). All the settings of parameters were kept as default for all the RF models developed in this study. Figure 9a,b present the comparison between predicted values obtained from RF I and experimental results for training and testing sets. An excellent value of R 2 (i.e. 0.99), which is approaching ideal condition 1, for both data sets indicates that RF is an excellent tool for the prediction of ′ of concrete with MK and has a high capability to forecast output based on un-seen data set. The high prediction capability of RF I is also obvious from the slope of the regression line which is 0.95 and 0.96 for training and testing data sets, respectively. For both training and testing sets, the R 2 values given by RF I are 22.2% and 5.3% higher than that of GEP I and ANN I, respectively. The percentages of data of AE below 10 MPa are almost 100% for both data sets as shown in Figure 9c,d. This depicts that the difference between actual and predicted results is not high and is less as compared with GEP I and ANN I. This is also clear from the value of average error for both data sets which is 1.45 MPa for training and 1.32 MPa for the testing set.

Comparison of GEP I, ANN I, and RF I
It is clear from Table 4 that the RF I predicted values correlated excellently with actual results, with values of R equal to 0.997 and 0.996 for training and testing sets, respectively. With respect to the values of R for both data sets, the order of correlation between predicted and actual results for developed models was RF I > ANN I > GEP I. In case of statistical errors (i.e., RMSE, RRMSE, RSE) and ρ value, the order of developed models for f c was RF I < ANN I < GEP I for the training data set. This shows that the RF I model has high performance (as indicated by the low value of ρ) and predicted results are close to experimental data followed by ANN I and GEP I. In the case of the testing set, for values of RRMSE, RSE, and ρ, the order of developed models is similar to that observed in the training set.
In addition, Figure 10 shows that, for the training set, the accuracy of RF I, ANN I, and GEP I is 99.43%, 90.23%, and 75.43%, respectively, as measured by the percentage of DR values that fall in the range of −0.1 to 0.1. In the case of the testing set, the accuracy is 99.66%, 94.29%, and 70.14% for RF I, ANN I, and GEP I, respectively. The percentages of data of AE below 10 MPa are almost 100% for both data sets as shown in Figure 9c,d. This depicts that the difference between actual and predicted results is not high and is less as compared with GEP I and ANN I. This is also clear from the value of average error for both data sets which is 1.45 MPa for training and 1.32 MPa for the testing set.

Comparison of GEP I, ANN I, and RF I
It is clear from Table 4 that the RF I predicted values correlated excellently with actua results, with values of R equal to 0.997 and 0.996 for training and testing sets, respectively With respect to the values of R for both data sets, the order of correlation between pre dicted and actual results for developed models was RF I > ANN I > GEP I. In case of sta tistical errors (i.e., RMSE, RRMSE, RSE) and value, the order of developed models for ′ was RF I < ANN I < GEP I for the training data set. This shows that the RF I model has high performance (as indicated by the low value of ) and predicted results are close to experimental data followed by ANN I and GEP I. In the case of the testing set, for values of RRMSE, RSE, and , the order of developed models is similar to that observed in the training set.   The high performance and prediction capability of ANN I over GEP I is in agreement with Nazari et al. [93] and Yu et al. [94]. Nazari et al. [93] predicted the water absorption (%) of HSC containing TiO 2 nanoparticles by using two ANN-based models (they named ANN I and ANN II) and two GEP-based models (GEP I and GEP II). For the training set, ANN I and ANN II showed the values of R 2 equal to 0.99 and 0.97, respectively, and GEP I and GEP II gave the values of R 2 equal to 0.91 and 0.85, respectively. In the case of the testing set, ANN I and GEP I gave the value of R 2 equal to 0.96 and 0.9, respectively, while these values were 0.93 and 0.85 in the case of ANN II and GEP II, respectively. Yu et al. [94] predicted the degradation of elastic modulus induced by the alkali-silica reaction by using soft computing techniques, including ANN and GEP. For the training set, ANN and GEP showed values of R equal to 0.98 and 0.86, respectively. For the testing set, these values were 0.93 and 0.91 by using ANN and GEP, respectively.
In addition, Figure 10 shows that, for the training set, the accuracy of RF I, ANN I, and GEP I is 99.43%, 90.23%, and 75.43%, respectively, as measured by the percentage of DR values that fall in the range of −0.1 to 0.1. In the case of the testing set, the accuracy is 99.66%, 94.29%, and 70.14% for RF I, ANN I, and GEP I, respectively. The high performance and prediction capability of ANN I over GEP I is in agreement with Nazari et al. [93] and Yu et al. [94]. Nazari et al. [93] predicted the water absorption (%) of HSC containing TiO2 nanoparticles by using two ANN-based models (they named ANN I and ANN II) and two GEP-based models (GEP I and GEP II). For the training set, ANN I and ANN II showed the values of R 2 equal to 0.99 and 0.97, respectively, and GEP I and GEP II gave the values of R 2 equal to 0.91 and 0.85, respectively. In the case of the testing set, ANN I and GEP I gave the value of R 2 equal to 0.96 and 0.9, respectively, while these values were 0.93 and 0.85 in the case of ANN II and GEP II, respectively. Yu et al. [94] predicted the degradation of elastic modulus induced by the alkali-silica reaction by using soft computing techniques, including ANN and GEP. For the training set, ANN and GEP showed values of R equal to 0.98 and 0.86, respectively. For the testing set, these values were 0.93 and 0.91 by using ANN and GEP, respectively.
The high prediction and generalization capability of RF I as compared with GEP I is in agreement with the study by Mohsin et al. [95]. They modelled the ′ of fly ash based geopolymer concrete by using GEP and RF. In case of RF model, the values of R for training and testing set were 0.98 and 0.99, respectively, while in case of the GEP model, these values were 0.86 and 0.96 for the training and testing sets, respectively.

GEP II
The parameters of the optimal GEP model (GEP II) are shown in Table 2 and the developed expression tree is given in Figure 11. The empirical expression decoded from the expression tree is given as follows: The high prediction and generalization capability of RF I as compared with GEP I is in agreement with the study by Mohsin et al. [95]. They modelled the f c of fly ash based geopolymer concrete by using GEP and RF. In case of RF model, the values of R for training and testing set were 0.98 and 0.99, respectively, while in case of the GEP model, these values were 0.86 and 0.96 for the training and testing sets, respectively.

GEP II
The parameters of the optimal GEP model (GEP II) are shown in Table 2 and the developed expression tree is given in Figure 11. The empirical expression decoded from the expression tree is given as follows: The comparison of GEP II predicted values and actual results along with AE for both training and testing sets is shown in Figure 12. A high correlation between estimated and experimental values was found as indicated by R 2 values which were 0.86 and 0.9 for training and testing sets. A slight difference of 0.04 was also observed between R 2 values of training and testing sets. The slope of the regression line is high, i.e., 0.89 and 0.87 for training and testing sets, indicating a high correlation. The comparison of GEP II predicted values and actual results along with AE for both training and testing sets is shown in Figure 12. A high correlation between estimated and experimental values was found as indicated by R 2 values which were 0.86 and 0.9 for training and testing sets. A slight difference of 0.04 was also observed between R 2 values of training and testing sets. The slope of the regression line is high, i.e., 0.89 and 0.87 for training and testing sets, indicating a high correlation. In addition, AE is plotted and the horizontal line is drawn on AE of 0.75 MPa, and the percentage of data below this line is shown in order to give more insight into the capability of the model to predict output close to real values as shown in Figure 12c,d. It is clear from the figures that the majority of data are below the AE of 0.75 MPa for both data sets (i.e., 93% and 95% for training and testing sets, respectively).

ANN II
After several trials, the best ANN model (ANN II) was obtained with an architecture of 7-11-1. The plot of predicted and actual values for the training and testing sets are shown in Figure 13a,b. Values of R 2 are relatively close to 1 for training (R 2 = 0.92) and testing (R 2 = 0.96) sets. These higher values of R 2 for both data sets indicate that ANN can recognize the relationship between input and output variables well and produce output by using unseen data with high accuracy. As compared with GEP II, these values were 6.98% and 6.7% higher for training and testing sets, respectively. Moreover, the slope of the regression line was slightly higher for the training set and considerably higher for the testing set as compared with GEP II. In addition, AE is plotted and the horizontal line is drawn on AE of 0.75 MPa, and the percentage of data below this line is shown in order to give more insight into the capability of the model to predict output close to real values as shown in Figure 12c,d. It is clear from the figures that the majority of data are below the AE of 0.75 MPa for both data sets (i.e., 93% and 95% for training and testing sets, respectively).

ANN II
After several trials, the best ANN model (ANN II) was obtained with an architecture of 7-11-1. The plot of predicted and actual values for the training and testing sets are shown in Figure 13a,b. Values of R 2 are relatively close to 1 for training (R 2 = 0.92) and testing (R 2 = 0.96) sets. These higher values of R 2 for both data sets indicate that ANN can recognize the relationship between input and output variables well and produce output by using unseen data with high accuracy. As compared with GEP II, these values were 6.98% and 6.7% higher for training and testing sets, respectively. Moreover, the slope of the regression line was slightly higher for the training set and considerably higher for the testing set as compared with GEP II.
Very large data are below the AE of 0.75 MPa for the training set (97%) and testing set (100%), as shown in Figure 13c,d. This shows that ANN II was trained well and predicted the output which was very close to the experimental output. Very large data are below the AE of 0.75 MPa for the training set (97%) and testing set (100%), as shown in Figure 13c,d. This shows that ANN II was trained well and predicted the output which was very close to the experimental output.

M5P II
In this study, all M5P models were developed by using the Waikato Environment for Knowledge Analysis (WEKA) software version 3.9.5 (Hamilton, New Zealand). The model trees are generated as shown in Figure 14. The term LM at the tree leaves represents the linear model identified by the M5P algorithm. The corresponding coefficients for linear models developed by M5P II based on Equation (3) are shown in Table 5.

M5P II
In this study, all M5P models were developed by using the Waikato Environment for Knowledge Analysis (WEKA) software version 3.9.5 (Hamilton, New Zealand). The model trees are generated as shown in Figure 14. The term LM at the tree leaves represents the linear model identified by the M5P algorithm. The corresponding coefficients for linear models developed by M5P II based on Equation (3) are shown in Table 5. Figure 15a,b depict the correlation between M5P II predicted values and experimental results. The values of R 2 were 0.88 and 0.86 for training and testing sets, respectively. For the training set, this value was slightly higher as compared with GEP II but lower than that of ANN II. In the case of the testing set, it was lower than both GEP II and ANN II. In term of the percentage of data below AE of 0.75 MPa, M5P II perform slightly better as compared with GEP II for both training and testing sets. Figure 16 shows the comparison of actual and predicted values by RF II along with AE for training and testing sets. Excellent correlation was observed for both data sets, with R 2 equal to 0.98 and 0.99 for training and testing sets. This impressive value of R 2 for both data sets indicates that RF is a potential candidate for predicting f st of concrete with MK with high accuracy and generalization capability. Moreover, almost all the data were below the AE of 0.75 MPa for both data sets. The average error for both data sets was approaching zero, indicating that the difference between actual and predicted values was very small.   Figure 15a,b depict the correlation between M5P II predicted values and mental results. The values of R 2 were 0.88 and 0.86 for training and testing sets tively. For the training set, this value was slightly higher as compared with GE lower than that of ANN II. In the case of the testing set, it was lower than both GE ANN II. In term of the percentage of data below AE of 0.75 MPa, M5P II perform better as compared with GEP II for both training and testing sets.    Figure 16 shows the comparison of actual and predicted values by RF II along with AE for training and testing sets. Excellent correlation was observed for both data sets, with R 2 equal to 0.98 and 0.99 for training and testing sets. This impressive value of R 2 for both data sets indicates that RF is a potential candidate for predicting of concrete with MK with high accuracy and generalization capability. Moreover, almost all the data were below the AE of 0.75 MPa for both data sets. The average error for both data sets was approaching zero, indicating that the difference between actual and predicted values was very small. For the training set, the values of R 2 obtained by RF II were 13.95%, 6.5%, and 11.4% higher as compared with GEP II, ANN II, and M5P II, respectively. In the case of the testing set, R 2 of RF II was 10%, 3.1%, and 15.11% higher than that of GEP II, ANN II, and M5P II, respectively. Table 6 shows statistical measures of different models developed for predicting f st of concrete with MK. For the training set, the order of value of R for different models was RF II > ANN II > M5P II > GEP II. In the case of the value of RMSE, RRMSE, RSE, and ρ , this order was reversed, i.e., RF II < ANN II < M5P II < GEP II. A slight change was observed in the case of the testing set, in which M5P II and GEP II replaced each other's position, respective to their positions in the case of the training set. In addition, Figure 17 shows that the accuracy of GEP II, ANN II, M5P II, and RF II is 69.45%, 89.58%, 81.94%, and 97.23%, respectively for the training set, while it is 77.05%, 95.08%, 68.86%, and 100%, respectively for the testing set. The order of accuracy given by DR for all models in Figure 17 is aligned with the order of accuracy of all models by using different statistical measures as shown in Table 6. Overall, RF performed better as compared with the other three machine learning techniques for predicting f st of concrete with MK and these results are in agreement with the prediction results of f c database. For the training set, the values of R 2 obtained by RF II were 13.95%, 6.5%, and 11.4% higher as compared with GEP II, ANN II, and M5P II, respectively. In the case of the testing set, R 2 of RF II was 10%, 3.1%, and 15.11% higher than that of GEP II, ANN II, and M5P II, respectively. Table 6 shows statistical measures of different models developed for predicting of concrete with MK. For the training set, the order of value of R for different models was RF II > ANN II > M5P II > GEP II. In the case of the value of RMSE, RRMSE, RSE, and this order was reversed, i.e., RF II < ANN II < M5P II < GEP II. A slight change was observed in the case of the testing set, in which M5P II and GEP II replaced each other's position, respective to their positions in the case of the training set. In addition, Figure 17 shows that the accuracy of GEP II, ANN II, M5P II, and RF II is 69.45%, 89.58%, 81.94% and 97.23%, respectively for the training set, while it is 77.05%, 95.08%, 68.86%, and 100% respectively for the testing set. The order of accuracy given by DR for all models in Figure  17 is aligned with the order of accuracy of all models by using different statistica measures as shown in Table 6. Overall, RF performed better as compared with the other three machine learning techniques for predicting of concrete with MK and these results are in agreement with the prediction results of ′ database.  The optimal parameters found after several trials for predicting FS of concrete with MK are shown in Table 2 with the name GEP III. The ET developed by GEP III is given in Figure 18 and decoded mathematical equation from the figure is as follows:

GEP III
The optimal parameters found after several trials for predicting FS of concrete with MK are shown in Table 2 with the name GEP III. The ET developed by GEP III is given in Figure 18 and decoded mathematical equation from the figure is as follows: where A = The relationship between GEP III predicted values and actual values along with AE for both training and testing sets is shown in Figure 19. The values of R 2 were 0.88 and 0.86 for training and testing sets, respectively. In the case of the training set, 80% of data are below AE of 0.75 MPa, while it is 89% in the case of the testing set, as shown in Figure 19c,d.

ANN III
Similar to ANN II, the best accuracy for ANN III was obtained with the architecture of 7-11-1. As shown in Figure 20a,b, the value of R 2 = 0.95 for both training and testing data sets shows that ANN III trained well with given inputs and outputs and estimated outputs by using unseen input parameters with high accuracy. This value of R 2 was 7.95% and 10.5% higher as compared with training and testing sets of GEP III, respectively. Significantly less error was noted in the case of the testing set as 100% of data are below the AE of 0.75 MPa, showing the high generalization ability of ANN III, as shown in Figure 20d.

M5P III
By using the default setting for parameters of WEKA software, only one linear model was obtained (due to the small database) for FS whose empirical expression is given below: The prediction capability of M5P III is observed to be relatively poor as compared with GEP III and ANN III as indicated by Figure 21a,b which shows that R 2 values for training and testing sets are 0.73 and 0.76, respectively. For the training set, this value was 17% and 23.16% lower as compared with GEP III and ANN III, respectively, while for the testing set, this value was 11.62% and 20% lower as compared with GEP III and ANN III, respectively.
In addition, the relatively low performance of M5P III can be observed in Figure 21c, which shows that 48% of the data are greater than the AE of 0.75 MPa. The relationship between GEP III predicted values and actual values along with AE for both training and testing sets is shown in Figure 19. The values of R 2 were 0.88 and 0.86 for training and testing sets, respectively. In the case of the training set, 80% of data are below AE of 0.75 MPa, while it is 89% in the case of the testing set, as shown in Figure  19c,d.

RF III
As compared with other ML techniques for modelling FS of concrete with MK, excellent performance was observed by RF III. Figure 22a,b show an excellent correlation between predicted and experimental values, i.e., for both training and testing sets, the value of R 2 = 0.98. This shows that RF III is a highly accurate, reliable, precise, and robust model. Figure 22c,d further strengthen the high accuracy of RF III, which shows that 100% of data are below 0.75 MPa for both data sets with low average errors.

ANN III
Similar to ANN II, the best accuracy for ANN III was obtained with the architecture of 7-11-1. As shown in Figure 20a,b, the value of R 2 = 0.95 for both training and testing data sets shows that ANN III trained well with given inputs and outputs and estimated outputs by using unseen input parameters with high accuracy. This value of R 2 was 7.95% and 10.5% higher as compared with training and testing sets of GEP III, respectively. Significantly less error was noted in the case of the testing set as 100% of data are below the AE of 0.75 MPa, showing the high generalization ability of ANN III, as shown in Figure  20d.

M5P III
By using the default setting for parameters of WEKA software, only one linear model was obtained (due to the small database) for FS whose empirical expression is given below:

RF III
As compared with other ML techniques for modelling FS of concrete with MK, excellent performance was observed by RF III. Figure 22a,b show an excellent correlation between predicted and experimental values, i.e., for both training and testing sets, the value of R 2 = 0.98. This shows that RF III is a highly accurate, reliable, precise, and robust model. Figure 22c,d further strengthen the high accuracy of RF III, which shows that 100% of data are below 0.75 MPa for both data sets with low average errors. For both training and testing sets, the order of R between predicted and measured values by different models is RF III > ANN III > GEP III > M5P III, as shown in Table 7. The same order was observed by the value of R 2 as discussed previously and by the per-

Comparison of GEP III, ANN III, M5P III, and RF III
For both training and testing sets, the order of R between predicted and measured values by different models is RF III > ANN III > GEP III > M5P III, as shown in Table 7. The same order was observed by the value of R 2 as discussed previously and by the percentage of DR that fall in the range of −0.05 to 0.05, as shown in Figure 23. In the case of the value of RMSE, RRMSE, RSE, and ρ, the order of models was RF III < ANN III < GEP III < M5P III for both data sets.  Table 7. The same order was observed by the value of R 2 as discussed previously and by the percentage of DR that fall in the range of −0.05 to 0.05, as shown in Figure 23. In the case of the value of RMSE, RRMSE, RSE, and , the order of models was RF III < ANN III < GEP III < M5P III for both data sets.

Sensitivity and Parametric Analysis
In order to find out the relative contribution of input parameters on outputs, sensitivity analysis (SA) was conducted by using the model proposed by Gandomi et al. [96] and as given in Equations (16) and (17).
where, f min (x i ) and f max (x i ) are the minimum and maximum predicted outputs based on ith input variable, in which other input variables are kept constant at their mean values. In order to calculate the variation in f c , f st , and FS by changing MK content and days, parametric analysis (PA) was carried out. The PA was performed by observing a change in f c , f st , and FS by increasing MK content and age of specimen from its minimum to maximum value while keeping all other input parameters at their mean values. The comparative study showed that RF performed better than the other modelling techniques; however, in this study, both sensitivity and parametric analysis was carried out by using GEP due to its convenience [79]. It is important to note that the accuracy of GEP models for the prediction of all three mechanical properties was high; therefore, it can be used to explore materials characteristics through sensitivity and parametric analysis.
For both f c and f st , the w/b ratio seemed to be the most influential parameter, followed by days, SP, C, MK, FA, and CA, as shown in Figure 24. However, in the case of FS, the number of days is the most influential parameter, followed by SP, w/b ratio, C, MK, FA, and CA. on i th input variable, in which other input variables are kept constant at their mean In order to calculate the variation in ′ , , and FS by changing MK con days, parametric analysis (PA) was carried out. The PA was performed by obs change in ′ , , and FS by increasing MK content and age of specimen from mum to maximum value while keeping all other input parameters at their mean The comparative study showed that RF performed better than the other modelli niques; however, in this study, both sensitivity and parametric analysis was car by using GEP due to its convenience [79]. It is important to note that the accuracy models for the prediction of all three mechanical properties was high; therefore, used to explore materials characteristics through sensitivity and parametric analy For both ′ and , the w/b ratio seemed to be the most influential param lowed by days, SP, C, MK, FA, and CA, as shown in Figure 24. However, in the ca the number of days is the most influential parameter, followed by SP, w/b ratio FA, and CA.  The PA in Figure 25 shows the variation in mechanical properties of concrete by increasing MK content from its minimum to maximum value. The development of concrete strength with MK incorporation can be attributed to: (i) pozzolanic reaction of MK with calcium hydroxide, (ii) acceleration of cement hydration, and (iii) the filling effect due to MK particles [15]. In mortar with MK, the formation of alumina phases such as C 2 ASH 8 is responsible for higher strength at early ages [97]. The rapid and early pozzolanic reaction of MK with CH may decrease the initial and final setting times of concrete with MK [55]. Figure 25a depicts that for up to 35 kg/m 3 addition of MK (which is about 10% cement replacement), f c increases linearly. After further increasing MK content up to 69 kg/m 3 (about 15% cement replacement), the f c increases but nonlinearly. MK content from 69 to 105 kg/m 3 (about 25% cement placement) increases the f c nonlinearly but with slower rate as compared with the 69 kg/m 3 . Further increase in MK beyond 105 kg/m 3 does not increase f c significantly. Similar to our results, Rahmat et al. [15] investigated f c of different SCC with MK incorporation (0-20%) at different w/b ratios and at different curing ages and observed that the most remarkable strength developments were found with 10-15% MK replacement. Hamdy et al. [39] investigated six different proportions of 0, 10, 15, 20, 30, 40, and 50% cement replacement with MK for high-strength concrete and observed f c at 3 and 7 days. They concluded that maximum f c was observed at 15% MK replacement. After the 30%, the f c was observed to be lowered as compared with plain concrete. The decrease in f c by increasing MK content beyond optimum content may be due to the reason that higher MK content decreases the CaO/SiO 2 ratio [72], which results in the higher requirement of SP [34] and clinker dilution effect due to partial replacement of cement with MK [46]. Moreover, at a low w/c ratio, increasing MK% beyond 15% decreased the f c and f st , as, in this case, less calcium hydroxide is available for reaction with MK [39]. Figure 25b shows that from 10-40 kg/m 3 MK content (about 2.5-10% cement replacement), f st increases almost linearly. By further increasing MK content up to 60 kg/m 3 (about 15% cement replacement), f st increases non-linearly, but with a faster rate as compared with further increase in MK content. Above 100 kg/m 3 (or about 25% cement replacement), no significant improvement in f st is observed. Rahmat et al. [15] also observed that SCC showed better f st with 10-15% MK as cement replacement. Kannan et al. [53] incorporated five different proportions of 5, 10, 15, 20, 25, and 30% MK as particle cement replacement and observed a maximum value of f st at 20% and then a decrease in f st beyond this percentage. The trend of f st with MK was similar to that observed for f c with MK. Figure 25c depicts that from the incorporation of MK up to 50 kg/m 3 (about 10% cement replacement), FS increases significantly. From 50 to 100 kg/m 3 MK content (about 10 to 20% cement replacement), the increase in FS was not significant. Lenka et al. [4] also observed that concrete gave a maximum performance in terms of FS at 10% MK inclusion as cement replacement.  Figure 26a shows that the rate of strength development at early ages (up to 7 days) is fast, which can be attributed to the fast pozzolanic reaction of MK [98]. Bai et al. [99] observed an up to 92% increase in early age strength with the incorporation of 5% MK as compared with plain concrete. Erhan et al. [16] found that concrete strength with MK at early ages (1-7 days) was 5-23% greater as compared with plain concrete, depending on MK replacement level and w/b ratio. From 14-60 days, the ′ development of concrete with MK was slow and no significant strength development was observed after 90 days. Figure 26b shows that, similar to ′ , the development of at an early age (1 day) is very high as compared with further increase in time. From 3 to 90 days, increased but with a much slower rate as compared with day 1. Similar to ′ and , the rate of FS  Figure 26a shows that the rate of strength development at early ages (up to 7 days) is fast, which can be attributed to the fast pozzolanic reaction of MK [98]. Bai et al. [99] observed an up to 92% increase in early age strength with the incorporation of 5% MK as compared with plain concrete. Erhan et al. [16] found that concrete strength with MK at early ages (1-7 days) was 5-23% greater as compared with plain concrete, depending on MK replacement level and w/b ratio. From 14-60 days, the f c development of concrete with MK was slow and no significant strength development was observed after 90 days. Figure 26b shows that, similar to f c , the development of f st at an early age (1 day) is very high as compared with further increase in time. From 3 to 90 days, f st increased but with a much slower rate as compared with day 1. Similar to f c and f st , the rate of FS development is higher in the first 7 days as compared with the rest of the days, as shown in Figure 26c.

Conclusions
Compressive strength ( ′ ), splitting tensile strength ( ), and flexural strength (FS) of concrete are the parameters of the design in many codes. An accurate and reliable estimation of these parameters can save cost and time, and help in scheduling activities such as formwork removal. In this study, ′ , , and FS of concrete with the incorporation of metakaolin (MK) as partial cement replacement were modelled using four machine learning (ML) techniques: gene expression programming (GEP), artificial neural network (ANN), M5P model tree algorithm, and random forest (RF). For this purpose, a comprehensive database was gathered from peer-reviewed published documents. The database used in the modelling was comprised of 982 data points for ′ , 204 data points for , and 63 data samples for FS of concrete with MK. For all three databases, the input parameters were cement, MK, w/b, fine and coarse aggregates, superplasticizer, and age of a specimen in days. Many statistical metrics were used to compare the predictive performance of ML techniques used in this study. In the end, sensitivity and parametric analysis (PA) was performed. Based on the application of GEP, ANN, M5P, and RF for predicting ′ , , and FS of concrete with MK, the following conclusions can be drawn: 1. For modelling ′ of concrete with MK, RF I (R 2 = 0.99) showed excellent predictive capability followed by ANN I (R 2 = 0.94) and GEP I (R 2 = 0.81) for both training and testing sets. These results were also supported by other statistical metrics such as R, RMSE, RSE, MAE, DR, and .

Conclusions
Compressive strength ( f c ), splitting tensile strength ( f st ), and flexural strength (FS) of concrete are the parameters of the design in many codes. An accurate and reliable estimation of these parameters can save cost and time, and help in scheduling activities such as formwork removal. In this study, f c , f st , and FS of concrete with the incorporation of metakaolin (MK) as partial cement replacement were modelled using four machine learning (ML) techniques: gene expression programming (GEP), artificial neural network (ANN), M5P model tree algorithm, and random forest (RF). For this purpose, a comprehensive database was gathered from peer-reviewed published documents. The database used in the modelling was comprised of 982 data points for f c , 204 data points for f st , and 63 data samples for FS of concrete with MK. For all three databases, the input parameters were cement, MK, w/b, fine and coarse aggregates, superplasticizer, and age of a specimen in days. Many statistical metrics were used to compare the predictive performance of ML techniques used in this study. In the end, sensitivity and parametric analysis (PA) was performed. Based on the application of GEP, ANN, M5P, and RF for predicting f c , f st , and FS of concrete with MK, the following conclusions can be drawn: 1.
For modelling f c of concrete with MK, RF I (R 2 = 0.99) showed excellent predictive capability followed by ANN I (R 2 = 0.94) and GEP I (R 2 = 0.81) for both training and