SSA was inspired by the foraging behavior and anti-predation behavior of sparrows (51) and has been applied to solve many complex engineering optimization problems (52–55). A brief description of SSA is as follows.

Sparrow population members are divided into three categories: finders, entrants, and scouts. The role of finders is to find food, providing all entrants with a foraging area and direction. The position update of the finder during each iteration is described by Equation (1). (1)Xijt+1={Xijt·exp(-iα·Tmax),R2<STXijt+Q·L,R2≥ST where t represents the current iteration number, T_max is the maximum number of iterations, Xijt is the position of the i sparrow in the j dimension, α is a uniform random number between [0, 1], Q is a random number that follows a standard normal distribution, L is a 1 × d matrix whose elements are all 1. R₂ ∈ [0, 1], represents the warning value; ST ∈ [0.5, 1], represents the safety value.

Entrants will follow the finders for foraging in order to obtain food, and its position update is described by Equation (2). (2)Xijt+1={Q·exp(Xwt-Xijti2), i>N2Xpt+1+|Xjit-Xpt+1|A+L, otherwise where Xpt+1 represents the best position of the finder at the t + 1 iteration; Xwt represents the global worst position at the t iteration; A is a matrix of size 1 × d with elements randomly assigned to 1 or −1, and A⁺ = A^T(AA^T)⁻¹; N is the population size.

Scouts are responsible for reconnaissance and early warning, alerting the entire population to give up foraging when danger is detected. The position update of the scout is described by Equation (3). (3)Xijt+1={Xbt+β·|Xijt-Xbt|,fi>fgXijt+K·(|Xijt-Xwt|(fi-fw)+ε),fi=fg Where Xbt represents the global optimal position at the t iteration; β is the step size parameter, which obeys the normal distribution random number with mean 0 and variance 1; K is a random number in the range [−1, 1]. f_i, f_w and f_g are the individual fitness value, the global worst fitness value, and the global optimal fitness value of the current sparrow, respectively. ε is an extremely small constant that prevents the denominator from being zero.

GWO is a novel swarm intelligence optimization algorithm inspired by the group predation behavior of gray wolves (56), which has been applied in many engineering fields (57–60). There is a strict hierarchy within the gray wolf group, and the entire wolf group is divided into four levels, namely α, β, δ, and ω. As shown in Figure 1, α is the optimal gray wolf, which is the leader of the wolf group and has the right to decide all major issues of the entire wolf group; β is the second-best gray wolf, assisting the leader wolf to make decisions; δ is at the third level, responsible for sentry, reconnaissance and other tasks; ω is the lowest level wolf and is under the command of the first three levels of gray wolves in action. More descriptions of the principle of GWO can be found in the literature (56, 61).

The predation process of gray wolves is divided into three processes: search, surround, and attack. During the hunting process, the individual gray wolf realizes the update of the position according to the following formula to round up the prey. (4)D=|C→·X→p(t)-X→p(t)| (5)X→(t+1)=X→p(t)-A→·D (6)A→=2a→·r→1-a→ (7)C→=2·r→2 Where, t is the current number of iterations; D is the distance between the individual gray wolf and its prey; X→ is the position vector of the individual gray wolf; X→p is the position vector of the prey; A→ and C→ are coefficient vectors; a→ is the convergence factor that decreases linearly from 2 to 0 with the number of iterations; r→1 and r→2 are vectors of random numbers between [0,1].

BOA is a metaheuristic algorithm proposed based on the foraging and courtship behavior of butterflies (62). Assume that each butterfly releases a certain intensity of scent and the concentration of the released scent is related to its adaptation, while each butterfly senses the scent of other butterflies around it and moves toward those that emit more scent. The intensity of the fragrance produced by butterflies is represented by Equation (8). (8)fi=cTα where f_i is the scent intensity perceived by the i-th individual butterfly, c is the sensory modality of the butterfly, T is the stimulus intensity, α is a power exponential parameter that depends on the sensory modality.

When the individual butterfly feels that a certain butterfly emits more fragrance in this area, it will conduct a global search toward the source of the fragrance. This process is represented by Equation (9). (9)xit+1=xit+r∗(g*-xit)∗fi where xit is the solution vector of the i-th butterfly in the t-th iteration; g^* represents the optimal solution among all solutions in the current iteration; r is a random number between [0, 1].

When a butterfly cannot perceive a scent higher than its own, it will fly randomly and perform a local search. The local search can be expressed by Equation (10). (10)xit+1=xit+r∗(xkt-xjt)∗fi where xjt and xkt represent the solution vectors of the j-th and k-th butterflies in the solution space of the t-th iteration. If xjt and xkt belong to the same population, and r is a random number within (0, 1), it means a local random walk.

TSA is a metaheuristic algorithm based on the jet propulsion and swarm behavior of tunicates in the ocean during foraging (63). Tunicates use two of their own behaviors to find food sources, jet propulsion and swarm intelligence. A specific description of the principle of the algorithm can be found in the literature (63). For mathematical modeling of jet propulsion behavior, tunicates should satisfy the following three conditions.

Avoid conflicts between search populations. To avoid the conflict between tunicates, TSA calculates the new tunicate position through Equation (11). (11)X→=c2+c3-2c1|Vmin+c1·Vmax-Vmin| Where X→ is the new tunic position; V_min and V_max are the minimum and maximum initial speed of social interaction, respectively, and are generally set as [1, 4]; c₁, c₂, c₃ are random numbers between [0,1].

Move to the position of the best search individual. (12)D→=|F→-r·P→(x)| where D→ is the distance between the food and the searched individual; F→ is the position of the table food; x is the number of iterations; P→(x) represents the position of the tunicate individual in the x-th iteration. r is a random number between [0,1].

Converge to the optimal position. (13)P→(x*)={F→+X→·D→,r≥0.5F→-X→·D→,r<0.5 where P→(x*) is the updated position.

After modeling the individual jet propulsion behavior, the swarm behavior is modeled. TSA saves the first two optimal solutions and updates the positions of other search individuals according to the position of the best search individual. Group behavior can be defined using Equation (14). (14)P→(x+1)=P→(x)+P→(x+1)2+c1

MVO is a metaheuristic optimization algorithm, which has been widely used in many fields (64–67). It is based on the principle that matter in the universe is simulated by transferring from a white hole to a black hole through a wormhole, which can be described in the literature (64). Here are some brief introductions about MVO.

Assuming that the universe matrix exists in the search space is: (15)U=[x11x12…x1dx21x22…x2d…………xn1xn2…xnd] Where d is the spatial dimension, n is the number of universes.

Due to the different inflation rates of each individual universe, the individuals in the universe are able to transfer through white or black holes, chosen randomly using the roulette wheel method, as shown in Equation (16). (16)xij={xkjr1<NI(Ui)xkjr1≥Nl(Ui) Where xij is the j-th variable of the i-th universe, xkj is the j-th variable of the i-th universe selected by a roulette wheel selection mechanism, U_i shows the i-th universe, NI(U_i) is a normalized inflation rate of the i-th universe, r₁ is a random number in [0,1].

The random transport of matter between universes through wormholes to ensure population diversity and the movement of individual universes toward the current optimal universe to increase the inflation rate is given by the following equation. (17)xij={{Xj+TDR∗((ubj−lbj)∗r4+lbj)r3<0.5Xj−TDR∗((ubj−lbj)∗r4+lbj)r3≥0.5xji r2≥WEP where X_j is the j-th variable of the optimal position found so far, ub_j and lb_j respectively represent the upper and lower boundaries of the search space where the j-th parameter is located, r₂, r₃, and r4 are all random numbers between [0, 1]. WEP and TDR are the wormhole existence probability and the travelling distance rate, respectively, and their expressions are as follows. (18)WEP=WEPmin+l∗(WEPmax-WEPminL) (19)TDR=1-l1/pL1/p Where WEP_min and WEP_max are the minimum and maximum values of parameter WEP, l, and L are the current number of iterations and the maximum number of iterations, p is the exploit precision in the algorithm.

The back propagation neural network (BPNN) is a multi-layer network that forwards signals and propagates errors backwards. BP neural network has been studied and applied to solve many problems (68–72), and is one of the most widely used networks at present. It can simulate the information transmission mode of human brain neurons, perform non-linear transformation and regression processing on complex information variables, and obtain operation results with a high fitting degree. The structure of the BP neural network consists of three layers: input layer, hidden layer, and output layer (52). The number of hidden layers is not fixed, it can be one layer or multiple layers.

Suppose the number of nodes in the input layer, hidden layer, and output layer are l, m, n. The weight from the input layer to the hidden layer is ω_ij, the weight from the hidden layer to the output layer is ωjk ′, the threshold from the input layer to the hidden layer is q_j, and the weight from the hidden layer to the output layer is qk ′. The excitation function g(t) is taken as a Sigmoid function with the following equation. (20)g(t)=11+e-t Then the output of the neuron in each layer is: (21){Hj=g(∑i=1lωijxi-qj),j=1,2,…,mOk=f(∑j=1mωjk ′Hj-qk ′),k=1,2,…,n where H_j is the output of the jth neuron in the hidden layer, O_k is the output of the neural network, and x_i is the input data.

Just like the information transmission between neurons in the brain, after the input variables are input in the input layer, a linear combination of the input variables will be obtained according to the initially set weights. When the weights are continuously modified such that the linear combination value exceeds the threshold, the information is passed to the output layer (73). Figure 2 illustrates the BPNN model framework for predicting induced stress at different locations within the pillar in this study based on 5 input variables and 1 output variable.

The data used in this study refer to reference (2), and the database used is from several coal mines whose geological and geometrical conditions are basically the same and which are excavated by bord and pillar method. The main purpose of this study is to predict the induced stress at different locations within the same pillars. Figure 3 shows the plan of the sub-panel, in which Arabic numerals are used to represent different locations within the same pillars. The parameters affecting the induced stress level in the point are the angle of goafline (A), depth of the working coal seam (H), specific gravity of the coal (G), distance of the point from the center of the pillar (C), and distance of the point from goafline (D). These five parameters have been taken as input variables and the induced stress (S) within the coal pillar point as the output variable in this study.

A total of 149 datasets were collected in this study to predict stress magnitudes. Table 1 presents some statistics about the input and output parameters (mean, range, and standard deviation) along with their units and descriptions to better illustrate the used data samples. In addition, in order to train and test the model, the database is randomly divided into a training set and a test set, and their data distribution status is represented by a violin plot, as shown in Figure 4. The violin plot is a combination of a box plot and a kernel density plot, where the black boxes indicate the range from the lower quartile to the upper quartile and the red dots in the middle represent the mean (61). From the figure, it can be seen that the parameter data distributions of the training set and testing set are approximately the same in order to ensure that the results obtained in this study will be reliable.

The overall analysis and modeling process of the five hybrid models is shown in Figure 5. All data sets were normalized to the range [−1, 1] in order to improve prediction accuracy and avoid redundant computational costs. To build the mixture model effectively, the collected database was randomly divided into 80% training set and 20% test set based on the Pareto principle (74). The weights and thresholds of the BPNN were then optimized using five metaheuristic algorithms, and the mean absolute error (MAE) was used as the fitness value in the optimization process to determine whether to stop, and the optimal weights and thresholds were assigned to the network predictions. Finally, the predictive performance of the developed neural-metaheuristic hybrid model was evaluated for the training and test sets using five evaluation indicators (R², RMSE, VAF, MAE, and MAPE).

To determine the optimal weights and thresholds to obtain the best column stress prediction performance, five meta-heuristic algorithms (SSA, GWO, BOA, TSA, and MVO) were used to tune the weights and thresholds in the BPNN model, and the mean absolute error (MAE) was used as the fitness. In this study, six population sizes of 50, 100, 150, 200, 250, and 300 were designed for each of the five hybrid models. the maximum number of iterations was set to 500 to select the best optimal parameters. Figures 6A–E shows the fitness curves of the five hybrid models with different population sizes, and it can be seen that the optimal process is different for different population sizes, and the MAE value decreases with the increase of the number of iterations. The lowest fitness values of SSA-BPNN, GWO-BPNN, BOA-BPNN, TSA-BPNN, and MVO-BPNN were 1.0827, 0.92683, 1.03569, 1.10515, and 0.99775, corresponding to the mixed models with population sizes of 200, 300, 200, 200, and 250, in that order. It indicates that these hybrid models perform best in the training set. Of course, the fitness curve alone is not sufficient to determine the optimal population size. Therefore, to determine the optimal population size, we evaluated the model performance with five indicators and calculated the GPI values of each model based on the evaluation indicators, and the results are shown in Table 2. Better models have higher GPI values, so the optimal population size of SSA-BPNN is 50 (Training set: R² = 0.9984, RMSE = 0.2009, VAF = 99.84, MAE = 0.1262, and MAPE = 0.8477; Test set: R² = 0.9980, RMSE = 0.1846, VAF = 99.80, MAE = 0.1400, and MAPE = 1.0736), the optimal population size of GWO-BPNN is 300 (training set: R^{2 =} 0.9991, RMSE = 0.1535, VAF = 99.91, MAE = 0.0884, MAPE = 0.6107; test set: R^{2 =} 0.9983, RMSE = 0.1783, VAF = 99.83, MAE = 0.1230, and MAPE = 0.9253), the optimal population size of BOA-BPNN is 200 (training set: R² = 0.9988, RMSE = 0.1750, VAF = 99.88, MAE = 0.0934, and MAPE = 0.6607; test set: R² = 0.9973, RMSE = 0.2207, VAF = 99.71, MAE = 0.1543, and MAPE = 1.1900), and the optimal population size of TSA-BPNN is 150 (Training set: R² = 0.9989, RMSE = 0.1667, VAF = 99.89, MAE = 0.1041, and MAPE = 0.7292; Test set: R² = 0.9975, RMSE = 0.2234, VAF = 99.75, MAE = 0.1642, and MAPE = 1.2421), the optimal population size of MVO-BPNN is 250 (training set: R² = 0.9989, RMSE = 0.1692, VAF = 99.88, MAE = 0.1076, and MAPE = 0.7078; test set: R² = 0.9968, RMSE = 0.2498, VAF = 99.67, MAE = 0.1754, and MAPE = 1.2480).

The correlations between predictions and measured induced stresses for the five optimal mixed model training and test datasets are shown in Figure 7. It can be seen from the figure that the training and testing sample points of these hybrid models are basically distributed near the perfect fitting line (“measured stress value = predicted stress value”), and the R² values are above 0.99, indicating that the five optimization techniques based on BP neural network proposed in this paper can achieve high training and testing outcome results.

In order to compare and evaluate the prediction performance of the hybrid models, based on the previous discussion, three classical machine learning models, BPNN, random forest (RF) and radial basis function network (RBF), were also introduced in this study for comparison, and the GPI values of each model were calculated, and the detailed data are listed in Table 3. from the table, it can be seen that the GWO-BPNN model has the best prediction performance (on the training set, R² = 0.9991, RMSE = 0.1535, VAF = 99.91, MAE = 0.0884, and MAPE = 0.6107; on the test set, R² = 0.9983, RMSE = 0.1783, VAF = 99.83, MAE = 0.1230, and MAPE = 0.9253). Figure 6F gives an iterative graph of the objective optimization of the five hybrid models, from which it can be seen that the GWO-BPNN hybrid model shows the best results on the training set. In addition, to further demonstrate the superiority of the hybrid model, the Taylor graph is utilized to more visually show the performance of the model on the test set (as shown in Figure 8). The Taylor graph shows the correlation coefficient, RMSE, and standard deviation between the actual and predicted values of the model (65, 75, 91–93). The Taylor graph shows that the hybrid model significantly outperforms the un-optimized classical model. In conclusion, while all hybrid models performed well in predicting induced stress with high accuracy, the GWO-BPNN showed the best predictive performance overall. It is therefore recommended to apply GWO-BPNN to predict the induced stress at different locations within the pillar.

It is important to note that the data used in this article are from a study conducted with Monjezi et al. (2). They developed an ANN model to predict induced stress at different locations in pillar and pillar recovery. They report that in the testing phase, the optimal ANN model has an R² value of 0.988, an RMSE value of 0.056233, an MAE value of 0.309281, and a MAPE value of 2.337007. Compared with this study, the proposed GWO-BPNN model has an R² value of 0.9983, RMSE value of 0.1783, MAE value of 0.1230, and MAPE value of 0.9253, which is better than the ANN model developed by Monjezi et al. (2).