Enumeration optimization of open pit production scheduling based on mobile capacity search domain
The production capacity, mining sequence and mining lifetime are the three significant factors in the optimization of open pit mine production scheduling. While production capacity is related to technical conditions of the mine, mining sequence is related to mine organization management, and mining lifetime is restricted by mine scheduling and industry norms; furthermore, these three factors are affected by geological conditions of the ore body, reserves scale, marketing environment and other elements^{1,2,3,4}. In addition, compared with other external conditions, there is also interaction and influence among the three elements of production scheduling. For example, the more significant the production capacity is, the shorter the mining productive life. Different mining sequences (pit advancement position at the end of the year) directly affect the year’s production capacity and then affect the overall mining life. In practice, therefore, optimized mining scheduling determines the annual estimated amount of mined ore, the annual estimated amount of stripped waste rock (i.e., production capacity), which areas are mined and stripped each year or how each step is advanced (i.e., the mining sequence), and the mining period (i.e., the mining lifetime).
Production scheduling provides a future production strategy for a mine (new mine or mine in producing). Moreover, for a given deposit, the quality of the production schedule has a significant impact on infrastructure investment and cash flow, which is distributed on the timeline after production. This influence becomes an important factor in the investment return rate of the whole mine project. That is why international mining companies have a strong interest in optimizing production scheduling. Production scheduling optimization is always a popular research topic in mine system engineering. From the perspective of optimization, the open pit mine production schedule determines the mining time of each module in the massive deposit model to determine which module should be mined each year to maximize the total NET present value and meet the space–time relationship and technical and economic constraints of open pit mining.
The “mining increment sorting method” is the earliest computer optimization method applied to production scheduling. It was first proposed by engineers of the Kennecott Company and applied in the company. The mining increment is generated by constructing cones, and cone structure, evaluation and sequencing are carried out via human–computer interaction trial and error^{5}. Production increments for production scheduling can be obtained by the “parametric analysis” proposed by Lerchs and Grossmann^{6}. This method is further developed into the “reserve parameterization” method, and many scholars have carried out further research on the solution of reserve parameterization and its application in production scheduling^{7,8,9,10,11,12}.
One of the inherent defects of parameterization is the “notch” problem; in the generated limit sequence, the increments between some of the adjacent limits are so large that the limit sequence cannot be used for production scheduling optimization. Therefore, some researchers use heuristic algorithms to generate nested limit sequences to overcome the gap problem^{13,14,15,16,17,18,19,20} and dynamically order the resulting limit sequence^{18,21,22,23,24,25}.
In summary, the essence of the production scheduling optimization problem is to determine the optimal production time of each module on the premise of meeting the necessary constraint conditions to obtain the maximum total NET present value. It is a typical linear programming problem. Therefore, linear programming (its specific form includes mixed programming, pure integer programming and 0–1 programming) is one of the most commonly used mathematical optimization methods to solve the optimization problem of production scheduling; related studies were carried out as early as the late 1960s ^{26,27}. Many researchers have established linear programming models with different concrete forms for different aspects of production scheduling problems^{27,28,29,30,31,32,33,34,35}.
However, the quantity of variables and constraint equations in the linear programming model for optimizing production scheduling is too enormous when a single module in the massive deposit model is used as the decision unit. That is a situation that even today’s computers cannot solve directly; if it is integer programming, it can hardly be solved. Therefore, some researchers try to find solutions in the construction of mathematical model forms (mainly constraints) or solving algorithms (usually by feat of approximate algorithms) to boost the solving speed^{36,37,38,39,40,41,42}. Increasing the decisionmaking unit to reduce the quantity of variables and constraints is a more common way, such as combining the modules in the deposit model into a “unit tree” as the decision unit in optimization or taking steps or panels as the decision unit^{37,43,44,45,46,47,48,49,50,51}. However, due to the low scheduling accuracy (or resolution) in enormous decision units, the results are significantly different from the optimal plans, which also reduces the practicality of the results^{47,49,52,53}. Therefore, many researchers make use of the unique structure of the mathematical model to reduce the model size with the Lagrange relaxation method and solve the model with other measures and algorithms, such as iteration, decomposition, gradient method, and the Dantzig network flow method^{54,55,56,57,58,59,60,61,62,63,64,65,66,67,68}. The biggest obstacle to this approach is the “gap” problem. Researchers have tried various methods to solve this problem, but they have not sought out suitable means.
Open pit mine production scheduling is a typical multiperiod decisionmaking problem; moreover, the time is interrelated so that dynamic programming can be used to solve it. Therefore, dynamic programming is also one of the most used mathematical optimization methods to solve the problem. Researchers have used dynamic programming to study the optimization of different aspects of production scheduling^{1,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84}.
Because of the insurmountable difficulties in applying mathematical optimization models to obtain the exact solution of production scheduling, some researchers turn to approximate algorithms, such as genetic algorithms, random, local search, particle swarm algorithms, and simulations, to obtain one or more “perfect” plans^{85,86,87,88,89,90,91,92,93,94,95,96,97}.
In the process of production schedule optimization, in addition to exploring mathematical methods to solve the problem, another factor is to consider the specific influencing factors of the optimization process. In the study of the traditional optimization problem, first, it is assumed that a known factor could affect the optimization of the production schedule. For example, most studies determine the mine production capacity according to the Taloy formula, similar mine factors, or relevant experience and then optimize the mining sequence based on the determined production capacity. Finally, the relevant infrastructure is allocated according to the production scheduling^{98,99,100,101} set boundary grades in different periods with the maximum net present value as the objective function to ensure that the infrastructure meets the maximum comprehensive processing capacity of mining, beatification and smelting. This method can obtain higher economic benefits than the Lane method; nonetheless, the characteristic of uniform grade distribution is not applicable to most mines. On the side, the method described by Khan and Asad et al. focuses on the impact of geological uncertainty, marketing environment change and risk management on production scheduling^{83,101,102,103,104,105,106,107,108}.
Consequently, due to the high complexity of production scheduling optimization and the difficulty of finding an exact solution, it is still a popular research topic. Many problems in the process of production scheduling optimization of open pit mines exist, such as low operation efficiency caused by large amount of calculation, disorderly expansion of production capacity caused by the scale effect, and difficulties in simultaneously optimizing the three elements due to the interaction relationship. In view of these problems, taking the maximum comprehensive net present value (NPV) as the objective function, the method of mobile capacity search domain is proposed to improve the operation efficiency, the infrastructure investment function based on the maximum production capacity is constructed to restrict the production capacity, the facility idle threshold is designed to reduce the fluctuation range of production capacity, and the enumeration method is used to evaluate all workable paths to realize the simultaneous optimization of the production capacity, mining sequence and mining life of open pit mines.
Mathematical model of discrete body dynamic programming
The purpose of production scheduling optimization of open pit mines is to determine how much ore and rock are extracted and stripped every year, where they are advanced and how long the mining life is. Therefore, the deposit should be divided into a finite number of discrete bodies, which should be taken as the alternative objects of open pit mining at the end of the year. Due to the characteristics of open pit mining technology and safety requirements, mining engineering needs to be carried out gradually from top to bottom in accordance with the specified slope angle, and the division of discrete bodies needs to meet the geometric space relationship. According to the definition and generation principle of the “geological optimal pit” proposed by^{1}, a series of fully nested discrete bodies, namely, the “geological optimal pit sequence”, can be obtained. The specific process for generating the “geological optimal pit sequence” is as follows:
Started from the bottom of the ultimate pit (also the first geological optimal pit), floating cone exclusion method is used to find out all the parts whose ore volume is equal to the given increment (such as 2 million tons), and among these parts, the one that holding the lowest average grade removed from the ultimate pit, leaving a pit that is 2 million tons less than the ultimate pit. The pit has the highest metal content of all the pits, with 2 million tons less ores than the ultimate pit. So the second geological optimal pit (the first is the ultimate pit) is obtained. Then, in the remaining pit (the second geological optimal pit), remove the part with the ore volume equal to the given increment (2 million tons) and the minimum average grade by the same method, and a smaller pit (the third geological optimal pit) is obtained. According to the method analogy, the ultimate pit is discretized into a series of completely nested discrete body sequences, that is to say, the geological optimal pit sequence.
The discrete bodies are put into the dynamic sorting model as state variables, as shown in Fig. 1.
The circles in Fig. 1 represent a discrete body (geological optimal pit), which is also where the pit is advancing at the end of a year. The two circles linked by arrows, such as J_{1} in the first year and J_{3} in the second year, linked by dotted arrows, indicate that the open pit advances to J_{1} in the first year and J_{3} in the second year. In geometric space, J_{3} completely contains J_{1}, therefore the quantity of ore corresponding to J_{1} and J_{3} can be expressed as O_{1},_{1} and O_{2},_{3}, and the quantity of rock can be expressed as R_{1,1} and R_{2,3}. In fact, since the quantities of ore and rock of each discrete body have been determined, O_{1,1} and O_{2,3} are equal to O_{1} and O_{3}; R_{1,1} and R_{2,3} are the same as R_{1} and R_{3}, mainly for the purpose of relating to time. Then, the quantity of ore and rock of the jth discrete body corresponding to the ith year can be expressed as O_{i,j} and R_{i,j,} respectively. Therefore, the quantity of ore o_{i−1, k} (i, m) and rock produced from discrete body J_{k} in year i − 1 to discrete body J_{m} in year i where m > k, because the discrete body must be advanced from small to large, are r_{i−1,k} (i, m), respectively:
$$o_i – 1,k \left( i,m \right) = \, (\textO_i – 1,k – \textO_i,m ) \times \gamma /\left( 1 – \delta \right)$$
(1)
$$r_i – 1,k \left( i,m \right) = \, \left( \textR_i – 1,k – \text R_i,m \right) – o_i – \text1,k \left( i,m \right) \, \times \delta + \, \left( \textO_i – 1,k – \text O_i,m \right) \times \left( 1 – \gamma \right)$$
(2)
where γ is the ore recovery rate and δ is the mixing rate of waste rock.
The concentrate amount v_{i1},_{k}(i, m), which transitions from the discrete body J_{k} in year i − 1 to the discrete body J_{m} in year i, is:
$$v_i – 1,k \left( i,m \right) = \left[ \left( O_i – 1,k \times g_k – \, O_i,m \times g_m \right) \times \gamma + o_i – 1,k \left( i,m \right) \, \times \delta \times g_r \right]/g_v$$
(3)
where g_{k} and g_{m} represent the average grade of ore in the kth and mth discrete bodies; g_{r} represents the average grade of rock (R_{i1,k} − R_{i,m}) transiting from discrete J_{k} to discrete J_{m}; g_{v} refers to concentrate grade.
The profit u_{i − 1},_{k}(i, m), which transitions from the discrete body J_{k} in year i − 1 to the discrete body J_{m} in year i, is:
$$u_i – 1,k \left( i,m \right) = v_i – 1,k \left( i,m \right) \, \times p_v \times \left( 1 + \eta \right)^i – 1 – \left[ o_i – 1,k \left( i,m \right) \, \times \left( c_o + c_p \right) \, \times \left( 1 + \varepsilon \right)^i – 1 + r_i – 1,k \left( i,m \right) \, \times c_r \times \left( 1 + \varepsilon \right)^i – 1 \right]$$
(4)
where p_{v} refers to the price of concentrate, ¥/t; c_{o}, c_{p} and c_{r} represent mining cost, processing cost and stripping cost, respectively, ¥/t; and η and ε represent the price rise rate and cost rise rate, respectively.
The net present value NPV_{i1},_{k} (i, m), which transitions from the discrete body J_{k} in year i − 1 to the discrete body J_{m} in year i, is:
$$\textNPV_i – 1,k \left( i,m \right) \, = u_i – 1,k \left( i,m \right)/\left( 1 + \lambda \right)^i$$
(5)
where λ represents the discount rate.
The transition process mentioned above is the transition process between different states. Without considering the constraints of production capacity, any year in Fig. 1 may be the final life of the mine. It is also the year in which the mine is mined to the final state J_{N} (discrete body, which is also the ultimate pit of the open pit mine). In addition to the three production costs of mining, stripping and processing given in Eq. (4), another significant cost is infrastructure investment. It generally occurs at the initial stage of mine construction, i.e., the infrastructure investment does not affect the state transition shown in Fig. 1. However, it will affect the final economic benefit of the mine. Therefore, it is necessary to consider the infrastructure investment in each NPV that reaches the final state J_{N} through different paths. Although infrastructure investment has no influence on a state transition, infrastructure investment is affected by mine production capacity, and mine infrastructure investment usually needs to meet the maximum production capacity in the whole life cycle of the mine. In other words, the infrastructure investment needs to meet the maximum ore quantity difference between any two adjacent time points in the state transition process on a path to the final state J_{N}, as shown in Fig. 1; that is, the maximum production capacity on this path is maxo_{i1,k}(i, m). According to Eq. (1), the quantity of ore produced by state transition at two adjacent time points is o_{i1},_{k}(i, m). Therefore for a particular path L(D) with a mining life of D years, its maximum production capacity q_{L(D)} is:
$$q_L\left( D \right) = \mathop max\limits_k \in \left[ i – 1,m – 1 \right] \left\{ { \left( \textO_i – 1,k – \text O_i,m \right) \times \upgamma /\left( 1 – \updelta \right) } \right\}$$
(6)
The infrastructure investment of the mine can be approximated as a linear function of the maximum production capacity of the ore. For the path L(D) with a mining life of D years, the infrastructure investment c_{L(D)} can be approximated as:
$$c_\textL(\textD) = \texta + \textb \times q_{\textL(\textD)}$$
(7)
where a refers to the infrastructure investment base unrelated to the production scale, 10^{4} ¥, and b refers to the infrastructure investment per unit of mining amount, ¥/t.
The comprehensive economic benefit NPV_{L(D)} for path L(D) can be expressed by the following formula:
$$\text NPV_\textL\left( \textD \right) = c_{\textL\left( \textD \right)} + \mathop \sum \limits_i = 1^D \textNPV_i – 1,k \left( i, m \right).$$
(8)
Enumeration optimization algorithm of production scheduling
Set a wide enough range of ore production capacity [q_{low}, q_{up}], and the optimal production capacity must be within this range. This range is set by the user and is the input data. Let n represent the ordinal number of the constraint domain of ore production capacity, and [q_{L}^{n}, q_{U}^{n}] define the nth constraint domain (capacity search domain).

