*Article* **Mining Plan Optimization of Multi-Metal Underground Mine Based on Adaptive Hybrid Mutation PSO Algorithm**

**Yifei Zhao, Jianhong Chen, Shan Yang \* and Yi Chen**

School of Resource and Safety Engineering, Central South University, Changsha 410083, China; zhaoyifei@csu.edu.cn (Y.Z.); jhchen@csu.edu.cn (J.C.); dlx8529@csu.edu.cn (Y.C.)

**\*** Correspondence: yangshan@csu.edu.cn

**Abstract:** Mine extraction planning has a far-reaching impact on the production management and overall economic efficiency of the mining enterprise. The traditional method of preparing underground mine production planning is complicated and tedious, and reaching the optimum calculation results is difficult. Firstly, the theory and method of multi-objective optimization are used to establish a multi-objective planning model with the objective of the best economic efficiency, grade, and ore quantity, taking into account the constraints of ore grade fluctuation, ore output from the mine, production capacity of mining enterprises, and mineral resources utilization. Second, an improved particle swarm algorithm is applied to solve the model, a nonlinear dynamic decreasing weight strategy is proposed for the inertia weights, the variation probability of each generation of particles is dynamically adjusted by the aggregation degree, and this variation probability is used to perform a mixed Gaussian and Cauchy mutation for the global optimal position and an adaptive wavelet variation for the worst individual optimal position. This improved strategy can greatly increase the diversity of the population, improve the global convergence speed of the algorithm, and avoid the premature convergence of the solution. Finally, taking a large polymetallic underground mine in China as a case, the example calculation proves that the algorithm solution result is 10.98% higher than the mine plan index in terms of ore volume and 41.88% higher in terms of economic efficiency, the algorithm solution speed is 29.25% higher, and the model and optimization algorithm meet the requirements of a mining industry extraction production plan, which can effectively optimize the mine's extraction plan and provide a basis for mine operation decisions.

**Keywords:** multi-objective optimization; mining plan; metal mines; adaptive; hybrid mutation

**MSC:** 90C29

#### **1. Introduction**

The preparation of the extraction plan, as a basic link in the production operation of a mining enterprise, is one of the most critical tasks in the production decision of a mine, and the rationality of the plan preparation directly affects the efficiency of the subsequent production links and the overall economic efficiency of the mining enterprise [1–3]. The traditional manual preparation method is not only time-consuming and intensive but also has poor accuracy and is difficult to modify. The reason for this is mainly the complex underground mine conditions during the preparation of the plan, which requires comprehensive consideration of the spatial and temporal constraints between the production processes and the mines and their continuity. Therefore, how to develop the underground mine production plan quickly and accurately has been an urgent problem to be solved.

With the continuous advancement of computer technology and operations research theory, many researchers have started to try to use the powerful simulation computing power of computers to simulate the mine production process, so as to continuously optimize the mine extraction plan [4,5]. Several researchers recognize that integer programming

**Citation:** Zhao, Y.; Chen, J.; Yang, S.; Chen, Y. Mining Plan Optimization of Multi-Metal Underground Mine Based on Adaptive Hybrid Mutation PSO Algorithm. *Mathematics* **2022**, *10*, 2418. https://doi.org/10.3390/ math10142418

Academic Editors: Wen Zhang, Xiaofeng Xu, Jun Wu and Kaijian He

Received: 20 June 2022 Accepted: 8 July 2022 Published: 11 July 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

can solve discrete production scheduling decision problems in the mining industry [6,7]. Many studies on mine production planning related to integer programming theory were subsequently carried out [3,8–12]. Dimitrakopoulos and Ramazan [13] developed an optimization framework for stochastic mine production scheduling considering mine uncertainties based on an ore body model and an integer planning approach. Weintraub et al. [14] developed a large, aggregated integer planning model (MIP) based on cluster analysis for mine planning at CODELCO, a national copper mine in Chile, through which data information of all CODELCO mines can be obtained to optimize the mine extraction planning. Newman et al. [15] developed a mixed integer planning model for underground mining operations in Kiruna mine, Sweden. This optimization model identifies an operationally feasible recovery sequence that minimizes deviations from the planned production quantities. Terblanche and Bley [16] used a mixed integer planning approach to construct a theoretical model applicable to the optimization of extraction production plans for open pit and underground mines but did not verify the feasibility of the model. Nehring et al. [17] proposed an improved modified model formulation for the classical model of long-term mine planning, assigning different human resources and equipment to each mining area. In the classical model, only one binary variable is assigned to each activity in each mining area, whereas the improved model assigns one binary variable to all activities under a more stringent assumption.

