Application of Soft Computing, Statistical and Multi-Criteria Decision-Making Methods to Develop a Predictive Equation for Prediction of Flyrock Distance in Open-Pit Mining

: Blasting operations in open-pit mines generally have various management strategies relating to ﬂyrock. There are empirical models for calculating the ﬂyrock distance, but due to the complexity and uncertainty of rock properties and their interactions with blasting properties, there are still no models that can predict the ﬂyrock distance that may be applicable across mining operations in general. In this regard, the Jajarm bauxite mine complex was used as a case study. The purpose of this study was to develop and evaluate different methods that can predict ﬂyrock distance. For this purpose, soft computing models were developed using generalized regression neural network (GRNN), gene expression programming (GEP) and genetic-algorithm-based GRNN (GA-GRNN) methods. To obtain statistical models, multivariable regression was applied in the form of linear and nonlinear equations. A ﬂyrock index was introduced using a classiﬁcation system developed by incorporating fuzzy decision-making trial and evaluation methods (fuzzy DEMATEL). In order to achieve this goal, the data of 118 blasts in eight mines of the Jajarm bauxite complex were collected and used. Following this, four performance benchmarks were applied: the coefﬁcient of determination (R 2 ), variance accounted for (VAF), root-mean-square error (RMSE) and mean absolute percentage error (MAPE). The performance of the models was evaluated, and they were compared with each other as well as with the most common previous empirical models. The obtained results indicate that the GA-GRNN model has a higher performance in predicting the ﬂyrock distance in actual cases compared to the other models. At ﬁrst, data on factors that were the main cause of ﬂyrock (and had a direct impact on it) were collected and classiﬁed from different blasts. Then, using the collected data, 19 different combinations were established, which can be used to provide the appropriate predictive equation. The purpose of this work is to more accurately predict ﬂyrock and prevent heavy damage to buildings and mining machines across the mining complex.

. Types of flyrock in open-pit mines (blue indicates charge length of the hole) [12].
The three mechanisms are as follows [18,20]: Face bursting: This phenomenon takes place when the surface of the bench is not smooth or the blast hole is filled with explosives close to the earth's weak structures and the max. burden is lower than normal.
Rifling: This type of phenomenon, which is one of the most dangerous types, occurs when the stemming of the blast hole is not properly carried out or the materials used do not have adequate efficiency in trapping explosive gases.
Cratering: Various factors contribute to the occurrence of this phenomenon, such as: • The presence of pebbles around blast holes as a result of drilling; • The lack of consideration of the delay sequence in blast holes; • The presence of weak layers in the superficial layers of the ground. Table 1. Calculating the maximum distance covered by flyrock through an explosion [18,20]. In the equations in Table 1, is the drill hole angle, LMax is the maximum flyrock distance, m is the mass of the explosive in each hole (kg/m), B is the burden (m), St is the stemming height (m), g is the gravitational constant (9.81 m/s 2 ) and k is a constant parameter according to the ground's condition.  [12]. Table 1. Calculating the maximum distance covered by flyrock through an explosion [18,20].

Mechanism of Flyrock Prediction Model
Rifling sin 2θ Cratering Face bursting The three mechanisms are as follows [18,20]: Face bursting: This phenomenon takes place when the surface of the bench is not smooth or the blast hole is filled with explosives close to the earth's weak structures and the max. burden is lower than normal.
Rifling: This type of phenomenon, which is one of the most dangerous types, occurs when the stemming of the blast hole is not properly carried out or the materials used do not have adequate efficiency in trapping explosive gases.
Cratering: Various factors contribute to the occurrence of this phenomenon, such as: • The presence of pebbles around blast holes as a result of drilling; • The lack of consideration of the delay sequence in blast holes; • The presence of weak layers in the superficial layers of the ground.
In the equations in Table 1, θ is the drill hole angle, L Max is the maximum flyrock distance, m is the mass of the explosive in each hole (kg/m), B is the burden (m), S t is the stemming height (m), g is the gravitational constant (9.81 m/s 2 ) and k is a constant parameter according to the ground's condition.
Finally, two general results, concerning the properties of the rock mass and the flyrock distance, can be obtained by using the evaluation of different mechanisms of flyrock and its factors [5,6,21,22]: Choosing improper explosives by disregarding the geological conditions and properties of the rock mass; 2.
Using an improper drilling pattern and having inaccuracies in execution. Table 2 showcases some of the most significant works of a wide range of research that has been conducted on the flyrock phenomenon. Table 2. Summary of research on the flyrock phenomenon.

