Hybrid ensemble soft computing approach for predicting penetration rate of tunnel boring machine in a rock environment

Aidhn Brdhn, Nvid Krdni, Ansu GuhRy, Avijit Burmn, Pijush Smui,Ynmei Zhng

a National Institute of Technology Patna, Patna, Bihar, India

b Civil and Infrastructure Engineering Discipline, School of Engineering, Royal Melbourne Institute of Technology (RMIT), Victoria, 3001, Australia

c Department of Civil Engineering, BITS-Pilani Hyderabad Campus, Hyderabad, Telangana, India

d College of Aerospace Engineering, Chongqing University, Chongqing, 400044, China

Keywords:Tunnel boring machine (TBM)Rate of penetration (ROP)Artificial intelligence Artificial neural network (ANN)Ensemble modelling

ABSTRACT This study implements a hybrid ensemble machine learning method for forecasting the rate of penetration (ROP) of tunnel boring machine (TBM), which is becoming a prerequisite for reliable cost assessment and project scheduling in tunnelling and underground projects in a rock environment. For this purpose,a sum of 185 datasets was collected from the literature and used to predict the ROP of TBM.Initially, the main dataset was utilised to construct and validate four conventional soft computing (CSC)models, i.e. minimax probability machine regression, relevance vector machine, extreme learning machine, and functional network. Consequently, the estimated outputs of CSC models were united and trained using an artificial neural network (ANN) to construct a hybrid ensemble model (HENSM). The outcomes of the proposed HENSM are superior to other CSC models employed in this study.Based on the experimental results(training RMSE=0.0283 and testing RMSE=0.0418),the newly proposed HENSM is potential to assist engineers in predicting ROP of TBM in the design phase of tunnelling and underground projects.

1. Introduction

In the field of civil engineering, there is an endless growing demand due to rapid industrialisation. The construction industry plays a vital role in the development of the country in the field of transportation,hydropower,mining,and underwater construction.In several areas such as rail, metro, and hydropower, construction projects are underway across the world. In recent decades, rapid industrialisation is ongoing on a fast-track basis to reduce the overall development costs and timely completion of projects. In most heavy and civil engineering projects, tunnelling works are now being carried out to avoid open excavation and blasting(Mahdevari et al., 2014). Even though the drilling and blasting method is a well-developed technology, professionals are implementing mechanical excavation in underground construction(Armaghani et al., 2017). A tunnel boring machine (TBM) is an essential tool used for tunnelling in underground projects. TBMs offer copious advantages over the drilling and blasting method(Javad and Narges,2010).It should be noted that TBMs are gigantic machines and during the course of excavation,TBMs are extremely sensitive to conditions of rock mass throughout the path (Zhou et al., 2021). Though TBM is the best option for tunnelling in a rocky medium, insufficient information and unknown rock mass conditions can lead to a decline in the safety and efficiency of the tunnelling operation. Therefore, proper planning and accurate estimation are required to adopt the best construction techniques and to predict the performance of TBMs in a rock environment(Zhou et al., 2021).

Prediction of TBM performance is an important factor to estimate project duration, project cost, project schedule, and procurement schedule of bulk materials (Yagiz et al., 2009).Theoretically, the performance of a TBM is measured in terms of rate of penetration(ROP),which is the ratio of excavation distance to continuous operation time of TBM (Mahdevari et al., 2014).Several factors,such as intact rock properties,rock mass properties,cutter geometry, and specification of TBMs, influence the performance of TBMs(Zeng et al.,2021).Proper application of TBM to any tunnelling projects depends on the parameters mentioned above,and hence a detailed assessment is required to compute the ROP of TBMs.

Predicting the performance of TBMs in a rocky medium is a complex task compared to that in a soil medium.Thus,researchers investigated the inter-relationship between the rock properties and ROP of TBMs (Graham,1976; Farmer and Glossop,1980; Cassinelli et al.,1982). Based on the findings, it is clear that there is neither an exclusive way nor a fixed technique for predicting the ROP of TBM in practice(Zhang et al.,2020a).Methods implemented by the researchers can be broadly categorised into three clusters: (i)theoretical and empirical, (ii) statistical, and (iii) computational.These methods are developed from laboratory experiments and field observation, a set of mathematical rules, and several soft computing techniques, respectively. As mentioned by Gao et al.(2020), TBM’s ROP estimation using theoretical and empirical models give relatively low accuracy,and hence new methods need to be developed. On the other hand, a comprehensive review of previous studies reveals that the performance of statistical models(i.e.simple and multiple linear regression models)is relatively low compared to the theoretical and empirical models (Koopialipoor et al., 2019; Shahrour and Zhang, 2021; Zhang et al., 2021a). Thus,recent studies have resorted to soft computing techniques as potential alternatives to improve the reliability and exactness of the existing empirical and statistical models (Grima et al., 2000;Benardos and Kaliampakos, 2004; Yagiz et al., 2009; Javad and Narges, 2010; Yagiz and Karahan, 2011; Martins and Miranda,2013; Mobarra et al., 2013; Ghasemi et al., 2014; Mahdevari et al.,2014; Salimi et al., 2016; Armaghani et al., 2017, 2018; Fattahi,2019; Koopialipoor et al., 2019; Xu et al., 2019; Gao et al., 2020;Zhang et al., 2020a, b; Li et al., 2021; Zhou et al., 2021).