Although there are more existing studies on mathematical planning methods, most of these models achieve the solution of mine production planning with a single economic indicator; however, since mine extraction production planning is a complex system of engineering, it is difficult to highlight the production plan preparation and optimization effect by considering only a single economic indicator [18]. In order to overcome these difficulties, researchers have used a variety of computational intelligence methods to solve multi-objective prediction and optimization problems, and heuristic algorithms are effective methods to improve the solution speed and avoid involving local optimal solutions, while having unique advantages for multi-objective optimization problems (MOP). Little and Topal [19] investigated the methodology of whole life of mine (LOM) production planning generated using a simulated annealing technique and stochastic simulation representation of the ore body with the objective of maximizing the net present value (NPV) of the mine. Hou et al. [20] addressed the production planning for the next three years of the mine using an artificial bee colony optimization algorithm. Otto and Bonnaire [21] developed a "Greedy randomized adaptive search" program to help solve models for copper mining development and improve the speed of solution. O'Sullivan and Newman [22] proposed an optimization heuristic algorithm to set a complex set of constraints based on an optimization algorithm in a mining operation model for an underground lead–zinc mine in Ireland. Wang et al. [23] proposed a multi-objective optimization model formulation by taking the grade of mined and processed ore as the main constraint, maximizing mining returns and efficient use of natural resources as the objective function, and used a genetic algorithm to find the optimal solution to the multi-objective optimization problem. Nesbitt et al. [24] considered the uncertainty faced in the economic value of minerals faced by mines with long operating cycles. For a hard rock mine, the method of creating a stochastic integer program helps the mine to customize a mining schedule with a high degree of feasibility.

However, the use of traditional heuristic algorithms in the preparation of mining industry extraction production plans leads to slow solution speed and reduced global convergence performance, which brings many difficulties to the preparation and optimization effect of the actual production operation plan in real time and has a profound impact on the production management and economic efficiency of the enterprise. Therefore, in order to solve these problems, which are due to the complexity of multi-metal mine production planning, and since it is difficult to achieve the optimal results by traditional methods, this paper takes the best economic efficiency, grade, and ore volume as the goal and integrates constraints such as ore grade fluctuation, ore output of mining sites, mine production capacity, and mineral resource utilization. A production planning model is established, and the model is optimally solved by an improved particle swarm optimization (PSO) algorithm with nonlinear inertia weights and adaptive mutation probability (NAMPSO) in the context of an engineering example to improve the model's solving speed and increase the global convergence performance of the algorithm, thereby verifying the feasibility of the model solution method.

#### **2. Balanced Mining Plan Model**

#### *2.1. Model Analysis*

There are many factors that need to be considered comprehensively in the preparation of an underground mine extraction plan, and by analyzing the current situation faced by underground mining and the existing mining technology conditions [25–30], the main factors that need to be considered in the process of extraction plan preparation are shown in Figure 1.


**Figure 1.** Factors affecting the preparation of the mining plan.

In principle, each underground metal mine plan should consider all the factors in Figure 1, but in practice some factors can be combined and simplified according to the actual situation of the specific mine. For example, only one of the planned mining quantities and planned metal quantity are required, and only the smaller of the mine transport and hoisting capacity and mine processing plant capacity should be considered. At the same time, it is possible to add other influencing factors.

The process of preparing an extraction plan based on balanced mining is essentially a process of seeking an extraction solution that meets the requirements of various stakeholders in the process of comprehensive coordination of various constraints. Therefore, with the help of the objective optimization method in operations research theory, the complex extraction planning can be transformed into an MOP problem model by converting various restrictions into constraints and the final objective into an objective function, thus turning the complex programming process into a mathematical optimization problem.

#### *2.2. Model Building*

#### 2.2.1. Model Assumptions

This modeling is carried out based on the following three basic assumptions:


#### 2.2.2. Model Parameters Definition

This model parameters are defined at the end of the paper.

#### 2.2.3. Objective Functions

Equation (1) indicates the maximization of corporate revenue: from the perspective of corporate profitability, one of the purposes of mining enterprises conducting extractive operations is to obtain economic benefits, so the maximization of corporate revenue is one of the main objectives of the preparation of the extraction plan.