Step one: Set n = 1 and juxtapose the Bohr variable LastPlan = false. In the discrete body sequence J_{N}, the discrete body whose ore quantity is greater than or equal to and closest to q_{low} is found, and its number is denoted as H.

Step two: Set the lower bounds q_{L}^{n} and upper bounds q_{U}^{n} of the constraint domain [q_{L}^{n}, q_{U}^{n}]:
$$q_\textL^n = \textO_\textH \Delta \textP/2$$
(9)
where ΔP is the ore increment of the pit set when the discrete body sequence J_{N} is generated^{1}, and O_{H} is the ore quantity of the Hth discrete body in the sequence J_{N}.
(10)
where O_{N} is the total ore quantity in the ultimate pit and x_{1}, x_{2}, x_{3} are the search domain constraints, which are determined according to the size of the mine. If G > N obtained by Eq. (10), set G = N. The upper bound of the domain is
$$q_\textU^n = \textO_\textG + \Delta \textP/2$$
(11)
where O_{G} is the ore quantity of the Gth discrete body in sequence J_{N}.
If q_{U}^{n} > q_{up}, LastPlan = true.

Step three: Find all production scheduling whose annual ore output meets the constraints of the domain [q_{L}^{n}, q_{U}^{n}]. That is, find all subsequences in which the annual ore output falls into the domain [q_{L}^{n}, q_{U}^{n}]. Calculate the NPV of each schedule, record the schedule with the highest NPV, and call the schedule the “local best schedule” in the domain [q_{L}^{n}, q_{U}^{n}]. The specific method is as follows:

Step 1: Set i = 1 (first year). A discrete body is found in the discrete body sequence J_{N}. The ore quantity is no less than and closest to q_{L}^{n}, the lower limit of the constraint domain [q_{L}^{n}, q_{U}^{n}] of the set ore production capacity, and the total quantity of ore and rock is no more than q_{up} of the set annual mining and stripping capacity. If such a discrete body is found, its number in J_{N} is H(1); that is, the pit at the end of the first year on the planned path L under construction is advanced to pit J_{H(1)}*. According to Eqs. (1), (3), calculate the quantities of produced ore o_{1}, concentrate v_{1} and stripped waste rock r_{1} in this year. Then, calculate the annual profit u_{1} and its net present value NPV_{1} discounted to time 0 according to Eqs. (4), (5). Move on to the next step. If such a discrete body is not found, there is no workable schedule, and the algorithm terminates.

Step 2: set i = i + 1 (next year).

Step 3: The discrete body number H(i) is equal to H(i − 1) + 1 in year i, and H(i1) is the discrete body number of the previous year of the planned path L under construction.

Step 4: According to Eqs. (1), (2), calculate the ore mining quantity o_{i} and waste rock stripping quantity r_{i} from i − 1 to i.

Step 5:

(A)
If o_{i} < q_{L}^{n} and o_{i} + r_{i} ≤ q_{up}.