Soft computing algorithms, with their ability in nonlinear modelling,provide feasible tools for understanding and simulating diverse complicated problems (Chen et al., 2020; Goh et al., 2020;Zhang and Goh, 2013, 2016), including predictions of the performance of TBMs (Goh et al., 2017, 2018; Xiang et al., 2018). In the recent past, several soft computing algorithms, such as neural network (NN), artificial NN (ANN),adaptive neuro-fuzzy inference system (ANFIS), support vector regression (SVR), support vector machine (SVM), gene expression programming (GEP), relevance vector regression(RVR),relevance vector machine(RVM),deep NN(DNN), long-short-term memory (LSTM), classification and regression trees (CART), k-nearest neighbour (KNN), and chisquared automatic interaction detection (CHAID) (Grima et al.,2000; Benardos and Kaliampakos, 2004; Mobarra et al., 2013;Ghasemi et al., 2014; Mahdevari et al., 2014; Salimi et al., 2016;Armaghani et al., 2018; Koopialipoor et al., 2019; Xu et al., 2019;Gao et al.,2020;Zhang et al.,2020a;Li et al.,2021),have been used successfully for predicting the ROP of TBMs.Successful application of these methods can also be found in the literature.Besides,many other conventional soft computing(CSC)models,such as minimax probability machine regression (MPMR), functional network (FN),Gaussian process regression,and extreme learning machine(ELM),have been used successfully in different engineering disciplines for predicting the desired output(s)(Rajasekaran,2004;Liu et al.,2015;Jagan et al., 2019; Kumar and Samui, 2019; Samui, 2019; Ahmed et al., 2020; Mahmoud et al., 2020; Raja and Shukla, 2020, 2021a;Kumar et al., 2021a, b).

For estimating the performance of TBM, several CSC models have been used, including ANN/NN (Benardos and Kaliampakos,2004; Yagiz et al., 2009; Javad and Narges, 2010; Mobarra et al.,2013; Armaghani et al., 2017; Koopialipoor et al., 2019; Xu et al.,2019), ANFIS (Salimi et al., 2016), SVM/SVR (Martins and Miranda,2013; Mahdevari et al., 2014; Salimi et al., 2016; Xu et al., 2019),least-squares support vector machine (LS-SVM) (Zhang et al.,2020a), GEP (Armaghani et al., 2018), RVM/RVR (Fattahi, 2019;Zhang et al., 2020a), and DNN (Koopialipoor et al., 2019). Xu et al.(2019) also used KNN, CHAID, and CART to predict the ROP of TBM and attained a significant accuracy level (based on the determination coefficient (R2) value) in many cases. These CSC models have also been successfully employed to solve problems in different engineering fields (Zhang and Goh, 2013, 2016; Chen et al., 2020;Goh et al., 2020; Ghani and Kumari, 2021; Ghani et al., 2021a;Kardani et al., 2021a; Kumar et al., 2021c). Despite the good performance of ANN,SVM,ANFIS,and many of the other conventional methods, they are considered black-box models (Alavi and Gandomi, 2011; Mohammadzadeh et al., 2014; Bui et al., 2018).Also, the overfitting and local minima trapping issues limit their applicability for further analysis,which are considered as the main disadvantages of CSC algorithms (Armaghani et al., 2014; Roy and Singh, 2020; Bardhan et al., 2021a).

Therefore, to circumvent the limitations of the CSC models,recent studies have resorted to hybrid soft computing model as potential alternatives to predict the desired output(Bui et al.,2018;Bardhan et al., 2021b). An amalgamation of meta-heuristic optimisation algorithms and CSC models yields high-performance computational models, which balance the exploration and exploitation processes during the course of optimisation, and offers a bendy and powerful method for solving high-dimensional and complex problems(Mohammadzadeh et al.,2014;Bui et al.,2018).For predicting the performance of TBM, Armaghani et al. (2017)used two ANN-based hybrid models, i.e. particle swarm optimisation (PSO)-ANN (ANN optimised with PSO) and imperialism competitive algorithm(ICA)-ANN(ANN optimised with ICA),Zhang et al. (2020a) proposed a hybrid method of RVM and PSO; Zeng et al. (2021) used six hybrid PSO-ELM techniques; Zhou et al.(2021) used three hybrid SVM models, i.e. grey wolf optimisation(GWO)-SVM,whale optimisation algorithm(WOA)-SVM,and moth flame optimisation (MFO)-SVM, and found notable results.Armaghani et al. (2017) also compared the performance of PSOANN and ICA-ANN models with that of the conventional ANN,and reported that the hybrid intelligent models are superior to the simple ANN method. On the other hand, Zeng et al. (2021)concluded that the hybridisation of ELM and PSO results in higher accuracy than the single ELM model,and Zhou et al.(2021)showed that the hybrid SVM model (i.e. MFO-SVM) is more powerful than standard SVM in addressing TBM’s performance with a high accuracy level.A thorough examination of these studies reveals that most of the authors consider single CSC models and hybrid models in predicting the performance of TBM.However,the performance enhancement of the CSC models through ensemble modelling was not addressed in earlier studies.