$$\text{argmax}Q = \sum\_{h=1}^{H} \sum\_{b=1}^{B} x\_b u e\_{hh} s\_h p\_h + \sum\_{h=1}^{H} \sum\_{r=1}^{R} y\_r e\_{rh} s\_h p\_h + \sum\_{h=1}^{H} \sum\_{k=1}^{K} z\_k e\_{kh} s\_h p\_h \tag{1}$$

$$\min E\_h = |\sigma\_h - \sigma\_h^\*| = \left| \frac{\sum\_{b=1}^{B} x\_b u e\_{bh} + \sum\_{r=1}^{R} y\_r e\_{rh} + \sum\_{k=1}^{K} z\_k e\_{kh}}{\sum\_{b=1}^{B} x\_b u + \sum\_{r=1}^{R} y\_r + \sum\_{k=1}^{K} z\_k} - \sigma\_h^\* \right| \tag{2}$$

$$\max A = \sum\_{b=1}^{B} x\_i u + \sum\_{r=1}^{R} y\_r + \sum\_{k=1}^{K} z\_k \tag{3}$$

Equation (2) indicates the grade optimization: from the perspective of mineral processing, due to the variability of mineral resource endowment, the development of the extraction plan needs to take full account of such differences to achieve the minimum fluctuations in the average grade of the extracted ore.

Equation (3) indicates the ore production maximization: the more ore extracted during the planning period, the higher the production efficiency and the better the mining operation.

#### 2.2.4. Constraints

After establishing the target of the extraction plan, it is necessary to transform each restriction into the constraint of the target function. Combined with the actual underground mining, the constraints are mainly the following six kinds.

$$\sum\_{b=1}^{B} \mathbf{x}\_b \boldsymbol{\mu} \boldsymbol{\omega} \boldsymbol{e}\_{\text{hh}} \mathbf{s}\_{\text{lh}} + \sum\_{r=1}^{R} \mathbf{y}\_r \mathbf{e}\_{r\text{h}} \mathbf{s}\_{\text{li}} + \sum\_{k=1}^{K} z\_k \mathbf{e}\_{k\text{h}} \mathbf{s}\_{\text{li}} \succ S\_{\text{li}} \tag{4}$$

$$\begin{cases} \ x\_b \le X\_b \\ \ y\_r \le Y\_r \\ z\_k \le Z\_k \end{cases} \tag{5}$$

$$\begin{cases} \sum\_{b=1}^{B} \boldsymbol{x}\_b \le D\_1\\ \sum\_{r=1}^{R} \boldsymbol{y}\_r \le D\_2\\ \sum\_{k=1}^{K} \boldsymbol{z}\_k \le D\_3 \end{cases} \tag{6}$$

$$\sum\_{b=1}^{B} lL\mathbf{x}\_{b} \ge (\sum\_{r=1}^{R} y\_{r} + \sum\_{k=1}^{K} z\_{k}) \times \eta \tag{7}$$

$$\sum\_{k=1}^{B} \varkappa\_{i}\mu + \sum\_{r=1}^{R} \mathcal{y}\_{r} + \sum\_{k=1}^{K} z\_{k} \le A\_{1} \tag{8}$$

$$\begin{cases} \ x\_b \ge 0\\ \ y\_r \ge 0\\ z\_k \ge 0 \end{cases} \tag{9}$$

Equation (4) indicates the metal quantity constraint: the actual quantity of metal produced cannot be lower than the minimum amount of metal required during the planning period.

Equation (5) indicates the ore reserve constraint: the ore reserves of each mine room and pillar, as well as the ore reserves of the mining preparation works, are limited and cannot be exceeded.

Equation (6) indicates the production capacity constraints: the development of extraction plans needs to consider the existing technical conditions of the mine's production and production equipment, etc., and cannot exceed the maximum production capacity.

Equation (7) indicates the production continuity constraints: once the mine is fully operational, it must ensure that production activities can be carried out continuously and steadily; therefore, the mine's mining plan must be prepared to consider mining preparation blocks. Back mining blocks in the ore reserves can be articulated reasonably as the general requirements of the planned period for mining preparation blocks, and the ore reserves are not less than the planned period for back mining of ore to ensure that the mine production activities can be carried out continuously and steadily.

Equation (8) indicates the mine hoisting capacity constraint: the mining capacity cannot exceed the maximum hoisting capacity of equipment during the planning period.

Equation (9) represents the variable non-negative constraint: recovery of ore quantity in all mine rooms and pillar and mining preparation block quantity are non-negative quantities that are greater than or equal to zero.