(B)
H(i) = N, that is, the discrete body J_{H(i)} is the last (ultimate pit) in the sequence J_{N}, then the planned path L under construction has reached the end point, and the final state J_{N} is the discrete body of the end point on the path. A complete workable planned path is obtained, and its mining life D = i. According to Eq. (3), the concentrate quantity v_{i} in the first year is calculated, and then according to Eq. (4–5), the annual profit u_{i} and its net present value NPV_{i} discounted to time 0 are calculated. Go to step 6.

(C)
H (i) < N, the ore quantity of year i is lower than the lower limit of the set annual ore production capacity, which is not workable. Set the discrete body number H(i) to H(i) + 1, that is, consider a larger discrete body, and return to step 4.

(D)
q_{L}^{n} ≤ o_{i}≤ q_{U}^{n}and o_{i} + r_{i} ≤ q_{up}.
Discrete body J_{H(i)} is the workable pit state at the end of i on the planned path L under construction. According to Eq. (3), the concentrated quantity v_{i} of year i was calculated. Then, according to Eqs. (4), (5), the annual profit u_{i} and its net present value NPV_{i}discounted to time 0 were calculated. If H(i) = N, then the planned path under construction has reached the final state, and a complete workable planned path has been obtained. Set the mining life of the intended path to D = i and carry out step 6. Otherwise, return to step 2.