Recently,some researchers have focused on ensemble machine learning models for predicting the desired outputs in different engineering disciplines(Basaran et al.,2019;Cai et al.,2020;Wang et al.,2020;Asteris et al.,2021;Kardani et al.,2021b,c;Zhang et al.,2021a, b). It is pertinent to say that ensemble models work on an advanced level to augment the performance of “unstable” models(Breiman,1999;Qiu et al.,2014;Chatterjee et al.,2015).It must be noted that an ensemble model constructed with two or more individual soft computing models exhibits better performance compared to a single model(Qiu et al.,2014;Chatterjee et al.,2015;Cai et al.,2020).Hence,taking these points as a reference,this study aims to implement a recently proposed approach, named hybrid ensemble technique(Asteris et al.,2021),for estimating the ROP of TBM. To this end, four common CSC techniques, i.e. MPMR, RVM,ELM, and FN, are used, analysed, and discussed. Among them, the former two models are categorised as regression-based models,while ELM and FN are NN-based models. However, the working principles of ELM and FN are different. As stated above, these CSC algorithms have also been successfully utilised in many engineering disciplines.Nonetheless,the application of MPMR,ELM,and FN in estimating TBM’s performance is relatively scarce.

In the subsequent step,a hybrid ensemble model(HENSM)was used for aggregating the outputs of the aforementioned CSC models via an ANN model. In particular, the ANN was utilised to combine the CSC model outputs, while the HENSM was formulated to predict the ROP of the TBM.The prediction performance of the HENSM and CSC models was investigated with respect to a number of performance indices followed by uncertainty analysis and statistical testing. The empirical results confirmed the applicability and benefits of the introduced HENSM in TBM performance estimation in a rock environment.

The remainder of this study is organised into the following sections, including the introduction section presented above. Section 2 presents the theoretical details of the utilised models,including the approach of HENSM, while Section 3 explains the details of the collected dataset, descriptive statistics, statistical analysis, and performance parameters. This is followed by a detailed discussion on the performance of the employed models in Section 4. Finally, in Section 5, the concluding remarks are provided.

2. Theoretical background of the employed models

This section discusses the theoretical background of the employed models in predicting the ROP of TBM. It begins with presenting the details of CSC models,i.e.MPMR,RVM,ELM,and FN employed. Finally, the details of the ANN and the methodological development of the proposed HENSM are presented.

2.1. Minimax probability machine regression (MPMR)

MPMR is a regression framework (Strohmann and Grudic,2003), which is established based on the minimax probability machine classification algorithm(MPMCA).A regression model can be formulated as maximising the minimum probability of predicted values where the prediction of the output lies within some bound of the true function. The regression formulation closely shadows the classification formulation (Marshall and Olkin, 1960; Huang et al., 2006), which can be given as

Based on kernel formulation, the MPMR model is depicted below, which shows the relationship between the independent variable (x) and response variable (y):

where K(xi，x) represents a kernel function, b and βiare the outputs,and n is the total number of observations of the input variable.In general,the radial basis function(RBF),i.e.K(xi，x) = exp[ - (xi，x)(xi，x)T/(2σ2)],is used,where σ is the width of the RBF.To use the above MPMR formulation for regression,original data with d inputs are used to create two classes of points (gi， hi) as follows:

The minimum probability machine classification boundary between points giand hiis shown below:

The above classification boundary between giand hicorresponds to a regression surface,termed as MPMR model.Eq.(4)can formulate using any binary classification algorithm, and it is not limited to MPMCA. Precisely, if a binary classifier derived from Eq.(1)has the above form where c=-1 for the first class and c=1 for the second class, then a crossing point needs to be found first. For the crossing point,the classifier separates these classes for a certain input (X = (x1， x2， …??， xd)) which is equivalent to finding the output of the regression model for a given input combination. For detailed methodology, the work of Strohmann and Grudic (2003)can be referred to.

2.2. Relevance vector machine (RVM)

RVM, introduced by Tipping (1999), is a supervised machine learning technique that uses Bayesian interpretation to find regression solutions and probabilistic classification.It has the same functional formulation as SVM, but provides probabilistic interpretation. In supervised machine learning, inputs (x1, x2, …, xn)come allied with their corresponding outputs (t1， t2， …??， tn). The outputs might be real values for regression and class labels for classifications.For a given training dataset,machine learns a model based on input and output values.SVM is also a supervised machine learning technique that makes predictions based on a function given by

where wnis the model weight.RVM uses training data for building regression model in the form of

where εiis the noise element of the predicted vector with the mean of 0 and variance of σ2,which is denoted as N(0，σ2).In the case of a linear model, the prediction function is a linear combination of basis functions αi(x) and depicted below:

The above equation can be rewritten in the form of y=Φw+ ε,where Φ is a matrix of N×M dimensions.Its ith column is formed with the values of the basis function αi(x)at all training points.The noise vector is ε=[ε1， ε2， …??， εN]. For a given dataset of independent and dependent pairs {xn， yn｝Ni=1, standard formulation is normally followed and assumes that p(y|x) is the multivariate Gaussian distribution N(y||y(x)，σ2).As defined in Eq.(6),the mean of this distribution is modelled by y(x)for a given x.The likelihood of the dataset can be given by

where t = (t1， t2， …??， tN), w = (w0， w1， …??， wN)， and Φ′is a design matrix of N × (N + 1) dimensions with Φ′nm= K(xn，xm-1)and Φn1′ = 1. Maximum likelihood estimation of w and σ2often results in overfitting. As suggested by Tipping (1999), the imposition of some prior constraints on w by adding a ‘complexity’ penalty term to the likelihood or error function can overcome the overfitting. This adjusts the simplification ability of the learning process.Generally,new parameters(high level)are used to compel an unambiguous zero-mean Gaussian prior probability distribution over the weights. Using Bayes’ rule, the posterior distribution of weights is given by

2.3. Extreme learning machine (ELM)

ELM,introduced by Huang et al.(2006),is a single hidden layer feed-forward NN used for solving regression, clustering, and feature learning problems. During the course of the training and predicting stages,weights and bias are assigned randomly between inputs and hidden nodes which remain constant throughout the operation. On the other hand, weights between the hidden and output layers can be learned quickly.In addition,the computational cost of ELM is on the lower side compared to another backpropagation algorithm network. A typical structure of an ELM is shown in Fig.1.

Fig.1. A typical network of an extreme learning machine.

where bk(k = 1， 2， …， q)is the bias of kth hidden neuron,wk?Rpis the weight vector of kth hidden neuron, and f(·) denotes the activation function(sigmoidal function).It is worth noting that the bias(bk)and weight(wk)vectors are produced in a random manner from which the hidden layer output matrix(H)can be formed.The weight matrix is then calculated using the‘Moore-Penrose pseudo inverse’ method, which is given by

where Y is the predicted class label. Detailed working principle of ELM can be referred to the literature (Zeng et al., 2021).

2.4. Functional network (FN)

FN, introduced by Castillo et al. (2000), is an extension of NNs.The main benefit of NNs is to learn from records with the help of various learning methods such as parametric and structural learning methods. In NN, the error function is minimised for the learning process. The main disadvantage is that the solution obtained from NN cannot be guaranteed as the best one due to the limited ability to develop the theory and uncertainties involved in the problem. While solving a complicated problem, such limitations can be overcome by using FN. FN is a more comprehensive version of NNs, and it combines domain knowledge and data. FN uses structural learning and parametric learning to develop a model. In structural learning, functional equations are used for simplification, whereas in parametric learning, neuron functions are estimated based on the combination of shape functions.

Note that FN is a combination of storing units, layers of the computing unit, and a set of links. Storing unit consists of input storing units,intermediate layer units,and output units.Neurons of the computing unit assess the set of input values approaching from the following layer. Functions that are not arbitrary but can be determined by the network structure are termed as links. Besides data, other informational and functional properties such as invariance, commutativity, and associativity are also considered during the construction of networks.In FN,weights of the functions are merged and learned.In some cases,the learning method leads to a global minimum quite easily. A graphical representation of an FN is shown in Fig. 2. In addition, the steps of FN working are depicted in Fig. 3.

The learning method of FN comprises neural functions obtained from a set of dependent and independent variables, and the learning process is based on minimising the Euclidean norm of the error function (Rajasekaran, 2004):

where θ is the shape function, and it can be associated with different functions such as algebraic expressions, trigonometric functions, and exponential functions. These associative optimisation functions lead the formulation to a system of linear or nonlinear equations.For a system with three inputs(x1，x2，x3)and one output (x4), the association of the FN can be constructed, as shown in Fig.4,the mathematical formulation of which is given by

where mscan be of any order; and θsican be polynomial of any degree, trigonometric, trigonometric polynomial, exponential, or any other mathematical functions,termed as shape functions.For a polynomial expression, the function can be expressed as

Fig. 2. Architecture of functional network.

Fig. 3. Flowchart showing the working steps of functional network.

Fig. 4. Associativity of functional network.

where xk0is the initial value. For the uniqueness of the above solution,the following equation must be satisfied:

To achieve more generalised functions, it is required to guarantee the individuality of the illustration of the network. Sometimes different neuron functions lead to a same output for a similar input which leads to an ill-conditioned estimation. To solve these problems,a generalised method is used for all types of FNs(Castillo and Gutiérrez,1998).

2.5. Artificial neural network (ANN)

ANNs,inspired by the biological neural system,consist of many single units and weighted artificial neurons (Asteris and Mokos,2020). An ANN is also recognised as a processing element (PE)since it processes information. Each PE contains one or more weighted inputs,an output,and a transfer function.Basically,a PE is an equation that balances the inputs and output. An ANN is also known as a connectionist model since the connection weights stand for the system’s memory.

An ANN involves three main sets of layers with neurons. The first and last layers of an ANN are known as the input and output layers, respectively. They contain many neurons as the input and output variables. Also, hidden layers exist between these two layers.Signals are transmitted through the input layer.The hidden layers serve as the computation engine, and predictions are performed within the output layer.The weights and biases are crucial parameters of an ANN.The weights represent the interconnections between a layer’s neurons, while the biases determine the degree of freedom.Each of the nodes,except for the input nodes,utilises a nonlinear activation function to produce the output. The output is then utilised as the input for the next node.This process continues until a proper solution to a given problem is found.To calculate the error by comparing the actual result (i.e. the objective of the problem)to the prediction(i.e.the outcome of the network),backpropagation is utilised.Then,the error propagates to a layer back at a time throughout the ANN structure,changing the weights on the basis of their error contributions.

2.6. Implementation procedure of hybrid ensemble model(HENSM)

Ensemble learning is defined as a machine learning process that is used to achieve superior results through the strategical combination of several learning algorithms (Qiu et al., 2014; Chatterjee et al., 2015), i.e. the principle of ensemble modelling. An ensemble model merges several standalone/individual models to deliver better prediction accuracy (Dietterich, 2000). Specifically,ensemble models have three main advantages: (i) the ability to decrease the risk of choosing the worst-performance model by the aggregation of multiple candidate models;(ii)efficiency in dealing with locally optimal solutions; and (iii) a representational advantage (in numerous cases, one model cannot represent a true function by a single hypothesis or formulation) (Qiu et al., 2014). It is worth noting that, for regression and time series forecasting problems, different machine learning models that use identical data produce different prediction results. Thus, all the outputs created by several machine learning models may be combined,providing a more robust prediction compared with those by separate machine learning models. An analysis of the relation between the actual and predicted outputs makes it promising to create a more effective model via different ensemble methods,such as randomness injection, voting, boosting, stacking, and bagging(Bar et al., 2013; Qiu et al., 2014).

Nonetheless, the present study uses an approach proposed by Asteria et al. (2021) instead of standalone ensemble techniques mentioned above. The presented HENSM is made up of four individual machine learning models (i.e. MPMR, RVM, ELM, and FN)and an ANN. In particular, the ANN was employed for aggregating the outputs of the four individual machine learning models to formulate the HENSM. The process details are as follows:

Fig. 5. Schematic diagram of the proposed hybrid ensemble model.

(1) The standalone machine learning models are trained via input data X and output data Y.

(2) The predicted values are generated.

(3) The predicted values of the individual machine learning models are amalgamated.

(4) New input data Xneware formed,and an ANN is trained using the real predicted values Y.

(5) An HENSM is constructed.

Fig. 5 shows the overall diagram of the proposed HENSM method used in this study for predicting the ROP of TBM.

3. Data description and computational analysis

3.1. Descriptive details of the collected dataset

A sum of 185 datasets was collected from the literature (Javad and Narges, 2010) and was used for predicting the ROP of TBM.The collected dataset comprised 150 data of Queens water tunnel project of USA,8 data of Karaj-Tehran tunnel project of Iran,and 27 data of Gilgel Gibe II hydroelectric project of Ethiopia. The main dataset consists of three input parameters,i.e.uniaxial compressive strength (UCS), rock quality designation (RQD), and distance between planes of weakness (DPW). These three parameters were used to predict the ROP (in m/h), i.e. the output variable. Table 1 represents descriptive statistics of the collected dataset with the minimum,the maximum,the mean,and so on,of each variable.In addition,Figs.6 and 7 plot the correlation matrix(based on Pearson correlation coefficient)and frequency histograms,respectively.The frequency histograms are presented in the form of a comparative histogram between the input and predicted variables.This diagram is useful to compare the input variables against the output. However,for a fair comparison,all the parameters were scaled in a predefined range between 0 and 1. From Table 1, it is noted that the UCS of rocks lies in the range of 30-199.7 MPa, RQD has the minimum and maximum values of 40.6% and 99.88%, respectively,while the DPW locates between 0.05 m and 2 m. The output variable, ROP, is distributed in the range of 1.27-14.43 m/h.

3.2. Statistical analysis

From the descriptive analysis presented above, it was understood that the present database comprises a wide range of data,and hence statistical analysis was performed to assess the linear association between the input parameters and ROP,the details of which are furnished in Fig.8a-c.The ROP values were correlated with the UCS, RQD, and DPW and presented in terms of determination coefficient (R2) value. As can be seen, the values of R2scatter in therange of 0.0944-0.2133, which in turn indicates that the parameters are marginally correlated with ROP. To visualise the nature of the input parameters,the collected dataset is presented in Fig.8a-c through scatter plots of each input parameter.dimensional effects. Therefore, before the development of the models, the main dataset was normalised in the range of 0-1.Subsequently,the normalised dataset was partitioned into training and testing subsets. For this purpose, 70% of the main dataset (i.e.130 samples) was selected at random for the training subset, and the left-over dataset, i.e. 30%of the main dataset(i.e. 55 samples),was used as the testing subset. The entire steps of developing the predictive models of ROP of TBM can be described as: (i) normalisation of the main dataset, (ii) selection of training and testing subsets, (iii) construction of models using training subset, and (iv)validation of the developed models using the testing subset.

Table 1 Descriptive statistics of UCS, RQD, DPW,and ROP.

Fig. 6. Correlation matrix.

Consequently, to evaluate the models’ performance, several widely used statistical indices, i.e. determination coefficient (R2),performance index (PI), variance account for (VAF), Willmott’s index of agreement(WI),mean absolute error(MAE),mean bias error(MBE), root mean square error (RMSE), and weighted mean absolute percentage error (WMAPE), were determined (Kardani et al.,2020, 2021a, d; Bardhan et al., 2021a, b; Ghani et al., 2021a,b;Kumar et al.,2021c;Ray et al.,2021).The mathematical expressions of the above-mentioned indices are given in the following equations:

3.3. Artificial intelligence-based analysis

In the soft computing field, data normalisation is a preprocessing task, which is performed to cancel out the multi-

Fig. 7. Frequency histograms.

Fig. 8. Liner relationship between the inputs and ROP.

Fig. 9. Illustration of RVM weights corresponding to the training dataset.

where yiandirepresent the actual and estimated ith values,respectively; n is the number of samples for the dataset under consideration;ymeanis the mean of the actual values;and adj?R2is the adjusted R2. Note that for an ideal model, the values of these indices should be equal to their ideal values, as given in Table 2.

Furthermore, score analysis was conducted to evaluate the developed models’ overall performance (Zorlu et al.,2008). It may be noted that score analysis is an easy and effective method for comparing the models’performances.In this technique,the model with the least value for each performance parameter is assigned a score of 1, and the model with the best value for the same performance parameter is assigned a score of m (the total number of considered models) for training and testing results separately. In the present work, m = 5, i.e. the considered models are MPMR,RVM, ELM, FN, and HENSM. The subsequent step involves calculating the total score by adding up their separate scores. The ultimate model score was computed by adding the total scores belonging to the training and testing phases.

4. Results and discussion

4.1. Realisations of conventional soft computing models

As stated earlier, a sum of 185 datasets was collected and partitioned into training and testing subsets, of which the training subset was employed to construct the models, while the testing subset was utilised to validate the developed models.Subsequently,different indices were determined to measure the prediction capabilities of the developed models, the details of which are presented and discussed below. Note that the same training dataset was used in all cases. However, before discussing the model’s performance, the parametric configurations of each model are presented below.

In MPMR modelling, two hyper-parameters, i.e.σ (width of the RBF) and ε (noise component of the measurement) play an important role, and hence the values of these parameters wererigorously checked within a predefined range of 0.5-5 for σ and 0.001-1 for ε,respectively.Following the trial-and-error approach,the best performance was obtained when the values of σ and ε were set to 1.1 and 0.003,respectively,which are the designed values of σ and ε for the developed MPMR model.

Table 2 Ideal values of performance indices.

For modelling the RVM, the widely used Gaussian kernel was used(Samui et al.,2011)for the same training dataset.The value of σ was selected through trial-and-error runs and the best performance was obtained at σ = 0.4. Therefore, the denominator term(2σ2)of the RBF,i.e.exp[ - (xi，x)(xi，x)T/(2σ2)],can be calculated as 2 × 0.42, i.e. 0.32. The details of the weight vector against the training dataset are shown in Fig. 9. The following equation represents the developed RVM model for calculating the ROP of TBM:

where xiis the training dataset of 130 observations, and x is the input for which ROP is to be determined.

In contrast, to determine the ELM structure, the number of hidden neurons varying from 2 to 30 was examined. In this study,sigmoid activation function was used to select the best possible ELM model. Through trial-and-error runs, the most appropriate value of hidden neurons was obtained as 12.

The fourth CSC model, i.e. FN, was also constructed using the same training dataset as those used in MPMR, RVM, and ELM modelling.A trade-off was made in the present study,and FN with degree 3 and polynomial basis function was adopted.The equation for the prediction of the ROP of TBM was obtained as follows:

where m′is the degree of variable.

4.2. Realisation of hybrid ensemble model (HENSM)

This sub-section discusses the realisation of the HENSM. Note that a model that achieves better accuracy in the testing phase is usually considered to be a robust model and accepted with higher conviction. Therefore, one may create an even more robust model by improving the testing phase performance. With this consideration,the superiority of the ensemble method was leveraged in the work to construct an HENSM.This will be more advantageous with respect to the variational inaccuracies of the utilised models.

Fig.10. Utilisation of standalone models for the development of hybrid ensemble model. TR- Training; TS - Testing.

It is worth mentioning that in many cases, computational models produce different results with different accuracies when there are insufficient samples. Hence, ensemble techniques can decrease the risk of choosing the less accurate model via the aggregation of multiple individual models. The present study attempts hybrid ensemble modelling via an HENSM combining four CSC models that are used to estimate the ROP of TBM.As mentioned in Section 2.6,the ANN was used for constructing the HENSM.This model combines the outputs from all four individual models as inputs and the real output, i.e. actual ROP of TBM. In Fig. 10, the process of implementing the HENSM is presented,from which the amalgamation procedure of CSC models can be visualised. It must be noted that a trial-and-error approach was followed to finalise the structure of ANN in HENSM modelling. In this case, the Levenberg-Marquardt back-propagation function was selected as the training function. The ultimate ANN structure is composed of three input neurons, two hidden neurons, and one output neuron.The activation function used in input to hidden and output layers are tan-sigmoid and purelin,respectively.Also,a total of 10%of the combined training dataset was utilised as the validation dataset to stop overfitting.

4.3. Performance indices and evaluation

A comprehensive assessment of the CSC and HENSM models based on multiple performance indices is presented in this subsection, the details of which are listed in Tables 3 and 4 for the training and testing datasets,respectively.Furthermore,the scatter plots between the actual and predicted outputs are presented in Figs.11 and 12.As can be seen,overall,all CSC models were able to predict the ROP with significant accuracy that can be evident by high R2value (ranging between 0.7214 for FN and 0.911 for ELM)and low RMSE(ranging between 0.0286 for ELM and 0.0505 for FN)in the training phase.For the testing dataset,the values of R2were obtained in the range of 0.7767 (for MPMR) to 0.8518 (for RVM),while the values of RMSE were obtained in the range of 0.0453(for RVM) to 0.0532 (for FN). Based on the score analysis, the best standalone model is ELM (total score of 30) in the training phase and RVM (total score of 30) in the testing phase. Contrarily, the constructed HENSM attained the most accurate prediction with the highest R2(0.9126 for training and 0.865 for testing dataset) and lowest RMSE value (0.0283 for training and 0.0418 for testing dataset) in both cases. The outcomes of score analysis also show that the constructed HENSM outperformed the best performing CSC model and attained a final score of 36 and 38 in the training and testing phases, respectively. The results of score analysis are presented in Fig. 13 (for training results) and Fig. 14 (for testing results) in the form of stacked bar plots.The comparisons of score analyses presented in Table 5 indicate that RVM and ELM are thetop two standalone models (attained overall scores of 54 and 52,respectively); however, when the outputs of each model are combined with the ANN,they exhibit better performance than the best standalone model in predicting the ROP of TBM.

Table 3 Outcomes of the training dataset.

Table 4 Outcomes of the testing dataset.

Fig.11. Scatter plots between actual and estimated ROP values for the training (TR) dataset.

In addition, the objective function (OBJ) criterion was used to assess the predictive performance of all the employed models.The details of this parameter can be found in the literature (Gandomi et al., 2013; Samadi et al., 2020; Raja et al., 2021; Raja and Shukla,2021b). However, the expression presented by Samadi et al.(2020) was used. It should be noted that the best model has a minimum OBJ value(ideal value=0).The values of OBJ for MPMR,RVM,ELM,FN,and HENSM were calculated as 0.0366,0.03,0.0307,0.0455,and 0.0284,respectively.These results indicate that HENSM attained the most desired prediction with the lowest OBJ value,and hence the proposed hybrid model is most accurate from this viewpoint.

To better demonstrate the values of performance indices, a recently proposed heat map-shaped graphical assessment, called accuracy matrix,is tested to visualise the model efficiency(Kardani et al.,2021a).This matrix displays multiple statistical parameters to measure the model’s predictive performance for the training and testing datasets. Fig. 15 displays the accuracy matrix for the performance indices determined in this study.It indicates the accuracy of the performance parameters(in percentage)by comparing them with their corresponding ideal values.For example,the ideal value of MAE is 0 and, in the present work, the value of MAE for the training subset was determined as 0.0227 for the MPMR model(see Table 3). Thus, it can be estimated that the MPMR model attained 97.73%((1-0.0227)×100%)accuracy in terms of MAE.On the other hand,the values of R2and PI were obtained as 0.865 and 1.6715 in the testing phase,respectively,for the HENSM(see Table 4),which shows that HENSM attained 86.5% (0.865/1 ×100%) and 83.575%(1.6715/2 × 100%) accuracy in terms of R2and PI, respectively.Similar procedure was followed for the other parameters as well.However, it may be noted that the parameters such as VAF, which are determined in percentage terms, should be converted in their decimal form before implementing the above procedure. The supremacy of the developed HENSM based on uncertainty analysis(UA) and statistical testing is presented in the following subsection.

Fig.12. Scatter plots between actual and estimated ROP values for the testing (TS) dataset.

4.4. Uncertainty analysis and statistical testing

This sub-section presents the quantitative evaluation of the developed models through UA in the prediction of the ROP of TBM.For this purpose, the testing dataset, which contains 55 real-time observational data points, was used. Thus, it can be useful to logically compare predictive outputs for evaluating the reliability of the predictive models (Gholami et al., 2019). To carry out the UA, the absolute error (εi) between observational (yi) and predicted (i)values,mean of absolute error(MOE),and standard deviation(SD)of error were determined using the following expressions:

where N is the number of observations(here N=55).Initially,MOE and SD were determined for the testing datasets using the above equations. Subsequently, the width of confidence bound (WCB)value was determined by calculating the margin of error(ME)at a 95% confidence interval. The standard error (SE), lower and upper bounds (LB and UB) were calculated as follows:

Note that in UA, WCB is an error range where about 95%of the entire data are located. Table 6 provides UA details, including the MOE,SD,SE,ME,LB,UB,and WCB.These indices help to assess the efficiency of the models. The smaller the value of WCB, the better the model certainty, meaning that the model will suffer from less error, and it will predict more precise expected output. Based on the values of WCB, all the five models were ranked. It is seen that the proposed HENSM attained the lowest WCB value and secured the first rank, rendering it the most reliable predictive model. The RVM model scored the second rank, which is the top performing model among the CSC models, followed by ELM, FN, and MPMR.However,the MPMR model exhibits the highest WCB value and the lowest uncertainty.Fig.16a-c provides UA graphical representation as a bar chart displaying MOE, WCB, and ME values for a better comparison.

Furthermore, a one-tailed t-test was executed to estimate the significant difference of the proposed HENSM’s performance to the CSC models in predicting the ROP of TBM in a rock environment.This statistical test was executed on the MAE values with the hypothesised mean difference HMD = 0. At HMD = 0 and confidence interval of 95%(i.e.α=0.05),the hypotheses of one-tailed ttest are H0: MAEHENSM- MAECSCmodels= 0 and HA: MAEHENSM-MAECSC-models<0,where H0and HArepresent null hypothesis and alternate hypothesis, respectively. The outcomes are presented in Table 7.Note that the rejection(failed to accept)of H0(i.e.t-stat

4.5. Sensitivity analysis

Fig.13. Bar plot of score analysis (training results).

Fig.14. Bar plot of score analysis (testing results).

Table 5 Overall score analysis for the developed models.

4.6. Discussion

From the experimental results presented in the above subsections, it is observed that RVM and ELM are the top two conventional models. However, detailed scrutiny reveals that RVM attained the most accurate prediction in the testing phase only.Although overall analysis shows that RVM outperformed other CSC models with R2=0.8518 and RMSE=0.0453,a lower performance in the testing phase (R2= 0.9099 and RMSE = 0.0287) against the performance of ELM (R2= 0.911 and RMSE = 0.0286) encourages the authors to analyse the prediction capability in another way.Therefore, a hybrid ensemble method was implemented to augment the CSC models’ performance, especially for the testing dataset. The performance of the proposed HENSM was then assessed and compared for predicting the ROP of TBM in a rock environment.Overall,the proposed HENSM attained the maximum final score of 74 (36 (training) + 38 (testing)) by far, followed by RVM(54),ELM(52),MPMR(35),and FN(25).These results indicate that RVM and ELM are the top two models among the 4 conventional models (i.e. MPMR, RVM, ELM, and FN) employed in this study. However, the FN is the worst performing model with the lowest final score of 25.

Analogous to CSC models,the proposed HENSM was developed utilising the outputs of the CSC models in a MATLAB (2015a version) environment. The computation costs were recorded as 0.311139 s for MPMR, 0.823362 s for RVM, 0.325061 s for ELM,7.864904 s for FN,and 8.1252 s for HENSM.All the programmes run in MATLAB with i3-8130U CPU @ 2.20 GHz,12 GB RAM, which indicates that the constructed HENSM attained higher prediction accuracy with very low computing cost. This is one of the main advantages of the proposed approach due to the fact that this approach can be implemented quickly.Also,the present technique can reduce the risk of picking the least performing model by combining individual models.

5. Concluding remarks

Based on the data mining from the literature,185 data records of ROP of TBM were collected, which cover three rock properties, i.e.UCS, RQD, and DPW. These parameters were used as the input parameters to predict the ROP of TBM in a rock environment. The collected dataset was used to rigorously train and validate the employed models to enable accurate predictions of ROP of TBM.Four standalone models, i.e. MPMR, RVM, ELM, and FN, and an HENSM were employed in this study. The accuracies of these models in predicting the ROP of TBM were analysed and compared from different aspects.Subsequently,score analysis was performed to select the best-performing model.The UA and statistical testing were also performed to estimate the reliability of the proposed HENSM in predicting the ROP of TBM. The experimental results reveal that the constructed HENSM attained higher performance compared to the standalone models and can be used effectively to estimate the ROP of TBM in the design phase of engineering projects.

To sum up, a hybrid ensemble approach is successfully applied in this study to predict the ROP of TBM in a rock environment. It appears to be an accurate yet computationally efficient method due to the following reasons: (i) easy implementation procedure; (ii)higher prediction accuracy; (iii) very low computational cost; (iv)produced greater performance compared to CSC models, and (v)being representational. Experimental results exhibit that the HENSM outperformed other models at all levels (in terms of performance indices,score analysis,UA,and statistical testing),which in turn indicates that the proposed model has high generalisation capability in handling overfitting related issues of CSC techniques.Despite these advantages,the proposed hybrid ensemble technique can be improved in the future as follows:(i)incorporation of large datasets to predict the desired output(s) more accurately; (ii) a detailed assessment of results of the testing dataset through crossvalidation of different conventional machine learning models; (iii)adding expert endorsed data to cover much broader diversities;(iv)extending the application of the proposed approach in different fields; and (v) amalgamation of meta-heuristic optimisation algorithms and in-depth assessment of the proposed model over the hybrid ANN and meta-heuristics models. Besides, more input parameters should be considered to predict the ROP of TBM more accurately,particularly in a rock environment.However,this study presents the application of hybrid ensemble techniques for predicting the ROP of TBM in a rock environment for the first time.

Fig.15. Accuracy matrix. TR - Training; TE - Testing.

Table 6 Results of uncertainty analysis.

Table 7 Results of one tailed t-test.

Fig.16. Bar plots of UA showing: (a) MOE, (b) WCB, and (c) ME values.

Fig.17. Bar plots of sensitivity analysis for the employed models.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.