#### *2.3. Optimized Solution*

The optimal solution (*xb*, *yr*, *zk*) to the objective is the optimal operation plan for the planning period. However, due to the large number of decision variables, it is difficult to find the optimal solution by manual calculation. In order to solve the above MOP problem, this paper does not use the traditional single-objective processing methods, such as weight coefficient summation or objective constraint, but finds a set of feasible solutions based on the Pareto optimal solution set in order to provide more decision options for the decision maker. Therefore, it is necessary to find an appropriate algorithm and then solve the problem with the powerful computing facilities. In this paper, an improved particle swarm algorithm is adopted to solve the above MOP problem [31].

#### **3. Materials and Methods**

#### *3.1. Basic Particle Swarm Optimization Algorithm*

The particle swarm algorithm was proposed by Kennedy, an American psychologist, and Dr. Eberhart, an electrical engineer, in 1995 [32]. This algorithm is used to guide the updating of the velocity and position of the next generation of particles by the individual optimal position and the global optimal position of the particles in the population, so that the particles always move toward the individual optimal position and the global optimal position and finally converge to the global optimal position [33,34].

In the *M*-dimensional search space, there exists a population *xij* consisting of *N* particles, where *i* ∈ {1, . . . , *N*}, *j* ∈ {1, 2, ··· , *M*}; the velocity and position components of the *i* particle at the *t* iteration in the *j* dimension are *vij*(*t*) and *xij*(*t*), respectively; *Pbest*,*ij*(*t*) denotes the optimal position component of the particle; *gbest*,*j*(*t*) denotes the optimal position component of the population; *w* is the inertia weight coefficient; *c*1, *c*<sup>2</sup> is the learning factor; *r*1,*r*<sup>2</sup> is a random number obeying uniform distribution within (0, 1) [35].

The velocity and position update of the basic particle swarm algorithm are as indicated by Equations (10) and (11) [36].

$$w\_{\overline{i}\overline{j}}(t+1) = wv\_{\overline{i}\overline{j}}(t) + c\_1r\_1(P\_{\text{best},\overline{j}}(t) - \mathbf{x}\_{\overline{i}\overline{j}}(t)) + c\_2r\_2(G\_{\text{best},\overline{j}}(t) - \mathbf{x}\_{\overline{i}\overline{j}}(t)) \tag{10}$$

$$
\pi\_{ij}(t+1) = \pi\_{ij}(t) + v\_{ij}(t+1) \tag{11}
$$

#### *3.2. Optimization Strategies*

#### 3.2.1. Nonlinear Dynamic Decreasing Inertia Weighting Strategy

The inertia weight *w* is one of the important parameters affecting the performance of PSO algorithm. The larger the inertia weight is, the faster the particle swarm moves, the ability to search the local space is reduced, and the ability to detect the global space is enhanced; the smaller the inertia weight is, the lower the speed of the particle moves, the ability to search the local space is enhanced, and the ability to detect the global space is reduced. Therefore, appropriate improvement of inertia weights can improve the performance of the PSO algorithm. For this reason, Shi [37] proposed a linear decreasing particle swarm optimization (LDPSO) algorithm, with the following linear decreasing inertia weight strategy:

$$w = (w\_{\text{max}} - w\_{\text{min}}) \times \left(\frac{t\_{\text{max}} - t}{t\_{\text{max}}}\right) + w\_{\text{min}} \tag{12}$$

where *tmax* denotes the maximum number of iterations; *t* denotes the current number of iterations; *wmax* and *wmin* denote the set maximum inertia weights and minimum inertia weights, respectively. Through several experimental analyses, the authors found that as the number of iterations increases, when *wmax* = 0.9, *wmin* = 0.4, the search performance of the algorithm is greatly improved [38,39].

In this paper, based on the LDPSO algorithm and inspired by the ideas in the literature [39–42], we propose a nonlinear dynamic decreasing weight strategy particle swarm optimization (NDIWPSO) algorithm in which the inertia weights are set as nonlinear exponential functions, as shown in Equation (13):

$$w = (w\_{\text{max}} - w\_{\text{min}}) \times \left(1 - \frac{1}{1 + \exp((-kt) \div t\_{\text{max}})} \right) + w\_{\text{min}} \tag{13}$$