(E)
If o_{i} > q_{U}^{n} or o_{i} + r_{i} > q_{up}.
The
quantity of ore mining or the total quantity of mining and stripping in year i exceeds the set upper limit, there is no workable plan, and the algorithm terminates.

(A)

Step 6: Thus far, a workable planned path with “minimum ore output” is obtained, that is, the quantity of ore extracted in each year of the path except the last year is just enough to meet the set minimum annual ore production capacity q_{L}^{n}. Calculate the infrastructure investment c _{L(D)} of the shcedule according to Eqs. (6), (7). Calculate the total NPV _{L(D)} of this path according to Eq. (8). Take the path as the current path and save it as the best path.

Step 7: Set time i = D − 1, where D is the mining life of the current path L.

Step 8: Build a new workable plan path starting from year i. The new path will be the same as the current path in year 1 to (i − 1). Add 1 to the current path’s discrete body number of year i, that is, set H(i) = H(i) + 1.

Step 9: According to Eqs. (1), (2), calculate the ore quantity o_{i} and waste rock quantity r_{i} from i − 1 to i.

Step 10: Distinguish feasibility and economic evaluation, divided into two cases:

(A)
If o_{i}< q_{U}^{n} and o_{i} + r_{i} ≤ q_{up}.
Discrete body J_{H(i)} is the workable pit state at the end of year i of the new planned path being constructed, and it becomes the discrete body of last year i of the current path (that is, the original discrete body is replaced). According to Eq. (3), the concentrated quantity v_{i} of year i is calculated. Then, according to Eqs. (4), (5), the annual profit u_{i} and its net present value NPV_{i}discounted to time 0 are calculated. If H(i) = N, the new planned path under construction has reached the final state and a complete workable planned path has been obtained, and its mining life is D = I; go to step 15. Otherwise, go to step 11.

(B)
If o_{i} > q_{U}^{n} or o_{i} + r_{i} > q_{up}.
If the quantity of ore mining or the total quantity of mining and stripping in year i exceeds the upper limit set, it is not workable. Set i = i − 1, that is, go back one year along the current path. If i > 0, go back to step 8. Otherwise, all workable plan paths are constructed and evaluated; proceed to step 16.

(A)

Step 11: Set i = i = i + 1.

