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⁵. Production increments for production scheduling can be obtained by the “parametric analysis” proposed by Lerchs and Grossmann⁶. 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 decision-making 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 multi-period decision-making 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¹, 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.

Discrete body dynamic sorting.

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₁ in the first year and J₃ in the second year, linked by dotted arrows, indicate that the open pit advances to J₁ in the first year and J₃ in the second year. In geometric space, J₃ completely contains J₁, therefore the quantity of ore corresponding to J₁ and J₃ can be expressed as O₁,₁ and O₂,₃, 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₁ and O₃; R_1,1 and R_2,3 are the same as R₁ and R₃, 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:

where γ is the ore recovery rate and δ is the mixing rate of waste rock.

The concentrate amount v_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:

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_i-1,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:

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_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:

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_i-1,k(i, m). According to Eq. (1), the quantity of ore produced by state transition at two adjacent time points is o_i-1,_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:

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:

where a refers to the infrastructure investment base unrelated to the production scale, 10⁴ ¥, 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:

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ⁿ, q_Uⁿ] 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ⁿ and upper bounds q_Uⁿ of the constraint domain [q_Lⁿ, q_Uⁿ]:

  $$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¹, 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₁, x₂, x₃ 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ⁿ > q_up, LastPlan = true.

• Step three: Find all production scheduling whose annual ore output meets the constraints of the domain [q_Lⁿ, q_Uⁿ]. That is, find all subsequences in which the annual ore output falls into the domain [q_Lⁿ, q_Uⁿ]. 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ⁿ, q_Uⁿ]. 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ⁿ, the lower limit of the constraint domain [q_Lⁿ, q_Uⁿ] 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₁, concentrate v₁ and stripped waste rock r₁ in this year. Then, calculate the annual profit u₁ and its net present value NPV₁ 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(i-1) 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:

    1. (A)

       If o_i < q_Lⁿ and o_i + r_i ≤ q_up.

    2. (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.

    3. (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.

    4. (D)

       q_Lⁿ ≤ o_i≤ q_Uⁿ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_idiscounted 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.

    5. (E)

       If o_i > q_Uⁿ 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.

  • 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ⁿ. 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:

    1. (A)

       If o_i< q_Uⁿ 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_idiscounted 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.

    2. (B)

       If o_i > q_Uⁿ 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.

  • 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(i-1) 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:

    1. (A)

       if o_i< q_Lⁿ and o_i + r_i ≤ q_up, in two cases:

    2. (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.

    3. (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ⁿ ≤ o_i≤ q_Uⁿ 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.

    4. (D)

       If o_i > q_Uⁿ 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.

  • 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 best-planned path on the interval [q_Lⁿ, q_Uⁿ].

• 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ⁿ 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ⁿ, q_Uⁿ) is, the more schedules satisfy the constraint of (q_Lⁿ, q_Uⁿ). If the domain is set wide, the number of schedules to be evaluated in step three is too large and time-consuming. 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.

Flow chart of enumeration optimization algorithm for production scheduling.

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.