where *k* is the control factor. Compared with LDPSO, when *w* is close to *wmax*, the value of *w* increases, and NDIWPSO is dominated by the first term of Equation (10) during most of the iterations. Then, the ability of the particle to expand the search space is enhanced, and it is more advantageous in searching the global space, while the contraction ability of the particle to the location of the optimal value decreases, and the ability to search the local space is then weakened; when *w* is close to *wmin*, the value of *w* decreases and NDIWPSO is dominated by the last two terms of Equation (10) during most of the iterations, and then the particle is more advantageous in searching for the optimal value in the local interval, and the ability to search the global space is weakened. So, the nonlinear exponential function, Equation (13), can coordinate the algorithm to achieve better between the ability of local search and global search.

#### 3.2.2. Dynamic Learning Factor Strategy

In addition to the inertia weights *ω*, the learning factors *C*1, *C*<sup>2</sup> also need to be improved to make them more suitable for system optimization search.

$$c\_1 = 1 - \ln 2 \left(\frac{t}{t\_{\text{max}}}\right) \tag{14}$$

$$c\_2 = 1 + \ln 2 \left(\frac{t}{t\_{\text{max}}}\right) \tag{15}$$

Equations (14) and (15) are the improvement of the learning factor [42]; as the iterative process proceeds, *c*<sup>1</sup> decreases while *c*<sup>2</sup> increases, and in the early iterative process, the particles are mainly influenced by the individual information, which is beneficial to increase the population diversity. In the later iterative process, the particles are mainly influenced by the population information, which is beneficial to the particles to rapidly approach the global extremes and obtain the optimal solution.

3.2.3. Adaptive Variation Probability and Global Optimal Hybrid Mutation Strategy

In order to increase the diversity of the population, a variation strategy is applied to it based on the previous optimization, which leads the particles to jump out of the local optimal value points and search in a more global space. The specific implementation is as follows.

If *f <sup>t</sup>* avg, *f*max*<sup>t</sup>* , *f* min*<sup>t</sup>* represent the mean, maximum and minimum values of particle fitness in the *t* generation, respectively, then Equation (16) for the aggregation degree *δ* of particles in the *t* generation is as follows.

$$\delta = \frac{1}{N} \sum\_{i=1}^{N} \left| \frac{f\left(\mathbf{x}\_i^t\right) - f\_{\text{avg}}^t}{f\_{\text{max}^t} - f\_{\text{min}^t}} \right| \tag{16}$$

When the deviation of the fitness of individual particles from the overall average fitness is larger, the diversity of particles is better, and vice versa, the diversity of particles is worse. Therefore, when *δ* is larger, the diversity of particles is better, and when *δ* is smaller, the diversity of particles is worse. Therefore, the aggregation degree *δ* can be used to dynamically adjust the probability of variation of particles in each generation, and let the probability of variation in the *t* generation be *p<sup>t</sup> <sup>m</sup>*, as shown in Equation (17). *α* is used to regulate how fast the variance probability changes and is a constant that takes values in the range [2,4].

$$
\psi\_m^t = \delta e^{-a\left(1 + \frac{t}{l\_{\text{max}}}\right)} \tag{17}
$$

In this paper, a mutation on the optimal position of the particle is used, and if *rand* ∈ " 0, *p<sup>t</sup> m* # , a mixture of Gaussian (Equation (18)) and Cauchy (Equation (19)) distributions is used to vary the optimal position of the particle. *pg* denotes a random one of the global optimal and global suboptimal positions; *gbest* denotes the global optimal position; *randn* is a random number of Gaussian distribution; *Cauchy* is a random number of Cauchy distribution.

$$g\_{best} = p\text{g} \times (1 + 0.5 \times randn) \tag{18}$$

$$\text{Cauchy} = \tan(\pi \times (rand - 0.5))\tag{19}$$

$$\lg\_{\text{best}} = p\lg \times (1 + 0.5 \times \text{Cauchy}) \tag{20}$$

The pseudo-code for the global optimal hybrid mutation strategy is **Mut1** (Algorithm 1).

#### **Algorithm 1: Mut1.**


#### 3.2.4. Worst Personal-Best Position Adaptive Wavelet Mutation Strategy

In the above hybrid variation strategy, the global optimal and suboptimal extremes are utilized for variation, while the worse individual extremes are not utilized. In fact, in the population, the worse individual extremes have limited guiding effect on the particles; thus, variation on the worse individual extremes is beneficial to accelerate the convergence of the population.

The worst individual optimal adaptive wavelet variation strategy [43] (**Mut2**):