Step 12: The discrete body number H(i) is equal to H(i − 1) + 1 in year i, where H(i1) is the discrete body number of the previous year of the new planned path being constructed.

Step 13: According to Eqs. (1), (2), calculate the ore quantity o_{i} and waste rock stripping quantity r_{i} from i − 1 to i.

Step 14:

(A)
if o_{i}< q_{L}^{n} and o_{i} + r_{i} ≤ q_{up}, in two cases:

(B)
H(i) = N, that is, the discrete body J_{H(i)} is the last (ultimate pit) in the sequence J_{N}, the new planned path under construction has reached its endpoint, and the ultimate pit J_{N} is the end discrete body on this path. Finally, a complete workable planned path is obtained, and its mining life D = i. According to Eq. (3), the concentrated quantity v_{i} of year i was calculated, and then according to Eqs. (4), (5), the annual profit u_{i} and its net present value NPV_{i} discounted to time 0 were calculated. Perform step 15.

(C)
H(i) < N, the ore quantity of year i is below the lower limit of the set ore production capacity, which is not workable. Set discrete body number H(i) = H(i) + 1, that is, consider a larger mining discrete body, and return to step 13.
If q_{L}^{n} ≤ o_{i}≤ q_{U}^{n} and o_{i} + r_{i} ≤ q_{up x}.
Discrete body J_{H(i)} is the workable state at the end of year i on the new plan path under construction. According to Eq. (3), the concentrated
quantity vi of year i was calculated, and then according to Eqs. (4), (5), the annual profit u_{i} and its net present value NPV_{i} discounted to time 0 were calculated. If H(i) = N, then the new planned path under construction has reached the final state, and a complete workable planned path has been obtained. Its mining life is D = i, and step 15 is carried out. Otherwise, return to step 11. 
(D)
If o_{i} > q_{U}^{n} or o_{i} + r_{i} > q_{up}.
The quantity of ore mining or the total quantity of mining and stripping in year i exceed the set upper limit, which is not workable. The algorithm gives up halfway and displays error information. Go to step 16 and output the best path thus far.

(A)

Step 15: After a new workable scheduling path is established, calculate the infrastructure investment c_{L(D)} of the schedule according to Eqs. (6), (7). Calculate the total NPV _{L(D)} of this path according to Eq. (8). If the total NPV _{L(D)} of this path is greater than the total NPV _{L(D)} of the saved best path, the path is reserved as the best path (that is, the original best path is replaced); otherwise, the original optimal path remains unchanged. Take this path as the current path and return to step 7.

Step 16: Output the bestplanned path on the interval [q_{L}^{n}, q_{U}^{n}].


Step four: If LastPlan = false, go to the next step. Otherwise, if LastPlan = true, go to step six.

Step five: Set n = n + 1 and H = H + 1. Return to step two and set the next constraint domain, and continue the iteration. In this iteration process, every time H increases by 1, the lower bound q_{L}^{n} of the domain calculated according to Eq. (9) increases once. The whole domain also moves to a higher production capacity.

Step six: The constraint domain covers the entire range of ore production capacity [q_{low}, q_{up}]. The best schedule is obtained by finding the one with the highest NPV among all the local best schedules recorded; Output the best schedule. The algorithm can also output all local best schedules to see how NPV varies with ore production capacity. The algorithm ends.
In the above algorithm, the larger the width of the constraint domain (q_{L}^{n}, q_{U}^{n}) is, the more schedules satisfy the constraint of (q_{L}^{n}, q_{U}^{n}). If the domain is set wide, the number of schedules to be evaluated in step three is too large and timeconsuming. If the range is too narrow, the probability of missing the optimal schedule is high. The conditions listed in Eq. (10) in step two are to control the width of the constraint domain within such a range. In addition, for a given domain width, the number of schedules s satisfying the domain constraint increases rapidly with increasing mining life. The detailed algorithm flow is shown in Fig. 2.
Among the workable scheduled paths, the difference between some paths’ total NPV and maximum total NPV is very small, which can be ignored. However, some paths may be more reasonable than those with the maximum total NPV, such as more stable ore yield and lower peel peak. Therefore, in the above algorithm, multiple optimal paths can be reserved and output for users to choose. When designing optimization software, the number of optimal paths to be reserved should be set by users as input data on the interface.
link