Authors Years Findings
Ladegaard-Pedersen and Holmberg [23] 1973 The charging geometry affects cratering flyrock more than the other two types.
Bajpayee et al. [28] 2002 A safe zone can be provided to prevent casualties and equipment damage.
Verakis and Lobb [29] 2003 The lack of consideration of geological conditions when choosing blasting patterns and explosives, improper stemming and charging and an unfit blasting sequence (blast hole's delay) are influencing factors of flyrock. Workman and Calder [30] 1994 Kecojevic and Radomsky [4] 2005 Flyrock depends on factors such as geological structure, improper blasting pattern, the erroneous choice of burden, explosive aggregation, poor stemming and inaccuracy in choosing the correct blasting delay.
Monjezi et al. [31] 2007 The distance of flyrock in a blast was controlled using the Topsis method.
Aghajani-Bazzazi et al. [9] 2009 Controllable parameters were used to develop an empirical model to predict flyrock by using the multivariable regression method.
Rezaei et al. [32] 2011 Flyrock distance predictions based on the fuzzy method (FIS and artificial neural network) and the statistical method were compared.
Monjezi et al. [33] 2011 The ANNS method and its functionality were applied to predict the flyrock distance.
Monjezi et al. [34] 2012 The genetic neural network model was used to predict flyrock and backbreak.
Amini et al. [35] 2012 The results of flyrock distance prediction calculated through the SVM (support vector machine) method were compared to those obtained using the ANN (artificial neural network) method.
Ghasemi et al. [36] 2014 The functionality of the developed ANN method was compared to that of fuzzy logic in predicting the flyrock distance.
Armaghani et al. [10] 2014 A new combination method, Bp-Ann, was used to predict the flyrock distance and decrease the error rate. This method is a combination of the PSO (particle swarm optimization) algorithm and the ANN (artificial neural network) function.
Faramarzi et al. [6] 2014 The functionality of the multivariable regression method was compared with the rock engineering system (RES) in predicting the flyrock distance.

Authors Years Findings
Saghatforoush et al. [37] 2016 Ant colony and optimization algorithms were used to predict the flyrock distance and backbreak, which eventually led to a new ACO method for minimizing flyrock distance and backbreak.
Esen [38] 2017 The flyrock distance was predicted with the aim of determining the safe zone in open-pit mines using the effects of parameters on flyrock.
Hasanipanah et al. [15] 2017 The PSO (particle swarm optimization) method was applied, and its results were compared to those of Multiple Linear Regression (MLR) in predicting the flyrock distance.
Armaghani et al. [39] 2020 Three different methods of machine learning techniques, i.e., PCR (principal component regression), SVM (support vector regression) and BN (Bayesian network), were applied to predict the flyrock distance, and the SVR method was chosen as the best prediction model; this model was also optimized with GWO (Gray Wolf Optimization) to decrease the flyrock distance.
Han et al. [12] 2020 The Random Forest Technique was used to select the effective parameters, which were employed in BN (Bayesian network technique) to predict the flyrock distance.
Hasanipanah, and Bakhshandeh Amnieh [40] 2020 Risk analysis was conducted and the flyrock distance was predicted using different kinds of artificial intelligence, and the fuzzy rock engineering system (FRES) was chosen as the best model for risk analysis and prediction in the studied mine.
Lu, Xiang, et al. [41] 2020 The best flyrock prediction model was determined by comparing the results of the extreme learning machine (ELM) and outlier-robust ELM (ORELM) methods to the ANN and multiple regression methods.
Nikafshan et al. [42] 2020 The Recurrent Fuzzy Neural Network (RFNN) and genetic algorithm (GA) were combined in order to establish a combination model (RFNN-GA method) for predicting the flyrock distance.

Authors Years Findings
Nguyen et al. [46] 2021 A data-driven model was introduced and used to predict flyrock. A combination of the whale optimization algorithm (WOA), support vector machine (SVM) and kernel functions was used. Four linear functions (L), radius basis function (RBF), polynomial (P) and hyperbolic tangent (HT) were used for embedding in the SVM model. Then, the WOA model was applied to optimize kernel-based SVM models. Additionally, a variety of models based on conventional data were used to predict the flyrock distance. The results showed that the WOA-SVM-RBF model had the highest accuracy in predicting the flyrock distance. Finally, the WOA-SVM model was proposed as a data-driven model for estimating fly rock with high reliability in mining.
Shakeri et al. [47] 2022 In the Sungun copper mine, Iran, the method of linear multivariate regression (LMR), imperialist competitive algorithm (ICA), adaptive neuro-fuzzy inference system (ANFIS) and artificial neural network (ANN) methods were used to predict flyrock. According to the results obtained from these methods, the authors chose Levenberg-Marquardt as the learning algorithm, log-sigmoid (logsig) as the transfer function and ANN as the optimal network. It can also be concluded that the ICA technique is more accurate in predicting the flyrock distance than LMR and ANFIS models. Finally, the sensitivity analysis revealed that the powder factor and blast hole diameters are very important in flyrock distance.
Hudaverdi [48] 2022 The variable reduction method was applied to predict flyrock. For this purpose, the dominant parameters in flyrock were selected by a multivariate statistical method. Two parallel ANFIS models were then developed. Using the results of stepwise regression, the first model was created. The second ANFIS model was then obtained based on the results obtained from the factor analysis of the model. Alternative accuracy criteria were also investigated to evaluate the prediction performance of the presented model. The results showed that standardized errors, normalized errors and Nash-Sutcliffe Efficiency were very useful for model validation. Finally, by analyzing the pre-statistical method of reducing variables, the performance of the predictive model can be increased.
Hosseini et al. [49] 2022 An artificial neural network and the fuzzy cognitive map (FCM) were integrated with z-number reliability information to predict the flyrock distance. The developed model was called causality-weighted artificial neural networks based on reliability (ACWNNsR). The reliability information of the z-number was used for uncertainty elimination in the initial matrix of FCM. Additionally, the integration of nonlinear Hebbian and differential evolution algorithms was used to calculate the weights of the input neurons. The performance of the proposed ACWNNsR model was compared with a Bayesian regularized neural network and a multilayer perceptron neural network. The results showed that this comparison leads to the accurate prediction of flyrock distance estimation. Finally, using sensitivity analysis, the burden was determined as the most important factor in flyrock.
Barkhordari et al. [50] 2022 Ensemble learning approaches such as simple averaging ensemble, weighted averaging ensemble, integrated stacking model, separate stacking model and Bayesian-extreme gradient boosting were used to predict the flyrock distance, which finally led to the presentation of a separate stacking model with a bagging meta-learner that performed better than other models. In addition, the Shapley Additive Explanations (SHAP) method was used in order to reveal the relative relevance of parameters affecting flyrock distance prediction.
Numerous researchers have presented empirical equations regarding the flyrock distance in recent years. However, some of these equations are not applicable to all mines. Others have therefore modified them to match the specific mine they were studying. Some of these are shown in Table 3. Table 3. Present empirical and modified methods in previous research.

Authors Years Finding and Equation Type of Equation
Lundborg et al. [25] 1975 Empirical Gupta [51] 1980 Ls B = 155.2 × D −1.37 Empirical Pal Roy [52] 2005 Flyrock can cover a distance ranging from a few meters to 1000 m Empirical Ghasemi et al. [53] 2012 On the basis of what has already been stated, flyrock can be considered a dynamic phenomenon that occurs as a result of the interaction between the rock mass and the explosive material. Based on previous research, despite various empirical models presented by several scientists to predict flyrock distance, there is no general model (or equation) that is applicable to all mines. This is due to the numerous parameters affecting the flyrock distance with regard to geological conditions, all of which cannot be evaluated at the same time and occur parallel to one another. Therefore, these equations cannot provide a proper and comprehensive solution to calculate the flyrock distance in different mines. As such, in order to achieve this aim, a new framework for calculating the flyrock distance is presented in this study.
This study can be divided into three stages. Firstly, data on numerous blasts in the Jajarm bauxite mines were gathered to calculate the flyrock distance. In this regard, 13 parameters-12 of which were effective and independent parameters and 1 of which was a dependent parameter (flyrock distance)-were measured and collected from the mine site. The reason for choosing these data is that except for RQD/J n , the rest of the data are effective and controllable parameters of the mine and can be used to control the flyrock distance in blasting operations, which prevents damage to structures, equipment and mineral formations in the mine. Then, 19 different combinations of these parameters were established and were used to provide the appropriate prediction equation. Secondly, the flyrock distance in the mine was predicted using statistical models (multivariable regression in the form of linear and nonlinear), equations based on fuzzy decision-making trial and evaluation laboratory methods (fuzzy DEMATEL) and soft computing methods using a generalized regression neural network (GRNN), gene expression programming (GEP) and genetic algorithm-based GRNN (GA-GRNN). Eventually, a complete model to predict the flyrock distance from the blast with regard to the mine's conditions, such as its geological structure and the blast's controllable parameters, was obtained by evaluating the results of the mentioned methods.

Site Description
The Jajarm bauxite mine complex is located 19 km away from the Jajarm city in the Northern Khorasan province in Iran ( Figure 2) and spans across a significant footprint of the mountainous, arid landscape. The geological structure of the Jajarm bauxite mine is layered (Figure 3), and it is located on the east-west mountain range in the north of the Jajarm desert at a 1000 m altitude. According to Figure 3, the layer of bauxite minerals has a slope between 40 and 60 degrees. The extraction of this mineral takes place in four or five stages in order to reach the layer of bauxite minerals, with two stages of extraction to be completed (waste-shale and sandstone with coal beds-and upper kaolinite).
The Jajarm bauxite mine complex is located 19 km away from the Jajarm city in the Northern Khorasan province in Iran ( Figure 2) and spans across a significant footprint of the mountainous, arid landscape. The geological structure of the Jajarm bauxite mine is layered (Figure 3), and it is located on the east-west mountain range in the north of the Jajarm desert at a 1000 m altitude. According to Figure 3, the layer of bauxite minerals has a slope between 40 and 60 degrees. The extraction of this mineral takes place in four or five stages in order to reach the layer of bauxite minerals, with two stages of extraction to be completed (waste-shale and sandstone with coal beds-and upper kaolinite).

Site Description
The Jajarm bauxite mine complex is located 19 km away from the Jajarm city in the Northern Khorasan province in Iran ( Figure 2) and spans across a significant footprint of the mountainous, arid landscape. The geological structure of the Jajarm bauxite mine is layered (Figure 3), and it is located on the east-west mountain range in the north of the Jajarm desert at a 1000 m altitude. According to Figure 3, the layer of bauxite minerals has a slope between 40 and 60 degrees. The extraction of this mineral takes place in four or five stages in order to reach the layer of bauxite minerals, with two stages of extraction to be completed (waste-shale and sandstone with coal beds-and upper kaolinite).

Data Sets
In this study, the required data were gathered from 118 blasts in 8 mines or anomalies in order to present a new model to predict the flyrock distance in the Jajarm bauxite mine complex. From this data, flyrock was selected as dependent (output) data, and the rest of the data were used as inputs (independent) for the flyrock distance prediction. The basic descriptive statistics of this database are summarized in Table 4. The rock mechanics information of the mine was measured based on the International Society for Rock Mechanics (ISRM) standards. Additionally, based on field observations, flyrock in most blasts was the face-bursting or rifling kind. To measure the flyrock distance, the blast face bench and the space in front of it were cleared of any rocks so that they did not interrupt the measurements, and measuring was conducted using a meter between the blast face bench and the horizontal landing place of the flyrock. A sample of the drilling pattern in these blasts and the start delay (imitation sequence) in each blast are presented in Figure 4.

Data Sets
In this study, the required data were gathered from 118 blasts in 8 mines or anomalies in order to present a new model to predict the flyrock distance in the Jajarm bauxite mine complex. From this data, flyrock was selected as dependent (output) data, and the rest of the data were used as inputs (independent) for the flyrock distance prediction. The basic descriptive statistics of this database are summarized in Table 4. The rock mechanics information of the mine was measured based on the International Society for Rock Mechanics (ISRM) standards. Additionally, based on field observations, flyrock in most blasts was the face-bursting or rifling kind. To measure the flyrock distance, the blast face bench and the space in front of it were cleared of any rocks so that they did not interrupt the measurements, and measuring was conducted using a meter between the blast face bench and the horizontal landing place of the flyrock. A sample of the drilling pattern in these blasts and the start delay (imitation sequence) in each blast are presented in Figure 4.  Due to the existence of different mechanisms of flyrock in the mine and the inability to separate the value of the flyrock distance that is specific to each type, we chose the general empirical equations (Lundborg and Gupta) listed in Table 3 for the prediction of flyrock distance; these formulas consider all states of flyrock and hence are called empirical Equations (1) and (2), respectively. As shown in the diagram ( Figure 5), the flyrock distance calculated for each blast was either more or less than the value predicted through  Due to the existence of different mechanisms of flyrock in the mine and the inability to separate the value of the flyrock distance that is specific to each type, we chose the general empirical equations (Lundborg and Gupta) listed in Table 3 for the prediction of flyrock distance; these formulas consider all states of flyrock and hence are called empirical Equations (1) and (2), respectively. As shown in the diagram ( Figure 5), the flyrock distance calculated for each blast was either more or less than the value predicted through empirical equations. This variance shows that empirical equations use site-dependent parameters that can vary in each location based on its condition and geological structure. Therefore, a specific prediction model is needed for the Jajarm bauxite mine complex.

Methodology
To conduct this study, three methods were used to present a model for the prediction of flyrock in the Jajarm bauxite mine complex: statistical models (multilinear and nonlinear regression), soft computing methods (GRNN, GA-GRNN and GEP) and fuzzy DE-MATEL. The reason for this selection is that the first and second methods are conducted according to real gathered data from each blast, and the third method is based on the knowledge and experience of experts regarding the area's features (mine site). All three methods are highly effective in predicting flyrock, and the accuracy of each can be measured and used in equations. Figure 6 shows the flowchart of a predictive model using the mentioned methods.

Methodology
To conduct this study, three methods were used to present a model for the prediction of flyrock in the Jajarm bauxite mine complex: statistical models (multilinear and nonlinear regression), soft computing methods (GRNN, GA-GRNN and GEP) and fuzzy DEMATEL. The reason for this selection is that the first and second methods are conducted according to real gathered data from each blast, and the third method is based on the knowledge and experience of experts regarding the area's features (mine site). All three methods are highly effective in predicting flyrock, and the accuracy of each can be measured and used in equations. Figure 6 shows the flowchart of a predictive model using the mentioned methods.
Mining 2023, 3, FOR PEER REVIEW 9 Figure 6. Flowchart for proposed model to predict flyrock distance. Figure 6. Flowchart for proposed model to predict flyrock distance.

Regression Analysis
The regression analysis method shows the relationship between independent and dependent variables in four stages in a linear and nonlinear (Equations (1) and (2)) manner. These stages include choosing the variables, gathering data, performing pattern recognition and fitting and validating the model [55][56][57][58].
In these equations, Y is the dependent variable; X 1 , X 2 . . . X n are the independent variables; β 1 , β 2 , . . . β n are the coefficients of independent variables; β 0 is the constant of the equation; a is the width of the origin; and ε is the total random error.

Generalized Regression Neural Network (GRNN)
Donald F. Specht presented artificial regression neural networks for the first time in 1991 [59]. These neural networks process input data using radial basis functions (such as a Gaussian function (Figure 7)) and show better performance than standard neural networks, such as feed-forward neural networks. The advantage of these networks is their less time-consuming design, quick data learning and simple operation compared to other networks in predicting data. These advantages have led to the increasing use of these networks in engineering, pharmaceutical, computing and so on [60][61][62].

Regression Analysis
The regression analysis method shows the relationship between independent and dependent variables in four stages in a linear and nonlinear (Equations (1) and (2)) manner. These stages include choosing the variables, gathering data, performing pattern recognition and fitting and validating the model [55][56][57][58].
In these equations, Y is the dependent variable; X1, X2 … Xn are the independent variables; β1, β2, … βn are the coefficients of independent variables; β0 is the constant of the equation; a is the width of the origin; and ε is the total random error.

Generalized Regression Neural Network (GRNN)
Donald F. Specht presented artificial regression neural networks for the first time in 1991 [59]. These neural networks process input data using radial basis functions (such as a Gaussian function (Figure 7)) and show better performance than standard neural networks, such as feed-forward neural networks. The advantage of these networks is their less time-consuming design, quick data learning and simple operation compared to other networks in predicting data. These advantages have led to the increasing use of these networks in engineering, pharmaceutical, computing and so on [60][61][62]. These networks operate in a way ( Figure 8) such that data first enter the network's input layer and then move to the radial basis layer for processing. In this layer, data are first classified using a data clustering method. Then, using transitional functions to learn the network, the weight relationship between the variables of the input layer and output layer is determined by decreasing the slope of the function and using linear regression. As seen in Figure 7, in this method, the Gaussian function (Equation (3)) is generally used to calculate the output of the next neuron. As the distance from the center of the function increases, the equations' response tends toward zero, which is the reason for the wide and common usage of the Gaussian function in this method. As a result, this function can divide the obtained equations radially into concentric circles and put vectors, which have equal distances from the center, into the same groups and decrease the prediction error [60][61][62][63]. Following this, the learned data enter the summation neurons and are finally sent to the output layer, which operates like a linear function, and are calculated through Equation (4). In this equation, λ is a constant, bjk is the weight constant of neuron j in the radial layer and neuron k of the output layer, and yj is the radial layer neuron's output [64]. These networks operate in a way ( Figure 8) such that data first enter the network's input layer and then move to the radial basis layer for processing. In this layer, data are first classified using a data clustering method. Then, using transitional functions to learn the network, the weight relationship between the variables of the input layer and output layer is determined by decreasing the slope of the function and using linear regression. As seen in Figure 7, in this method, the Gaussian function (Equation (3)) is generally used to calculate the output of the next neuron. As the distance from the center of the function increases, the equations' response tends toward zero, which is the reason for the wide and common usage of the Gaussian function in this method. As a result, this function can divide the obtained equations radially into concentric circles and put vectors, which have equal distances from the center, into the same groups and decrease the prediction error [60][61][62][63]. Following this, the learned data enter the summation neurons and are finally sent to the output layer, which operates like a linear function, and are calculated through Equation (4). In this equation, λ is a constant, b jk is the weight constant of neuron j in the radial layer and neuron k of the output layer, and y j is the radial layer neuron's output [64].

Genetic Algorithm (GA)
The genetic algorithm (GA) was first developed by Go function optimizer that is used for complex and nonlinear one of the problematic algorithms because of its limitation rameters of the algorithm, namely, the population size, th netic operator rate. Therefore, the operation of this algo when choosing the values of parameters to enter the algo directly affect the algorithm's response and its convergen is called a chromosome that has a fixed length and solves (0-1) strings. These chromosomes are chosen randomly bas

Genetic Algorithm (GA)
The genetic algorithm (GA) was first developed by Goldberg [66]. This algorithm is a function optimizer that is used for complex and nonlinear problems. Additionally, GA is one of the problematic algorithms because of its limitations in choosing the different parameters of the algorithm, namely, the population size, the proper function and the genetic operator rate. Therefore, the operation of this algorithm requires great attention when choosing the values of parameters to enter the algorithm because these values can directly affect the algorithm's response and its convergence [67,68]. Here, each response is called a chromosome that has a fixed length and solves problems in the form of binary (0-1) strings. These chromosomes are chosen randomly based on their different characteristics and are then evaluated so that they create a population as parents. These so-called parents create offspring (new chromosomes). These children with their new characteristics are then selected and evaluated based on their function and environmental requirements and, in a repeating cycle, create a new generation based on the remaining chromosomes (children). This cycle or genetic operator consists of three stages, called crossover, mutation and selection. In the selection step, chromosomes are selected using various transitional functions, such as the roulette wheel method. Then, in the next step, the chromosomes are joined together in pairs of fixed length. Finally, mutation plays a key role in this algorithm, as it selects chromosomes to produce a combination of chromosomes with minimal error. The cycle is continued until the best generation based on its functionality is selected. The functionality of the selected generation (the best generation) in this algorithm depends on the population of the chromosomes (usually between 50 and 100) [11,[68][69][70][71]. The aforementioned steps are schematically shown in Figure 9.

Gene Expression Programming (GEP)
Gene expression programming (GEP) is an evolutionary algorithm first introduced by Feriera in 1999 [72]. This method is a genotype-phenotype genetic algorithm with high accuracy due to its tree-like structure, and it has made a great difference in comparison to genetic programming. Additionally, it has GA's simplicity and GEP's capability; consequently, it compensates for the shortcomings of these two methods [73,74]. Chromosome and tree structures are two entities of GEP algorithms in which chromosomes contain several genes. Each gene has two parts: (1) the head, which includes any function or terminal symbols, is obtained using the trial-and-error approach and (2) a tail portion that only has a terminal symbol, and its length is calculated through Equation (5). In this equation, t is the tail's length, h is the length of the head and nMAX is the maximum argument of the functions.
The stages of this model are schematically shown in Figure 10. As can be observed, firstly, chromosomes are randomly selected and make up the primary population. Then, these chromosomes appear as a tree-like structure and are evaluated using fitness functions (such as roulette wheel sampling) and are then combined to create the new (gene) population. The chromosomes in the newly created generation are again selected based on their performance relative to the environment and, after correction (amendment), create a new generation. These new children are developed through the same cycle as their parents, consisting of mutation, transposition and recombination [72,75,76].

Gene Expression Programming (GEP)
Gene expression programming (GEP) is an evolutionary algorithm first introduced by Feriera in 1999 [72]. This method is a genotype-phenotype genetic algorithm with high accuracy due to its tree-like structure, and it has made a great difference in comparison to genetic programming. Additionally, it has GA's simplicity and GEP's capability; consequently, it compensates for the shortcomings of these two methods [73,74]. Chromosome and tree structures are two entities of GEP algorithms in which chromosomes contain several genes. Each gene has two parts: (1) the head, which includes any function or terminal symbols, is obtained using the trial-and-error approach and (2) a tail portion that only has a terminal symbol, and its length is calculated through Equation (5). In this equation, t is the tail's length, h is the length of the head and n MAX is the maximum argument of the functions.
The stages of this model are schematically shown in Figure 10. As can be observed, firstly, chromosomes are randomly selected and make up the primary population. Then, these chromosomes appear as a tree-like structure and are evaluated using fitness functions (such as roulette wheel sampling) and are then combined to create the new (gene) population. The chromosomes in the newly created generation are again selected based on their performance relative to the environment and, after correction (amendment), create a new generation. These new children are developed through the same cycle as their parents, consisting of mutation, transposition and recombination [72,75,76]. Mutation plays a key role in these stages due to its internal correction ability. According to the fixed length of the chromosomes, any function or terminal symbols can be replaced with each other in the head, but in the tail, only the terminal can be replaced. Numerous researchers have suggested mutation values ranging from 0.01 to 0.1 in their studies [77,78]. At the same time, in the transposition stage, transportable particles are moved from one chromosome to another. This transposition of particles takes place in three states: (a) moving by replacing the head, which is called insertion sequence transposition, (b) moving to the root of the chromosome, which is called root insertion sequence transposition, and (c) moving a gene to the beginning of the chromosome, which is called gene transposition. The value of this transposition is considered to range between 0.01 and 0.1 according to previous researchers [72,75]. In the final stage, a recombination stage (also known as crossover) is executed on the chromosome in three different ways: one-point, two-point and gene-point recombination. According to Ferreira's studies, recombination's value is considered to be 0.7. This process is repeated for a certain number of generations until a suitable solution is obtained [72,75,76,79].

Fuzzy DEMATEL Technique
The decision-making trial and evaluation laboratory method is a practical scientific method for showing the cause-and-effect relationship between variables [80]. In this method, the quantities recognized by experts in a survey matrix are represented as a fuzzy set using a fuzzy triangle number (TFN) ( Table 5) [81][82][83][84][85][86][87]. Mutation plays a key role in these stages due to its internal correction ability. According to the fixed length of the chromosomes, any function or terminal symbols can be replaced with each other in the head, but in the tail, only the terminal can be replaced. Numerous researchers have suggested mutation values ranging from 0.01 to 0.1 in their studies [77,78]. At the same time, in the transposition stage, transportable particles are moved from one chromosome to another. This transposition of particles takes place in three states: (a) moving by replacing the head, which is called insertion sequence transposition, (b) moving to the root of the chromosome, which is called root insertion sequence transposition, and (c) moving a gene to the beginning of the chromosome, which is called gene transposition. The value of this transposition is considered to range between 0.01 and 0.1 according to previous researchers [72,75]. In the final stage, a recombination stage (also known as crossover) is executed on the chromosome in three different ways: one-point, two-point and gene-point recombination. According to Ferreira's studies, recombination's value is considered to be 0.7. This process is repeated for a certain number of generations until a suitable solution is obtained [72,75,76,79].

Fuzzy DEMATEL Technique
The decision-making trial and evaluation laboratory method is a practical scientific method for showing the cause-and-effect relationship between variables [80]. In this method, the quantities recognized by experts in a survey matrix are represented as a fuzzy set using a fuzzy triangle number (TFN) ( Table 5) [81][82][83][84][85][86][87]. Table 5. The most common TFNs between linguistic terms and linguistic values [82].

Linguistic Terms Linguistic Values
Very Here, an n × n primary non-negative matrix is created as the primary matrix (x k ij = (l ij , m ij , u ij )), in which i, m and u, respectively, indicate the lower, medium and higher TFN markers. Then, the average of the available matrixes is calculated and normalized. Then, a total matrix is obtained (Equation (6)), and each of its members is in the form of t ij = l ij , m ij , u ij . The weight of factors (W i ) is then calculated by taking the conditions of the problem into consideration [85]. All these steps are illustrated in a flowchart ( Figure 11).
Mining 2023, 3, FOR PEER REVIEW 14 Table 5. The most common TFNs between linguistic terms and linguistic values [82]. Here, an n × n primary non-negative matrix is created as the primary matrix ( = ( , , )), in which i, m and u, respectively, indicate the lower, medium and higher TFN markers. Then, the average of the available matrixes is calculated and normalized. Then, a total matrix is obtained (Equation (6)), and each of its members is in the form of = ( , , ). The weight of factors (Wi) is then calculated by taking the conditions of the problem into consideration [85]. All these steps are illustrated in a flowchart ( Figure  11).

Linguistic Terms Linguistic Values
Preparation of initial matrixes Gather the survey matrix (l , , )

Averaging the initial matrices
Normalize the mean matrixes Calculate the weight of factors Figure 11. The flowchart of fuzzy DEMATEL technique.
During these stages, all obtained numbers are fuzzy numbers. Therefore, a procedure is required to turn these numbers into non-fuzzy numbers. Here, the best non-fuzzy performance (BNP) method (Equation (7)) can be used to turn the values into non-fuzzy numbers [82,87].

Development of Predictive Models
In order to propose the above-explained models in this study, the gathered data were randomly divided into a testing and training set (90 blasts were chosen randomly from 118 blasts as a training set). On this basis, data from 90 blasts were used for the training set, and the rest were used for the testing and validation of the models. It should be noted that among the 90 blasts selected for the development of predictive models, 65 blasts were used for model training, and 25 blasts were used to validate the models for flyrock prediction. It must be mentioned that, in accordance with Table 4, the flyrock distance was selected as the output data, and other parameters affecting this phenomenon were selected as input data. For this aim, first, the data that were randomly selected as training data were classified. Second, using SPSS software, PCA was performed on the data, and During these stages, all obtained numbers are fuzzy numbers. Therefore, a procedure is required to turn these numbers into non-fuzzy numbers. Here, the best non-fuzzy performance (BNP) method (Equation (7)) can be used to turn the values into non-fuzzy numbers [82,87].

Development of Predictive Models
In order to propose the above-explained models in this study, the gathered data were randomly divided into a testing and training set (90 blasts were chosen randomly from 118 blasts as a training set). On this basis, data from 90 blasts were used for the training set, and the rest were used for the testing and validation of the models. It should be noted that among the 90 blasts selected for the development of predictive models, 65 blasts were used for model training, and 25 blasts were used to validate the models for flyrock prediction. It must be mentioned that, in accordance with Table 4, the flyrock distance was selected as the output data, and other parameters affecting this phenomenon were selected as input data. For this aim, first, the data that were randomly selected as training data were classified. Second, using SPSS software, PCA was performed on the data, and those whose variance was less than 0.5 were removed. Then, statistical analysis (normalization, skewness and kurtosis) was performed on the remaining parameters to examine the data. Finally, different permutations were formed between the selected parameters (12 parameters), and 19 permutation models (different combinations) that had no correlations with each other and could show the best possible solution for the mentioned methods were selected and are presented in Table 6. In order to choose a developed equation or model with a better performance in predicting the flyrock distance, the root-mean-square error (RMSE) and Correlation Coefficient (R 2 ) were used (Equations (8) and (9)) [88].  In the equations above, X msr is the actual value of data i, X pre is the model's output, x pre and x msr are the average values of predicted and actual data, respectively, and n is the total number of data points.
Here, 19 different combinations (Table 6) without any correlations or direct influence on each other were obtained to develop predictive models. These combinations were then used as input data. Finally, using the 19 different combinations (Table 6) mentioned, as well as the methods mentioned in Section 3, the flyrock prediction was performed.

Regression Model
To develop a regression model, 2 methods of linear and nonlinear regression were used to predict flyrock. Data that were randomly selected as training data were assessed for correlations and normal distributions using the SPSS software. Finally, among all the developed models (Table 7), model number 11 for linear regression and nonlinear regression showed the best performance compared to other trained linear and nonlinear models because of its higher and lower R 2 and RMSE, respectively (Figures 12 and 13). Therefore, this model was chosen and suggested in order to predict the flyrock distance in the Jajarm bauxite mine complex.

GRNN Model
In order to employ the GRNN for the flyrock distance prediction, coding was carried out in the Python programming language. Additionally, the Euclidean distance, which is the distance between the central vector of the middle-layer neuron and the input neuron, was determined for the different models. This distance is one of the most important controlling parameters used in the Gaussian transfer function in the GRNN method. Here, this distance was changed in 0.1 intervals from 0.1 to 10 in order to obtain the optimal network among different combinations. In this method, the optimal network is the one in which the values of RMSE and R 2 are the lowest and highest among all distances, respectively. According to Table 8 and Figure 14, model no. 16, with a Euclidean distance of 9.9, has the best functionality based on the GRNN method among the different combinations.

GRNN Model
In order to employ the GRNN for the flyrock distance prediction, coding was carried out in the Python programming language. Additionally, the Euclidean distance, which is the distance between the central vector of the middle-layer neuron and the input neuron, was determined for the different models. This distance is one of the most important controlling parameters used in the Gaussian transfer function in the GRNN method. Here, this distance was changed in 0.1 intervals from 0.1 to 10 in order to obtain the optimal network among different combinations. In this method, the optimal network is the one in which the values of RMSE and R 2 are the lowest and highest among all distances, respectively. According to Table 8 and Figure 14, model no. 16, with a Euclidean distance of 9.9, has the best functionality based on the GRNN method among the different combinations. Figure 13. Comparison of R 2 and RMSE for nonlinear regression models.

GRNN Model
In order to employ the GRNN for the flyrock distance prediction, coding was carried out in the Python programming language. Additionally, the Euclidean distance, which is the distance between the central vector of the middle-layer neuron and the input neuron, was determined for the different models. This distance is one of the most important controlling parameters used in the Gaussian transfer function in the GRNN method. Here, this distance was changed in 0.1 intervals from 0.1 to 10 in order to obtain the optimal network among different combinations. In this method, the optimal network is the one in which the values of RMSE and R 2 are the lowest and highest among all distances, respectively. According to Table 8 and Figure 14, model no. 16, with a Euclidean distance of 9.9, has the best functionality based on the GRNN method among the different combinations.

GA-GRNN Model
Since the genetic algorithm is an optimization algorithm, it was used to optimize the best GRNN model (model number 16, whose parameters are Mc, S/B, H/B, TD, St and BI) in order to reduce the flyrock prediction error. The optimization operator (genetic algorithm) used to optimize the Euclidean distance of the GRNN network was selected for model number 16. As stated earlier, in order to design a prediction model using the GA-GRNN method, similar to the coding of GRNN in the PYTHON programming language, the effective variables in this method must be obtained. The trial-and-error method is the best method for determining the network's input variables [42,89]. For this, the size of the primary population was set between 100 and 500. Additionally, after determining the primary population, the values for mutation, recombination and crossover were obtained to start the network cycle according to the upper and lower limits. The value for mutation equaled 0.2, and a crossover in the range of (0-1) was used as the input data with a 0.8 probability. On the other hand, the maximum generation (Gmax) was considered 1000 because this number led to the minimum RMSE value in calculating the optimal network and, therefore, represents the best model.
The best model and its respective results when varying the population size are shown in Table 9 and Figure 15. The desired network was trained for 0-450 iterations. As can be seen in Table 9 and Figure 15, this network was stabilized with different populations from 250 iterations onward, with some overfitting. As a result, the trained neural network

GA-GRNN Model
Since the genetic algorithm is an optimization algorithm, it was used to optimize the best GRNN model (model number 16, whose parameters are Mc, S/B, H/B, TD, S t and BI) in order to reduce the flyrock prediction error. The optimization operator (genetic algorithm) used to optimize the Euclidean distance of the GRNN network was selected for model number 16. As stated earlier, in order to design a prediction model using the GA-GRNN method, similar to the coding of GRNN in the PYTHON programming language, the effective variables in this method must be obtained. The trial-and-error method is the best method for determining the network's input variables [42,89]. For this, the size of the primary population was set between 100 and 500. Additionally, after determining the primary population, the values for mutation, recombination and crossover were obtained to start the network cycle according to the upper and lower limits. The value for mutation equaled 0.2, and a crossover in the range of (0-1) was used as the input data with a 0.8 probability. On the other hand, the maximum generation (G max ) was considered 1000 because this number led to the minimum RMSE value in calculating the optimal network and, therefore, represents the best model.
The best model and its respective results when varying the population size are shown in Table 9 and Figure 15. The desired network was trained for 0-450 iterations. As can be seen in Table 9 and Figure 15, this network was stabilized with different populations from 250 iterations onward, with some overfitting. As a result, the trained neural network (model number 16) with a population of 350 and best fitness = 0.0203 had the best performance in iteration 250. A Euclidean distance of 9.906 was selected for the GRNN network when using the genetic algorithm because it resulted in the lowest RMSE and the highest R 2 ( Figure 15).

GEP Model
In order to design a prediction model using the GEP network, a set of regulations are needed in the modeling process when coding in the Python programming language to obtain the best model. These regulations are usually determined through the trial-anderror method. The number one priority is determining the number of genes that can affect the network's accuracy, as a higher number of genes leads to lengthier functions with a lack of performance. It is important not to choose very low numbers as well so that the accuracy and performance of the model do not decline. The number of genes and other parameters used in this network are shown in the summary in Table 10.

GEP Model
In order to design a prediction model using the GEP network, a set of regulations are needed in the modeling process when coding in the Python programming language to obtain the best model. These regulations are usually determined through the trial-and-error method. The number one priority is determining the number of genes that can affect the network's accuracy, as a higher number of genes leads to lengthier functions with a lack of performance. It is important not to choose very low numbers as well so that the accuracy and performance of the model do not decline. The number of genes and other parameters used in this network are shown in the summary in Table 10.  Based on the information in Table 10, the GEP network was applied for all different combinations to predict the flyrock distance. According to Table 11 and Figure 16, model no. 17 is the best equation among all tested equations due to its lower RMSE and higher R 2 value. In addition, Figure 17 is the tree representation of model no. 17, which demonstrates the relationship between effective variables in predicting the flyrock distance. Table 11. The best GEP models to predict flyrock distance with values of RMSE and R 2 .

Fuzzy DEMATEL Model
In order to predict flyrock distance using the fuzzy DEMATEL method, 15 surveys concerning variables directly affecting flyrock were composed and distributed among mining experts. The experts were asked to grade the effect of each variable on the flyrock distance from 0, the least, to 4, the most. The given grades were converted into fuzzy numbers using Table 5, and the primary equations' direct matrix was formed. Then, averaging was performed from the primary matrix, and the average matrixes ( ) were obtained. In the next step, the average matrixes were normalized, and based on these, the general fuzzy relation matrixes for the upper, mean and lower limits were constructed.

Fuzzy DEMATEL Model
In order to predict flyrock distance using the fuzzy DEMATEL method, 15 surveys concerning variables directly affecting flyrock were composed and distributed among mining experts. The experts were asked to grade the effect of each variable on the flyrock distance from 0, the least, to 4, the most. The given grades were converted into fuzzy numbers using Table 5, and the primary equations' direct matrix was formed. Then, averaging was performed from the primary matrix, and the average matrixes ( A) were obtained. In the next step, the average matrixes were normalized, and based on these, the general fuzzy relation matrixes for the upper, mean and lower limits were constructed.
After analyzing the importance of variables using the upper, mean and lower limit matrixes in terms of quality, it is vital for them to be analyzed in terms of quantity, which is necessary for predicting the flyrock distance. Therefore, the weight of each variable was calculated and is shown in Table 12 as both the fuzzy weight and deterministic weight.
In addition, Figure 18 shows the ranking of each variable compared to others based on its weight.

Proposed Fuzzy DEMATEL Method
In order to introduce a prediction model based on the fuzzy DEMATEL method, all parameters affecting flyrock were classified and ranked using Equation (10), and the ranges of parameters were proposed based on other researchers' obtained results [6,55,90]. This ranking is shown in Table 13.

Proposed Fuzzy DEMATEL Method
In order to introduce a prediction model based on the fuzzy DEMATEL method, all parameters affecting flyrock were classified and ranked using Equation (10), and the ranges of parameters were proposed based on other researchers' obtained results [6,55,90]. This ranking is shown in Table 13. Based on this, the value of the flyrock distance was determined through Equation (11): (11) in which W i is the ith parameter's weight, P i is the ith parameter's range from 0 to 4 and P max is the maximum range of each parameter in its respective classification. Now, in order to accurately assess the flyrock distance using the fuzzy DEMATEL method, the measured flyrock distance from the Jajarm bauxite mining complex was determined through the FRI function (Equation (12)). Using the Matlab software, the best curves for the data obtained from Equation (12) were drawn. Among the obtained equations (Table 14), the one with a higher R 2 and lower RMSE was chosen as the equation to measure the flyrock distance using the fuzzy DEMATEL method. According to Table 14, the Gaussian equation has the highest value of R 2 and the lowest value of RMSE and was chosen as the prediction model for the flyrock distance.

Evaluation of Proposed Models
In the past, knowledge-driven methods such as DEMATEL have been used. Nowadays, it is common to use data-driven methods in order to reduce the uncertainty and obtain more accurate results. In this case study, in addition to using the knowledge-driven method, we used data-driven methods for the first time to obtain more accurate results by reducing the uncertainty. In addition, there are optimization algorithms that can be used to increase the performance of data mining methods. In this study, genetic algorithms were used for optimization.
The data from 28 blasts that had no role in the development of the proposed models (Section 4) were used to predict and evaluate the functionality of the developed equations (Table 15) in order to evaluate the performance of the best proposed model; a comparison table (Table 15) is shown that demonstrates the flyrock distances calculated by the proposed models and the actual values measured in the field. Due to a lack of reliability and a comparison of the results in table data, in order to choose a practical equation with a better performance in predicting the flyrock distance, except for Equations (8) and (9), the variance accounted for (VAF) and mean absolute percentage error (MAPE) were used (Equations (13) and (14)) [88].
In the equations above, X msr is the actual value of data i, X pre is the model's output, and n is the total number of data points.
According to Table 16 and Figure 19, it can be understood that a complete and ideal prediction model is one that has RMSE and MAPE values equal to zero and R 2 and VAF values equal to 100 percent. Therefore, higher values of R 2 and VAF and lower values of RMSE and MAPE show the appropriate and optimal performance of a model, respectively [55,79,88,91].  In the equations above, Xmsr is the actual value of data i, Xpre is the model's output, and n is the total number of data points.
According to Table 16 and Figure 19, it can be understood that a complete and ideal prediction model is one that has RMSE and MAPE values equal to zero and R 2 and VAF values equal to 100 percent. Therefore, higher values of R 2 and VAF and lower values of RMSE and MAPE show the appropriate and optimal performance of a model, respectively [55,79,88,91]. Based on Table 16 and Figure 19, it can be understood that the value predicted by the GA-GRNN model is closer to the actual value in the training and testing stages. Therefore, it can be concluded that the aforementioned model is acceptable in predicting the flyrock distance and can be considered the predictive model for the Jajarm bauxite mining complex.

Conclusions
Given the fact that flyrock in mines is the leading cause of danger and casualties and Based on Table 16 and Figure 19, it can be understood that the value predicted by the GA-GRNN model is closer to the actual value in the training and testing stages. Therefore, it can be concluded that the aforementioned model is acceptable in predicting the flyrock distance and can be considered the predictive model for the Jajarm bauxite mining complex.

Conclusions
Given the fact that flyrock in mines is the leading cause of danger and casualties and damage to equipment, accurately predicting it is vital in preventing possible harm. To substantiate this matter, in this study, data from 118 blasts were gathered and divided into training (90 blasts) and test (28 blasts) data. On the basis of these data, flyrock was considered the output, and different combinations of effective parameters of flyrock were selected as the input data. Three approaches, soft computing, statistical and multicriteria decision-making methods, were selected and evaluated according to four statistical evaluation criteria to predict the flyrock distance in Jajarm bauxite mines and obtained the following results:

3.
Fuzzy DEMATEL was the last approach applied to the parameters affecting flyrock. In this regard, the amount of flyrock was predicted based on the function of the flyrock, the best of which was a Gaussian function. The performance of this function on the training data was R 2 = 87.171, RMSE = 34.682. In addition, its performance was R 2 = 48.08 and RMSE = 31.415 on the test data.

4.
Finally, the best models obtained from each method were selected and evaluated and compared to each other. The results of these evaluations showed that the GA-GRNN model's performance was far better compared to others and has reasonable conformity to actual values due to its higher VAF and lower MAPE and RMSE, and it represents a direct practical equation.