The worst individual extreme value particle *m* is selected, and the search boundary of the selected particle in the *j* dimension is *Pbest*,*mj* ∈ [*x*min, *x*max]. *Pbest*,*mj* is varied according to Equation (21), and in this paper, adaptive wavelet variation is used to improve the worst individual extreme value to speed up the evolution of the population. *δ* is the wavelet function value, *ψ*(*x*) is the wavelet function, and Morlet wavelet is chosen as the wavelet base in Equations (22)–(24). *a* is the scale parameter, and more than 99% of the energy of the wavelet function is contained in (−2.5, 2.5), so the range of values of *ϕ* in the formula is the pseudo-random number of (−2.5a, 2.5a) [44]. The wavelet amplitude decreases continuously with the increase of parameter *ψ*(*x*). In order to adjust the wavelet amplitude *ψ*(*x*) adaptively, the adaptive parameter *a* is proposed in this paper, and its expression is Equation (25). *k* and *a*<sup>0</sup> are positive constants, and this paper sets *k* = 10, *a*<sup>0</sup> = 5.

$$P\_{\text{best},\text{mj}} = \begin{cases} P\_{\text{best},\text{mj}} + \sigma \left( \mathbf{x}\_{\text{max}} - P\_{\text{best},\text{mj}} \right), \text{if } \sigma > 0 \\\ P\_{\text{best},\text{mj}} + \sigma \left( P\_{\text{best},\text{mj}} - \mathbf{x}\_{\text{min}} \right), \text{if } \sigma \le 0 \end{cases} \tag{21}$$

$$
\sigma = \frac{1}{\sqrt{a}} \psi \left( \frac{\varrho}{a} \right) \tag{22}
$$

$$\psi(x) = e^{-\frac{\mathbf{x}^2}{2}} \cos(5x) \tag{23}$$

$$
\sigma = \frac{1}{\sqrt{a}} e^{-\frac{\left(\frac{\rho}{2}\right)^2}{2}} \cos\left(5\left(\frac{\rho}{a}\right)\right), \rho \in [-2.5a, 2.5a] \tag{24}
$$

$$a = kt + a\_0 \tag{25}$$

#### *3.3. Algorithm Steps*


#### *3.4. Construction of the Fitness Evaluation Function*

For the MOP problem model of the extraction plan, this paper uses the evaluation function method to construct the fitness function. Specifically, we set an arbitrary best value *Q*<sup>∗</sup> and *A*<sup>∗</sup> for the objective functions *Q*(*xi*, *yj*, *zk*) and *A*(*xi*, *yj*, *zk*) respectively, and construct the fitness evaluation function *F*(*xi*, *yj*, *zk*) of the extraction plan by calculating the difference between each objective function and it and pursuing the minimum of the sum. Equation (26). *A*<sup>∗</sup> is the planned quantity of ore to be produced, *Q*<sup>∗</sup> is the planned economic efficiency of the enterprise during the plan period.

$$F(\mathbf{x}\_{i\prime}y\_j, z\_k) = \sqrt{\left(\frac{\mathbf{Q} - \mathbf{Q}^\*}{\mathbf{Q}^\*}\right)^2} + \sqrt{\left(\frac{A - A^\*}{A^\*}\right)^2} + \sqrt{\left(\frac{\sigma\_l - \sigma\_l^\*}{\sigma\_l^\*}\right)^2} \tag{26}$$

#### **4. Engineering Applications and Results Analysis**

#### *4.1. Mine Engineering Overview and Basic Data*

In order to verify the effectiveness of the NAMPSO algorithm for industrial mining planning of multi-metal mines, the actual mining production data of a large underground metal mine rich in gold (Au), antimony (Sb), and tungsten (WO3) in Hunan Province, China, are taken as an example, and the numerical model of the mine is shown in Figure 2. In this paper, the annual production task index of this mine and the actual situation of the mine are combined, and the balanced mining extraction plan model established in Section 2 and the NAMPSO algorithm for solving multi-objective optimization problems proposed in Section 3 are applied to optimize the extraction plan of this mine for each quarter and obtain the annual optimal mining plan.

**Figure 2.** Digital model of underground multi-metal mine.

The production target of the mining company for this mine in the current year is shown in Table 1.

**Table 1.** Production task indicators.


Among them, after the completion of the mining tasks of the previous year, the ore reserves of the mine's preparation rooms (13) and preparation pillars (8), as well as the part of the works of the excavated blocks (9) that need to be completed in this year, are shown in Table 2.


