A heuristic method for production scheduling of an open pit mining operation

ABSTRACT

A mathematical model for production scheduling of open pit mines maximises the net present value and satisfies the reserves, pit slope angle, and production capacity constraints. A solution to this model aims to deliver an extraction sequence and the movement of materials across various stages within the operation over a defined planning horizon. However, given that the mineral reserves delineated into thousands of mining blocks form the geological input, the model falls in the category of large-scale optimisation problems, i.e. it is computationally intractable, and exact methods cannot solve realistic scenarios of the problem. Therefore, this paper contributes an alternative to the conventional mathematical model along with a corresponding heuristic method to solve this proposed model. An implementation of the proposed method at various realistic instances reveals better performance in terms of net present value, computation time and optimality gap as compared to the traditional methods.

KEYWORDS:

• Scheduling
• optimisation
• open pit mining
• heuristic methods
• exact methods

1. Introduction

The extraction of mineral resources in an open pit mine provides valuable (ore) and waste materials for the subsequent stages of the supply chain [1–3]. Figure 1 presents the framework of a typical open pit mining operation. A mixed integer programming (MIP) based mathematical model defines the movement of materials from sources to destinations within this framework at both long- and short-term planning stages [4]. The MIP model for a strategic or long-term planning horizon considers the economic (commodity price, operational costs, etc.), technical (mine, process, and refinery capacities, process recoveries, etc.), and geological data. A three-dimensional mineral resource (orebody) model forms the geological input. Given the size or extent of mineral resources, a typical orebody model constitutes thousands of uniform-sized blocks and the location (x, y, z coordinates), rock type, quality (metal content or grade), and quantity (tonnes) of material in each block become the basic geological inputs to the MIP model [3,5]. Knowing this information, a pit production scheduling formulation generates a block-by-block sequence of extraction and their destination within the system over a planning horizon.

Figure 1. The supply chain of a typical open-pit mining operation.

The location of a block within the three-dimensional orebody model establishes the extraction precedence, i.e. the extraction of a candidate block at lower levels may require a prior extraction of the blocks overlying this candidate block. This spatial relationship defines an extraction precedence of the blocks overlying a block located at lower level, and consequently, the mining industry practices a one-five pattern, i.e. for a block at level below, five overlying blocks at the next level above maintain an extraction precedence [6]. This precedence pattern guides the shape of the ultimate pit limit as the number of overlying blocks within the precedence matrix for a candidate block may increase or decrease with the depth of this block within the three-dimensional space and/or the rock type that controls the safe slope angle for geotechnical stability of the pit.

Given this background, an MIP model aims to generate an optimal schedule of block-by-block extraction that maximises net present value (NPV) and satisfies the mining reserves, block precedence and production (mine, process and refinery) capacity constraints over a planning horizon of an open pit mining operation. A list of sets, parameters and decision variables that describe the conventional MIP formulation [3] is presented below.

1.1. Sets and indices

$\mathit{t}\in \mathit{T}=$ set of periods or years,

$\mathit{b}\in \mathit{B}=$ set of blocks within the orebody model,

${\text{b}}^{\prime }$ set of blocks overlying (or have a mining precedence over) block $\mathit{b}$,

$\mathit{k}\in \mathit{K}=$ set of mineral (for example, copper, gold, silver, etc.) attributes,

$\mathit{d}\in \mathit{D}=$ set of material (ore and waste) destinations (for waste dump $\mathit{d}=1$; for ore processing streams $\mathit{d}>1$).

1.2. Parameters

$Q_b=$ quantity of material in block $\mathit{b}$ (tonnes),

$G_{\mathit{bk}}=$ grade or metal content of mineral attribute $\mathit{k}$ in block $\mathit{b}$ (% metal or troy ounce per tonne of ore or grams per tonne of ore),

$M^t=$ mining capacity during period $\mathit{t}$ (tonnes of material),

$P_d^t=$ processing capacity at ore destination $\mathit{d}$ during period $\mathit{t}$ (tonnes of ore),

$R_k^t=$ refining capacity of mineral attribute $\mathit{k}$ during period $\mathit{t}$ (tonnes or troy ounces or grams of metal),

$Y_{\mathit{kd}}=$ metallurgical recovery of mineral attribute $\mathit{k}$ at destination $\mathit{d}$ (%),

$S_k=$ commodity price of mineral attribute $\mathit{k}$ ($ per tonne of metal or $ per troy ounce of metal or $ per gram of metal),

$\bar{M}=$ mining cost ($ per tonne of material),

${\bar{P}}_{\mathit{d}}=$ processing cost at processing stream $\mathit{d},\mathrm{\forall d}>1$ ($ per tonne of ore),