**Table 2.** Parameters of the mine rooms, mine pillars, and excavated blocks.

#### *4.2. Balanced Mining Model Parameters*

According to the balanced mining extraction plan model constructed in Section 2.2 of this paper, the parameters involved in the mathematical model of the extraction plan are organized as follows, considering the actual situation of the mine and referring to the production experience of the mine in previous years:


#### *4.3. Algorithm Parameter Setting*

For the above engineering example data and the constructed model, the algorithm is applied on the MATLAB platform to solve the plan model optimally. There are 30 parameters of decision variables in this metal mine mining planning model, so the dimension of each particle in the algorithm *M* = 30; the population size of particles *N* = 150, the number of iterations *T* = *t*max = 3000; the learning factor *c*1, *c*<sup>2</sup> varies with the number of iterations *t*; the ore loss rate of the mining site is considered to be 5%.

#### *4.4. Engineering Example Simulation*

Based on the equilibrium mining model and the data in Tables 1 and 2, the PSO algorithm and the NAMPSO algorithm were applied to find the optimal extraction plan for each quarter of a year, and the iterative convergence process of the model solution was calculated to obtain four quarters, as shown in Figure 3. From the fitness function curves, we can see that the fitness function values of the two algorithms are close to the improved algorithm slightly better, but the convergence speed of the improved PSO algorithm is significantly faster in the first two quarters (PSO converges around the 1550th, 450th, 1400th, and 120th generations, and NAMPSO converges around the 1100th, 250th, 1200th, and 80th

generations), the convergence speed is improved by about 29.25%, and the specific results of the balanced extraction plan are shown in Tables 3–5.

**Figure 3.** Convergence curve of the optimal adaptation degree of the mining plan.


**Table 3.** Quarterly and annual mining plan.

**Table 4.** Quarterly and annual excavation plan.


From the convergence curve in Figure 3, it can be seen that when the NAMPSO algorithm is used to solve the extraction plan, the value of the objective function fluctuates between 3.2 and 1.2 when the number of iterations is between 0 and 500, indicating that the algorithm converges to the global optimal solution slowly at the early stage of iterative computation when the number of particles within the population is high under a certain initial population size and at the late stage of iterative computation, i.e., after 500 iterations, the algorithm starts to converge smoothly and rapidly to the global optimal solution of the objective function. It can also be seen that the NAMPSO algorithm, when solving the production planning model with complex constraints, performs nonlinear dynamic optimization of the inertia weights of the algorithm and introduces an adaptive variational probability strategy to make the particles outside the feasible domain enter the feasible domain quickly, which leads to a slow convergence computation at the early stage of the algorithm search, while the computation speed of the algorithm is faster at the later stage, and it is easy to jump out of the local optimal solution problem.


**Table 5.** Extraction plan optimization results by quarter and year.

#### *4.5. Analysis of Results*

Combined with the enterprise's target requirements for the mine's extraction production, a comparative analysis of the extraction plans given in this paper for each quarter shows that:

(1) As shown in Figure 4, in terms of mining quantity, the actual mining quantity of each quarter and year in the balanced mining plan is slightly higher than the planned mining quantity, which can meet the annual mining quantity target specified by the enterprise. At the same time, the proportion of mining quantity in each quarter is 26%, 25%, 24%, and 25%, respectively, which meets the requirement of balanced mining.

**Figure 4.** Convergence curve of the optimal adaptation degree of the mining plan.

(2) As shown in Figure 5, in terms of the production of gold, antimony, and tungsten, the actual production for each quarter and year in the extraction plan given in this paper is slightly higher or equal to the planned production, which fully meets the annual metal production target set by the enterprise. The percentages of gold production in each quarter are 24%, 23%, 24%, and 29%, respectively; the percentages of antimony production in each quarter are 26%, 24%, 23%, and 27%, respectively; and the percentages of tungsten production in each quarter are 25%, 25%, 25%, and 25%, respectively. From the absolute uniform distribution of tungsten production in each quarter, it can be seen that the algorithm takes the achievement of tungsten production as one of the key conditions when searching for the optimal mining plan, which fully meets the metal content constraint requirement of the model, and because the grade of tungsten is relatively the lowest among the three metals in the ore of each mining site, the value of tungsten is relatively the lowest among the three metals, and the recovery rate of tungsten beneficiation is also the lowest.

**Figure 5.** Comparison of actual and planned production of metal quantities.

(3) As shown in Figure 6, the average grade of gold and tungsten is higher than the planned ore grade, the average grade of gold and tungsten in the ore mined in each quarter is basically balanced, and the average grade of antimony in the ore mined in each quarter is basically the same as the planned grade. Therefore, the mining plan searched by the method of this paper achieves the balance of the ore grade, and the ore grade fluctuation is small.

**Figure 6.** Comparison of actual grade and planned grade of ore mined.

(4) As shown in Figure 7, in terms of economic benefits, the annual economic benefits of the extraction plan given in this paper increased by 41.88% compared to the planned economic benefits. The actual production of each quarter and year is slightly higher than the planned production, and the production of each metal also meets the requirements, so the economic benefits of each quarter and year are necessarily higher than the planned economic benefits, and the actual benefits of each quarter account for 24%, 23%, 24%, and 29%, respectively, and in general, the economic benefits of each quarter are balanced and stable, which can well fulfill the benefit targets given by the company.

**Figure 7.** Comparison of actual and planned economic benefits.

(5) As shown in Figure 8, in terms of mining balance, in order to seek the quarterly quantity of ore production that satisfies each constraint, the quantity of ore mined from each mineral deposit (including the mine room and pillar) varies each quarter, which reflects the equilibrium process of ore matching. Through four quarters of mining, the ore reserves of each mining field are basically depleted, and in actual production the ore reserves of the mining field can be considered as the end of mining when they are below a certain value. After one year of mining, most of the 21 recovery mining fields involved in the plan can be considered as completed. In order to ensure that mine production can be carried out continuously, it is necessary to adhere to the principle of balanced mining and excavation with excavation first. In this paper, the mining plan is formulated mainly through the balance between the recovery quantity and mining preparation quantity to ensure the continuity of production, as shown in Figure 9. This paper contains nine mining preparation fields, each of which has a different footage of excavation in each quarter, and the reason for this is that the

quantity of ore production from each mining preparation project also must meet the production constraints. After four quarters of mining preparation work, the nine mining preparation fields are basically finished and can be used as back mining fields for the next year, which can ensure the continuous progress of production work in the next year.

**Figure 8.** Distribution of ore production and remaining ore in each mining field by quarter.

#### **5. Conclusions**

Aiming at the multi-objective optimization problem faced by underground multimetal mine extraction plans, this paper proposes an improved particle swarm optimization algorithm, which makes the global relative optimal solution converge smoothly and quickly by nonlinear dynamic optimization of inertia weights of a particle swarm algorithm, while introducing an adaptive mutation probability strategy, mixed mutation strategy for relatively optimal individuals, and adaptive wavelet mutation strategy for the poorer optimal individuals to achieve the solution of the balanced mining plan. The optimal extraction plan searched was compared and analyzed with the production target proposed by the enterprise in five aspects: economic efficiency, ore production, metal quantity, ore grade, and extraction balance. The results prove that the balanced mining plan model constructed by using this paper and the quarterly and annual mining plans obtained by using the NAMPSO algorithm not only fully meet the production task requirements, but also the algorithm solution results are 10.98% higher than the mine plan index in terms of ore quantity, 41.88% higher in terms of economic efficiency, and 29.25% higher in terms of algorithm solution speed, which can well achieve the balanced production of the mine. Thus, the practicality and feasibility of the model and algorithm are highlighted.

The optimization of the engineering examples in this paper focuses on optimization with deterministic information and has limitations that can provide a reference for decision makers. Future research will focus on optimization under uncertainty, seeking quantitative methods for the uncertain information in the complex system of a mine. At the same time, the optimization algorithm solution can be adapted in the future not only in the mine quarry but in the whole mine system, as it is included in the optimization model, while updating the optimization algorithm to achieve the optimization of the mine system from a single link to the whole, thus helping the mining enterprises to obtain greater economic benefits and improve the efficient use of limited mineral resources.

**Author Contributions:** Conceptualization, methodology, formal analysis, and writing—review and editing, Y.Z.; software, investigation, and visualization, Y.C.; writing—original draft preparation and validation, Y.Z. and Y.C.; resources and data curation, S.Y.; supervision, project administration, and funding acquisition, J.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Natural Science Foundation Project of China, grant no. 72088101, no. 5140430, and no. 551374242, and the Graduated Students' Research and Innovation Fund Project of Central South University, grant no. 2021zzts0283.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** The authors would like to express thanks to the National Natural Science Foundation and Innovation Fund Project of Central South University.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**



#### **References**