${\bar{R}}_{\mathit{k}}=$ refining cost of mineral attribute $\mathit{k}$ ($ per tonne of metal or $ per troy ounce of metal or $ per gram of metal),

$v_{\mathit{bd}}=\left({\sum }_k\left(S_k-{\bar{R}}_{\mathit{k}}\right)\left(Q_b{G}_{\mathit{bk}}{Y}_{\mathit{kd}}\right)\right)-{\bar{P}}_{\mathit{d}}{Q}_b-\bar{M}{Q}_b$ = value ($) of block $\mathit{b}$ mined and moved to a processing destination $\mathit{d}>1$, such that, $v_{\mathit{bd}}=\left(\begin{array}{c}{v}_{\mathit{bd}};\\ -\bar{M}{Q}_b,\mathrm{if}{v}_{\mathit{bd}}\le 0\mathrm{\forall d};\end{array}\right.$

$\mathit{i}=$ annual discount rate (%).

1.3. Variables

$x_{bd}^t=\left\{\begin{array}{c}1,\mathrm{if}mining\mathit{block}\mathit{b}\mathrm{is}\mathrm{mined}\mathrm{and}\mathrm{moved}\mathrm{to}\mathrm{destination}\mathit{d}\mathrm{during}\mathrm{period}\mathrm{or}\mathrm{year}\mathit{t};\\ 0,\mathrm{otherwise}.\end{array}\right.$

Given set, parameter and decision variable definitions above, formulations (1)-(6) present the conventional model.

(1) $max\mathit{z}=\sum _{\mathit{b}\in \mathit{B}}\sum _{\mathit{d}\in \mathit{D}}\sum _{\mathit{t}\in \mathit{T}}\left[\frac{{\mathit{v}}_{\mathit{bd}}}{{\left(1+\mathit{i}\right)}^t}\right]x_{bd}^t$(1)

(2) $x_{bd}^t-\sum _{\mathit{\tau }=1}^t\sum _{\mathit{d}\in \mathit{D}}{x}_{b^{\prime }\mathit{d}}^{\mathit{t}}\le 0,\mathrm{\forall b},\mathit{d},\mathit{t}$(2)

(3) $\sum _{\mathit{d}\in \mathit{D}}\sum _{\mathit{t}\in \mathit{T}}{x}_{bd}^t\le 1,\mathrm{\forall b}$(3)

(4) $\sum _{\mathit{b}\in \mathit{B}}\sum _{\mathit{d}\in \mathit{D}}{Q}_b{x}_{bd}^t\le M^t,\mathrm{\forall t}$(4)

(5) $\sum _{\mathit{b}\in \mathit{B}}{Q}_b{x}_{bd}^t\le P_d^t,\mathrm{\forall t},\mathit{d}>1$(5)

(6) $\sum _{\mathit{b}\in \mathit{B}}\sum _{(\mathit{d}>1)\in \mathit{D}}\left(Q_b{G}_{\mathit{bk}}{Y}_{\mathit{kd}}\right)x_{bd}^t\le R_k^t,\mathrm{\forall k},\mathit{t}$(6)

Equation (1)(1) $max\mathit{z}=\sum _{\mathit{b}\in \mathit{B}}\sum _{\mathit{d}\in \mathit{D}}\sum _{\mathit{t}\in \mathit{T}}\left[\frac{{\mathit{v}}_{\mathit{bd}}}{{\left(1+\mathit{i}\right)}^t}\right]x_{bd}^t$(1) is the objective function and it maximises the discounted value of the cash flows realised during the planning horizon. Constraints (2) relate to the extraction precedence of mining blocks overlying other mining blocks located at levels (or benches) below. Similarly, Constraints (3) ensure that a block is mined once during the life of operation. Finally, Constraints (4) to (6) keep the mining, processing, and refining rates within the available capacities of these stages, respectively. Note that Constraints (5) consider blocks $\mathit{b}$ that correspond to $v_{\mathit{bd}}$ for a specific ore processing destination $\mathit{d}$. A simple form of the conventional open pit production scheduling model presented above (1)-(6) considers a single material attribute and a single ore processing stream (i.e. an ore processing stream and a waste dump) is given in Dagdelen and Johnson [7] and Ramazan [8]. This simple form of the conventional model is reproduced below.

(7) $max\mathit{z}=\sum _{\mathit{b}\in \mathit{B}}\sum _{\mathit{t}\in \mathit{T}}\left[\frac{v_b}{{\left(1+\mathit{i}\right)}^t}\right]x_b^t$(7)

(8) $x_b^t-\sum _{\mathit{\tau }=1}^t{x}_{b^{\prime }}^{\mathit{t}}\le 0,\mathrm{\forall b},\mathit{t}$(8)

(9) $\sum _{\mathit{t}\in \mathit{T}}{x}_b^t\le 1,\mathrm{\forall b}$(9)

(10) $\sum _{\mathit{b}\in \mathit{B}}{Q}_b{x}_b^t\le M^t,\mathrm{\forall t}$(10)

(11) $\sum _{\mathit{b}\in \mathit{B}}{Q}_b{x}_b^t\le P^t,\mathrm{\forall t},\mathrm{if}{v}_b>0$(11)

(12) $\sum _{\mathit{b}\in \mathit{B}}\left(Q_b{G}_b\mathit{Y}\right)x_b^t\le R^t,\mathrm{\forall t},\mathrm{if}{v}_b>0$(12)

Given the size of the realistic orebody models, an open pit mine production scheduling model translates into a computationally complex large-scale optimisation problem. Caccetta and Hill [5] confirm the computational complexity of the conventional model (7)-(12). Fathollahzadeh et al. [3] demonstrate the evolution in computational complexity as models transition from simple (7)-(12) to relatively realistic (1)-(6) versions.

Earlier studies have adopted numerous strategies to address the computational complexity of the pit production scheduling problem. In this context, Dagdelen and Johnson [7] share an exact method that applies Lagrangian relaxation with subgradient procedure for updating the multipliers to potentially solve the relatively large instances of the conventional model (7)-(12). Extensions of the pioneering work in Dagdelen and Johnson [7] are available in Akaike and Dagdelen [9] and Kawahata [10]. Akaike and Dagdelen [9] offer improvement in updating the values of Lagrangian multipliers such that the solution to the relaxed problem satisfies the production capacity constraints in the original problem. Similarly, Kawahata [10] applies Lagrangian relaxation to the expanded model (1)-(6) considering variable cut-off grades along with stockpiling and waste dump constraints. The algorithm in Kawahata [10] creates Lagrangian relaxation sub-problems and modifies the bounds on production capacity constraints to achieve a feasible solution to the relaxed problem and this procedure makes the algorithm suitable to solve large instances of the problem. Cullenbine et al. [11] solve the modified formulation of the problem (7)-(12) in Caccetta and Hill [5] through an application of a sliding time window heuristic that solves three sub-models over initial, middle, and later periods of the planning horizon and derive near optimal solution for relatively smaller instances of the problem. Later, Lambert and Newman [12] apply pre-processing and heuristics for variables elimination and an initial feasible solution, respectively, and then utilise Lagrangian relaxation of the resource or production capacity constraints to generate near-optimal solutions to relatively smaller (unrealistic) instances of the problem. Similarly, realising the importance of geological uncertainty, Ramazan and Dimitrakopoulos [13] propose a two-stage stochastic version (i.e. a stochastic integer programming model) of the formulation (7)-(12) and then solve this modified model using a heuristic procedure that generates a schedule sequentially by joining a number of periods within the planning horizon into computationally manageable groups. A few authors [14–16] have used dynamic programming as a solution approach to solve the original and modified versions of the model (7)-(12). However, dynamic programming is inefficient for large scale scheduling problems [17].

Caccetta and Hill [5] propose a branch and cut based exact approach that solves a modified version of the model (7)-(12) that defines the variable $x_b^t$ as if block $\mathit{b}$ is mined ‘by’ period $\mathit{t}$. The model in this study includes grade blending, stockpiling, and operational constraints in addition to the reserve, precedence, and production capacity constraints. The branch and cut algorithm in Caccetta and Hill [5] employ ‘best-first’ and ‘depth-first’ as search strategies along with a linear programming heuristic to generate solutions to large-scale problems in relatively short time. Later on, Boland et al. [18] and Bley et al. [19] rely on this pioneering work in [5] and solve instances of problems that consider thousands of blocks in the orebody model and up to 10 periods of scheduling horizon. In addition, Bienstock and Zuckerberg [20] share an LP relaxation approach for solving some of the largest instances of the problem to optimality in relatively short computation time, and so far, this work is the best in literature in the context of quality and timing of solutions. Similarly, Chicoisne et al. [21] propose a decomposition method that relies on the structure of the precedence constrained knapsack problem to solve the LP relaxation of the formulation (7)-(12), and then employ topological sorting as well as local search algorithms to offer feasible integer solution. Later, Vossen et al. [22] also propose a hierarchical Benders decomposition based procedure to solve LP relaxation with application to medium-sized instances of the problem. Samavati et al. [23] strengthen the LP relaxed formulation initially through block precedence and production capacity inequalities, and then use LP solutions to define a new weighting mechanism that facilitates their heuristic procedure. Recently, Letelier et al. [24] propose new pre-processing and cutting plane procedures to reduce the size of LP relaxed problem and generate good quality solutions for instances as large as 5 million blocks and 50 time periods. Contrary to LP relaxation approaches, Liu and Kozan [25] present network flow and conjunctive graph theory based algorithms for the formulation (7)-(12) and demonstrate the quality of their solutions against some of the heuristics reported in literature.

Aggregation of mining blocks with similar properties (grade, rock type, etc.) is another strategy that results in reducing the model size and leads to resolution of the computational complexity. Ramazan [8] applies linear programming approach to reduce the number of blocks by aggregation into ‘fundamental trees’, where each positive-valued fundamental tree that satisfies the slope or precedence constraints and may not be divided into smaller trees forms an input to the conventional model (7)-(12). Boland et al. [18] rely on binary variables only to ensure the slope or precedence requirement of the aggregated blocks, as continuous variables in their model capture the quantity of material moving from an aggregate as well as a block within an aggregate. Boland et al. [18] offer an aggregate partition strategy that refines the aggregates and facilitates a feasible solution to the original (without aggregation) problems with relatively large-scale data sets. Similarly, Eivazy and Askari-Nasab [26] utilise clustering algorithms for block aggregation and then applied the conventional model for large size models to develop medium-term production plans. Later, Jelvez et al. [27] propose an aggregation and disaggregation heuristic for the formulation (7)-(12) and offer better solutions for the data sets in public domain [28]. However, Fathollahzadeh et al. [3] express concerns that the matter of feasibility of the solutions based on block aggregates upon moving back to the original (without aggregation) models requires a focus in future studies.

Denby and Schofiled [29] present an application of genetic algorithms as a method to solve the classical pit production scheduling problem. However, this study confirms that the typical crossover and mutation strategies that improve the solution may result in excessive increase in solution time for large scale problems. Shishvan and Sattarvand [30] present several variants of the ant colony based metaheuristic algorithm to solve problem (7)-(12), and then select the best among these variants using valid comparisons to the production schedules generated by commercial software tools. Similarly, a number of recent studies [31–34] have contributed various metaheuristic methods to solve the modified model (7)-(12) and its extensions in the stochastic domain.

Osanloo and Gholamnejad [35] and Fathollahzadeh et al. [3] provide a comprehensive overview of the literature on formulations and solution methods for open pit production scheduling problem. Unlike earlier studies that rely on diverse aggregation strategies that employ a modified structure of the problem or its solution method, use new and robust procedures to solve LP relaxed version of problem or generate metaheuristic-based solution algorithms, the original contribution of this paper is a new formulation and a corresponding heuristic method that generates an integer solution to the pit production scheduling problem. More specifically, contrary to Equation (1)(1) $max\mathit{z}=\sum _{\mathit{b}\in \mathit{B}}\sum _{\mathit{d}\in \mathit{D}}\sum _{\mathit{t}\in \mathit{T}}\left[\frac{{\mathit{v}}_{\mathit{bd}}}{{\left(1+\mathit{i}\right)}^t}\right]x_{bd}^t$(1) in the conventional model, the objective function in the new formulation integrates the cost of not mining, processing, and refining the material in a specific period, i.e. a penalty component has been introduced in the objective function, which ensures that the solution process schedules the best (high grade) possible material during a specific period. This structure of the proposed method clearly differentiates it from the previous studies. In addition, the proposed method promises an integer solution to the problem. While the applications in Bienstock and Zuckerberg [20] and Letelier et al. [24] solve very large instances of the problem as reported in the literature so far, these studies are restricted to the LP relaxed version of the problem. Therefore, the derivation of the required integer solution to the relatively large instances of the problem remains an exclusive feature of the new formulation and the heuristic procedure. Moreover, the application of the proposed method at the data available in the public domain [28,36,37] confirms that the new formulation and the heuristic algorithm provide efficient solution to the realistic and relatively large instances of the problem. Next sections present the new formulation along with a discussion on how its structure may facilitate better solutions, heuristic method, numerical results, and conclusions.

2. New formulation and heuristic method

The proposed formulation in (13)-(23) is based on the the following sets, parameters, and decision variables. This list is in addition to some of the definitions shared previously.

2.1. Variables

$m^t=$ quantity of material mined during year $\mathit{t}$,

$m^{t^{\prime }}=$ quantity of material not mined during year $\mathit{t}$,

$p_d^t=$ quantity of ore processed at processing stream $\mathit{d}$ during year $\mathit{t}$,

$p_d^{t^{\prime }}=$ quantity of ore not processed at processing stream $\mathit{d}$ during year $\mathit{t}$,

$r_k^t=$ quantity of mineral attribute $\mathit{k}$ refined during year $\mathit{t}$,

$r_k^{t^{\prime }}=$ quantity of mineral attribute $\mathit{k}$ not refined during year $\mathit{t}$,

$x_{bd}^t=\left\{\begin{array}{c}1,\mathrm{if}mining\mathit{block}\mathit{b}\mathrm{is}\mathrm{mined}\mathrm{and}\mathrm{moved}\mathrm{to}\mathrm{destination}\mathit{d}\mathrm{during}\mathrm{current}\mathrm{period}\mathrm{or}\mathrm{year}\mathit{t};\\ 0,\mathrm{otherwise}.\end{array}\right.$

(13) $max\mathit{z}=\frac{\left[\left({\sum }_k\left(S_k-{\bar{R}}_{\mathit{k}}\right)r_k^t\right)-{\sum }_d{\bar{P}}_{\mathit{d}}{p}_d^t-\bar{M}{m}^t\right]}{{\left(1+\mathit{i}\right)}^t}-\frac{\left[\left({\sum }_k\left(S_k-{\bar{R}}_{\mathit{k}}\right)r_k^{t^{\prime }}\right)-{\sum }_d{\bar{P}}_{\mathit{d}}{p}_d^{t^{\prime }}-\bar{M}{m}^{t^{\prime }}\right]}{{\left(1+\mathit{i}\right)}^{t^{\prime }}}$(13)

(14) $x_{bd}^t-\sum _{\mathit{d}\in \mathit{D}}{x}_{b^{\prime }\mathit{d}}^{\mathit{t}}\le 0,\mathrm{\forall b},\mathit{d}$(14)

(15) $m^t-\sum _{\mathit{b}\in \mathit{B}}\sum _{\mathit{d}\in \mathit{D}}{Q}_b{x}_{bd}^t=0$(15)

(16) $p_d^t-\sum _{\mathit{b}\in \mathit{B}}{Q}_b{x}_{bd}^t=0,\mathrm{\forall d}>1$(16)

(17) $r_k^t-\sum _{\mathit{b}\in \mathit{B}}\sum _{\mathit{d}>1\in \mathit{D}}\left(Q_b{G}_{\mathit{bk}}{Y}_{\mathit{kd}}\right)x_{bd}^t=0,\mathrm{\forall k}$(17)

(18) $m^{t^{\prime }}-\left[\left(\sum _{\mathit{b}\in \mathit{B}}{Q}_b\right)-\left(\sum _{\mathit{b}\in \mathit{B}}\sum _{\mathit{d}\in \mathit{D}}{Q}_b{x}_{bd}^t\right)\right]=0$(18)

(19) $p_d^{t^{\prime }}-\left[\left(\sum _{\mathit{b}\in \mathit{B}}{Q}_b\right)-\left(\sum _{\mathit{b}\in \mathit{B}}{Q}_b{x}_{bd}^t\right)\right]=0,\mathrm{\forall d}>1$(19)

(20) $r_k^{t^{\prime }}-\left[\left(\sum _{\mathit{b}\in \mathit{B}}\sum _{\mathit{d}>1\in \mathit{D}}\left(Q_b{G}_{\mathit{bk}}{Y}_{\mathit{kd}}\right)\right)-\left(\sum _{\mathit{b}\in \mathit{B}}\sum _{\mathit{d}>1\in \mathit{D}}\left(Q_b{G}_{\mathit{bk}}{Y}_{\mathit{kd}}\right)x_{bd}^t\right)\right]=0,\mathrm{\forall k}$(20)

(21) $m^t\le M^t$(21)

(22) $p^t\le P^t$(22)

(23) $r_k^t\le R_k^t,\mathrm{\forall k}$(23)

Equation (13)(13) $max\mathit{z}=\frac{\left[\left({\sum }_k\left(S_k-{\bar{R}}_{\mathit{k}}\right)r_k^t\right)-{\sum }_d{\bar{P}}_{\mathit{d}}{p}_d^t-\bar{M}{m}^t\right]}{{\left(1+\mathit{i}\right)}^t}-\frac{\left[\left({\sum }_k\left(S_k-{\bar{R}}_{\mathit{k}}\right)r_k^{t^{\prime }}\right)-{\sum }_d{\bar{P}}_{\mathit{d}}{p}_d^{t^{\prime }}-\bar{M}{m}^{t^{\prime }}\right]}{{\left(1+\mathit{i}\right)}^{t^{\prime }}}$(13) presents the objective function maximising the discounted value of cash flows generated during the current period $\mathit{t}$ by penalising the value derived from the extraction of remainder of the orebody during the next period $t^{\prime }$. Here, $t^{\prime }$ equates to $\mathit{t}+1$, i.e. the remainder of the resource is assumed to be extracted, processed, and refined during the next period. Constraints (14) maintain extraction precedence of mining blocks overlying other mining blocks located at levels (or benches) below. Equation (15)(15) $m^t-\sum _{\mathit{b}\in \mathit{B}}\sum _{\mathit{d}\in \mathit{D}}{Q}_b{x}_{bd}^t=0$(15) -(17) define the production rates related to mine, processing streams, and refinery stages during current year, respectively. Similarly, Equations (18)(18) $m^{t^{\prime }}-\left[\left(\sum _{\mathit{b}\in \mathit{B}}{Q}_b\right)-\left(\sum _{\mathit{b}\in \mathit{B}}\sum _{\mathit{d}\in \mathit{D}}{Q}_b{x}_{bd}^t\right)\right]=0$(18) -(20) derive the quantity of material, quantity of ore, and the quantity of metal that is not mined, processed, or refined during the current period $\mathit{t}$, i.e. this material remains within the orebody for subsequent handling during time $t^{\prime }$ (i.e. after current period $\mathit{t}$). Finally, Constraints (21)-(23) ensure that the mining, processing, and refining rates remain within mining, processing, and refining capacities during current period $\mathit{t}$, respectively. Note that the proposed formulation (13)-(23) schedules the extraction of blocks and their destination for current period $\mathit{t}$, and accordingly, it is not generated for the whole planning horizon (i.e., $\mathrm{\forall }\mathit{t}$).

Figure 2 presents the steps of the heuristic procedure that implements this new formulation. The procedure creates and solves the new formulation successively for each period $\mathit{t}$, i.e. it generates a period-by-period production plan over a number of iterations of the algorithm. This sequential procedure confirms the heuristic nature of the algorithm. An update in the reserves leads to the removal of the mining blocks from the orebody block model, i.e. the mining blocks scheduled for extraction during current period $\mathit{t}$ will not be considered for subsequent periods. This structure allows reduction in number of binary variables over successive iterations of the algorithm. This confirms that the proposed strategy in this study is different from the heuristic procedure in Ramazan and Dimitrakopoulos [13]. Because, Ramazan and Dimitrakopoulos [13] solve the problem over a number of stages, where each stage develops schedules for a group of periods by allowing an overlap of a production period between two successive stages to account for available production capacities. In addition, the new formulation is structured to schedule the extraction of high-grade mining blocks as early as possible. This is realised through the second term in the objective function, because it incurs excessive penalty to the high-grade or high-value blocks allocated to the remaining reserves in time $t^{\prime }$ [38,39]. Note that Equations (15)(15) $m^t-\sum _{\mathit{b}\in \mathit{B}}\sum _{\mathit{d}\in \mathit{D}}{Q}_b{x}_{bd}^t=0$(15) -(20) facilitate the calculation of terms 1 and 2 in the objective function. This aspect differentiates the new formulation from the existing formulations, i.e. it confirms the originality of this research and delivers relatively high quality and efficient solutions to the open pit production scheduling problem (as reflected in the case studies presented in the next section).

Figure 2. Steps of the heuristic algorithm.

3. Numerical results

This section presents an application of the new ((15)-(23)) and the conventional ((7)-(12)) formulations. The application utilises kd, p4hd, and mclaughlin limit block models from the Minelib [28] and copper orebody block model in Asad and Dimitrakopoulos [36,37]. Tables 1 and 2 demonstrate the applicable economic, technical, and operational parameters.

Table 1. Economic and technical parameters for implementation of the new and the conventional methods.

Mineral	Case	Description	Value
Attribute	Number
Copper	1 to 14	Copper price ($ per tonne)	6000
		Refining, marketing, or sales cost ($ per tonne)	1000
		Mining cost ($ per tonne of material)	2.00
		Mill processing cost ($ per tonne of ore)	9.00
		Mill recovery (%)	90
		Leach processing cost ($ per tonne of ore)	2.00
		Leach recovery (%)	30
Gold, Silver, and Copper	15	Gold price ($ per oz.)	1200
		Gold refining, marketing, or sales cost ($ per oz.)	300
		Silver price ($ per oz.)	16
		Silver refining, marketing, or sales cost ($ per oz.)	8
		Copper price ($ per tonne)	6000
		Copper refining, marketing, or sales cost ($ per tonne)	1000
		Mining cost ($ per tonne of material)	2.00
		Mill processing cost ($ per tonne of ore)	19.00
		Mill recovery - Gold (%)	85
		Mill recovery - Silver (%)	80
		Mill recovery - Copper (%)	90
Gold	16	Gold price ($ per oz.)	1200
		Refining, marketing, or sales cost ($ per oz.)	300
		Mining cost ($ per tonne of material)	2.00
		Mill processing cost ($ per tonne of ore)	19.00
		Mill recovery (%)	90

Table 2. Operational parameters for implementation of the new and the conventional methods.

Case	Mining Capacity	Processing Capacity	Refining Capacity
Number	(tonnes)	(tonnes)	(tonnes/oz.)
1	2,000,000	1,500,000	15,000
2	4,000,000	3,000,000	15,000
3	6,000,000	4,500,000	30,000
4	9,000,000	6,000,000	50,000
5	11,000,000	7,000,000	60,000
6	12,000,000	8,000,000	60,000
7	14,000,000	9,000,000	75,000
8	15,000,000	10,000,000	75,000
9	17,000,000	11,000,000	80,000
10	18,000,000	12,000,000	90,000
11	22,000,000	14,000,000	100,000
12	22,000,000	12,000,000	90,000
13	100,000,000	82,000,000	1,000,000
14	100,000,000	82,000,000 (Mill)	60,000
		10,000,000 (Leach)
15	50,000,000	13,000,000	500,000 (Gold, oz.)
			11,000,000 (Silver, oz.)
			20,000 (Copper, tonnes)
16	11,000,000	7,000,000	1,000,000 (Gold, oz.)

Given the input parameters in Tables 1 and 2, 16 diverse instances (cases) of the proposed (heuristic) and traditional (Exact) methods have been implemented using Microsoft Visual Studio environment with IBM CPLEX Concert Technology [40] in an Intel^(R) Core^(TM) i7-4712MQ CPU at 2.3 GHz and 8GB RAM system. The detail of these cases is as follows:

1. Cases 1 to 4 - application of the new and conventional formulations at subsets of the kd block model for small size models,

2. Cases 5 to 12 - application of the new and conventional formulations at subsets of the kd block model for relatively large size models. Given the problem size, the new formulation generates an integer solution and the conventional formulation is able to generate only the LP-Relaxed solution.

3. Case 13 - application of new formulation and heuristic method at copper orebody model in Asad and Dimitrakopoulos [36,37] for a large size (realistic) instance of the problem considering single processing option as ore destination.

4. Case 14 - application of new formulation and heuristic method at copper orebody model in Asad and Dimitrakopoulos [36,37] for a large size (realistic) instance of the problem considering multiple processing options as ore destinations.

5. Case 15 - application of new formulation and heuristic method at p4hd model for a large size (realistic) instance of the problem considering multiple material attributes (copper, gold, and silver).

6. Case 16 - application of new formulation and heuristic method at mclaughlin limit model for the largest size (realistic) instance of the problem considering.

As indicated above, Cases 1 to 4 are designed to establish the performance of the new formulation through valid comparisons with the conventional formulation. Note that Cases 1 to 13 correspond to the simplest version of the new and conventional formulations, i.e. the destinations in the system include a waste dump (d = 1) and a single processing stream (d = 2) for a single material attribute (K = 1; copper in kd and Asad and Dimitrakopoulos [36,37] block models). Similarly, Cases 14 to 16 facilitate the application of new formulation at diverse (multiple destinations, multiple materials, etc.), large scale, and realistic instances of the pit production planning problem. Table 3 summarises the numerical results of these implementations. It indicates that the heuristic procedure that implements new formulation generates better solutions for Cases 1 to 4 in relatively less time and optimality gap as compared to the exact method that solves the conventional formulation. In addition, the exact method does not offer integer solutiosn for the Cases 5 to 12, while the proposed heuristic method even generates an integer solution to the largest instance (Case 12) within this category in merely 13.93 minutes at 6.92% gap. Equation (24)(24) $\mathit{Gap}=\left[\frac{\left(\mathit{BB}-\mathit{BI}\right)}{\mathit{BI}}\right]100$(24) defines the percent gap ($\mathit{Gap}$) for the heuristic method. Here, $\mathit{BB}$ represents the best bound and $\mathit{BI}$ represents the best integer solution. This gap analysis in Cases 5 to 12 relies on the LP-Relaxed solution from the conventional formulation. Accordingly, LP-Relaxed solution value (NPV) is used as $\mathit{BB}$ for calculation of the $\mathit{Gap}$ value.

Table 3. Performance evaluation of the new heuristic and the conventional exact methods.

Case Number	Number of Blocks	Number of Periods	Number of Binary Variables	Number of Constraints	Exact Method (CPLEX)			Heuristic Method
					Solution Time	NPV	Gap	Solution Time	NPV	Gap
					(Seconds)	($M)	(%)	(Seconds)	($M)	(%)
1	687	6	4,122	9,597	5,085	78.87	0.05	52	78.89	0.035
2	1,358	6	8,148	23,888	23,480	160.72	0.86	52	161.37	0.45
3	2,289	6	13,734	46,203	144,000	292.95	1.35	120	293.40	1.19
4	3,450	6	20,700	75,786	180,000	455.69	4.30	221	465.08	2.19
-					LP-Relaxed Infeasible Solutions			Integer Solutions
5	4,641	7	32,487	126,063	1,660	667.68	LP-Relaxed	213	634.60	5.12
6	5,638	8	45,109	181,286	212	828.75		240	791.07	4.76
7	6,550	8	52,400	216,158	412	1,013.54		385	962.91	5.25
8	7,601	8	60,808	255,425	455	1,103.75		316	1,055.48	4.57
9	8,376	9	75,384	318,732	814	1,183.29		307	1,130.35	4.68
10	9,568	9	86,112	367,678	1,050	1,280.96		273	1,214.01	5.50
11	10,582	9	95,238	409,048	1,614	1,376.31		355	1,296.33	6.16
12	11,610	10	116,100	500,390	3,114	1,366.62		836	1,278.16	6.92
13	53,725	14	752,150	3,395,133	No Solution			3,044	10,906.85	-
14	53,725	14	752,150	3,395,133	No Solution			4,237	10,926.05	-
15	40,947	7	286,630	1,335,996	No Solution			1,353	396.08	-
16	110,916	11	1,220,076	5,649,812	No Solution			18,612	1,323.44	-

(24) $\mathit{Gap}=\left[\frac{\left(\mathit{BB}-\mathit{BI}\right)}{\mathit{BI}}\right]100$(24)

Once the performance of the proposed heuristic procedure in terms of the quality as well as the efficiency in timing of solutions is confirmed through Cases 1 to 13 (simplest version of the formulations), then its application in Cases 14 to 16 is diversified to account for relatively larger instances of the problem along with multiple destinations (Case 14) and multiple materials (Case 15) options. Case 16 constitutes about 1.22 million binary variables and 5.65 million constraints, and therefore, it forms the largest instance of problem. The new formulation and heuristic method solved this case in record 5.17 hours.

It is important to note that both the exact and the proposed heuristic algorithms (see Figure 2) employ CPLEX solver for deriving solutions to the 16 instances illustrated in Table 3. This strategy ensures valid comparisons of the conventional and new formulations. The final results for conventional approach (Cases 1 to 4) presented here have been computed after several experiments with the two solution control parameters including the minimum optimality gap and the maximum solution time (whichever is achieved first). This forms the reason for a higher than zero optimality gap for the exact method solving the conventional formulation. These results are then utilised as a reference or benchmark to validate the quality of solutions derived from the proposed heuristic algorithm. The purpose here is to emphasise the fact that while the conventional formulation takes prohibitively longer computation times for small to medium size instances, the proposed formulation and the heuristic procedure derive near optimal solutions for the same instances of the problem in significantly short computation time with an acceptable optimality gap. Additionally, even the LP-Relaxed version of the conventional formulation for Cases 13 to 16 did not converge to a solution within 60 hours in the machine (Intel^(R) Core^(TM) i7-4712MQ CPU at 2.3 GHz and 8GB RAM) used for all 16 cases.

4. Conclusions

Open pit mine production scheduling problem is classified as computationally complex (NP-hard). This paper demonstrates that an increasing trend in input data (specifically, the geological input) size leads to an exponential growth in the number of binary variable and constraints, which precludes the solution to this problem. Therefore, deriving an optimal solution using exact methods has been a challenge. While earlier studies like Bienstock and Zukerberg [20] and Letelier et al. [24] provide quality solutions for very large instances of the LP-Relaxed version of the problem, the new formulation and the corresponding heuristic method have the capability to provide integer solutions for relatively large-scale realistic instances of this NP-hard problem within an acceptable time frame. An application at several diverse cases confirms the robustness and computational efficiency of the proposed method. However, the heuristic method relies on CPLEX solver [40] and it is a disadvantage that restricts its commercial applications. In addition, the formulation and the heuristic method are limited in application to the deterministic version of the problem and future work may focus on the stochastic domain by considering the commodity price and/or geological uncertainties that risk the achievement of defined objectives.

Disclosure statement

No potential conflict of interest was reported by the authors.

Data availability statement

The datasets used or analysed in this study are available in the public domain at Minelib espinoza2013minelib. The link to these datasets in Minelib is http://mansci-web.uai.cl/minelib/Datasets.xhtml.