A heuristic approach for the stochastic optimization of mine production schedules

Traditionally, mining engineers plan an open pit mine considering preestablished conditions of operation of the plant(s) derived from a previous plant optimization. By contrast,mineral processing engineers optimize the processing plants by considering a regular feed from the mine, with respect to quantity and quality of the materials. The methods implemented to optimize mine and metallurgical plans simultaneously are known in the mining industry as global or simultaneous optimizers. The development of these methods has been of major concern for the mining industry over the last decade. Some algorithms are available in commercial mining software packages however, these algorithms ignore the inherent geological uncertainty associated with the deposit being considered, which leads to shortfalls in production, quality, and expected cashflows. This paper presents a heuristic method to generate life-of-mine production schedules that consider operating alternatives for processing plants and incorporate geological uncertainty. The method uses iterative improvement by swapping periods and destinations of the mining blocks to generate the final solution. The implementation of the method at a copper deposit shows its ability to control mine and processing capacities while increasing the expected net present value by 30% when compared with a solution generated using a standard industry practice.


Introduction
Mining complexes contain multiple sequential activities that are strongly interrelated: (1) mining the materials from one or multiple sources; (2) blending the material considering stockpiling, (3) transforming the material in different processes or processing paths; (4) transporting the products to final destinations, etc. The quality of the input material of a metallurgical process may determine its corresponding throughputs, costs and metallurgical recoveries. Mill throughput can be sensitive to rock hardness, work index or the ratio of clay materials; costs and reaction times in an autoclave depend on sulphur content; recoveries are affected by deleterious materials (Wharton 2005). Multiple approaches have been developed to optimize the different parts of a mining complex in isolation: (Caccetta and Hill 2003;Lerchs and Grossman 1965;Picard 1976) for pit design and mine production scheduling. The process of optimizing all activities of a mining complex simultaneously is known in the mining literature as global optimization (Whittle 2007(Whittle , 2010. This is a problem with high complexity due to the link between time periods and discounting, the blending requirements, the flexibility generated from the stockpiles, the multiple processing alternatives, and the variability and uncertainty associated with grades and physical characteristics (Whittle 2007).
Over the last decade, several algorithms that seek for generating optimal solutions in mining and processing plans have been developed. Hoerger et al. (1999) formulate the problem of optimizing the simultaneous mining of multiple sources (pits and underground mines) and the delivery of ores to multiple plants as a mixed integer program. The model calculates the net present value of the mining complex by using variables that represent material sent from the mines to the stockpiles, material sent from the mines to the processes, and material sent from the stockpiles to the processes and their corresponding associated costs. The blocks are grouped into 'increments' based on the metallurgical properties, which belong to sequences (or mining phases). The integer variables are used to model mine sequencing constraints at a phase level and plant startups and shutdowns. Due to the use of phase sequencing constraints instead of block sequencing constraints to decrease the complexity of the problem, there is a loss of resolution in the solution generated from the method that may lead to the inability of meeting the blending and production requirements. Furthermore, the method does not consider multiple operating alternatives for each process and ignores the geological uncertainty associated with the ore deposits. Whittle (2007) presents the Prober optimizer for global optimization that aggregates the mining blocks into parcels of mine material type. These material types are defined from different grade bin categories; that is, for each relevant grade or attribute, cut-offs are defined to allow flexibility for blending purposes. The method considers stockpiles for each material type that may be combined with the material obtained directly from the mines to satisfy the different process requirements. Prober uses a random sampling and a local optimization approach to generate the solution. The random sampling consists of a search algorithm that samples the feasible domain of alternative LOM (life-ofmine) mining plans; the local optimization is an evaluation routine that determines the optimal cut-off grade, stockpiling, processing selection, blending and production plan and determines the NPV (net present value). The optimizer works by repeatedly creating a random feasible solution and then finding the nearest local maximum. The various NPVs that the algorithm finds are stored, and the run is usually stopped when the top ten values lie within 0.1% of each other. Although very flexible and able to handle complex blending operations, the algorithm has some drawbacks: it groups the parcels into panels and assumes that the parcels are consumed in the same proportion within a panel; good solutions may be found but it does not guarantee optimality; and, geological uncertainty is discarded. Several methods are available to model geological uncertainty by means of stochastic simulations (Chatterjee et al. 2012;Goovaerts 1997;Zhang et al. 2012).
Regarding geological uncertainty, some approaches have been developed in the last decade to account for grade and material type uncertainties into pit design and mine production scheduling. Ramazan and Dimitrakopoulos (2013) formulate the mine scheduling problem as a two-stage stochastic integer program (SIP) in which the first stage variables represent mining decision variables and the second stage variables represent deviation from grade and production targets evaluated on a set of orebody simulations. The formulation can be extended to include stochastically designed stockpiles, multiple processors and integrate short-term information (Bendorf and Dimitrakopoulos 2013;Leite and Dimitrakopoulos 2014;Dimitrakopoulos and Jewbali 2013). Menabde et al. (2007) develop and implement a method that accounts for geological uncertainty and simultaneously optimizes the sequence of extraction of the mining blocks and the cut-off grade policy. The authors aggregate blocks into panels to reduce the number of binary variables and obtain an increase of 26% in expected NPV when compared to a solution that uses a deterministic marginal cutoff grade policy. Boland et al. (2013) propose a multistage stochastic programming approach that considers the decision of processing block aggregates as posterior-stage variables. The approach provides a set of policies to follow according to the actual scenario (orebody) obtained with the advance of the extraction. The implementation of the approach using realistic mining data increases the expected NPV by 3% when compared to a conventional deterministic method. However, some drawbacks can be noted in the approach: continuous variables with aggregates do not guarantee slope constraints; it assumes all scenarios can be covered with orebody simulations; and, it does not penalize deviation from production targets.
Although, the SIP formulation generates substantial improvements in terms of NPV and meeting production targets, industry standard optimizers such as CPLEX are unable to solve big size problems due to the large amount of integer variables, thus alternative solution avenues are being sought (Lamghari and Dimitrakopoulos 2012). Many different approaches are available to solve large combinatorial optimization problems. Some of them have been implemented for solving complex mine scheduling optimization problems. Godoy (2003) and Godoy and Dimitrakopoulos (2004) develop a multi-stage method for mine production scheduling that integrates the joint local geological uncertainty and uses the simulated annealing (SA) algorithm. The method seeks to generate a risk-based mine production schedule that minimizes deviation from ore and waste production targets. Leite and Dimitrakopoulos (2007) apply the method at a copper deposit obtaining an expected NPV 20% greater than the ones obtained using conventional deterministic schedulers. Albor and Dimitrakopoulos (2009) implement the method at the same copper deposit and observe that the schedule obtained was not sensitive after 10 orebody simulations. A similar study on the numer of simulations is performed with the method proposed in this paper. Albor and Dimitrakopoulos (2009) point out that the stochastic final pit limit was 17% greater than the deterministic one, adding 9% to the expected NPV. Goodfellow and Dimitrakopoulos (2013) develop a simulated annealing implementation for pushback design to control deviation from pushback size targets considering different material types and processing plants. Lamghari and Dimitrakopoulos (2012) implement tabu serach (TS) and variable neighbourhood search (VNS) for the mine scheduling problem obtaining near-optimal solutions while outperforming CPLEX in terms of computational time. Lamghari et al. (2013) develop a hybrid approach that combines exact methods and metaheuristics for solving the LOM production scheduling problem.
In this paper, a heuristic approach that iteratively perturbs an initial solution to generate a mine plan that accounts for geological uncertainty is described and implemented in a copper deposit. The main contribution is the incorporation of geological uncertainty and operating alternatives in the processing plants to the mine production scheduling problem. Section 3 describes the proposed model and the solution approach, whereas Sect. 4 outlines its implementation at a copper deposit. Section 5 concludes the paper.

Preliminaries
In mining operations, the different types of material extracted from a deposit are sent to stockpiles or the different processing streams (destinations). Processing a particular mining block can be profitable; however, the decision of mining and processing that particular block in a given period not only depends on the individual characteristics of the block (grades, metallurgical properties, etc.) but in the compound characteristics of the total material sent to the destination, including both: the material sent directly from the deposit and the material sent from the stockpiles. At any destination, there may be multiple operating conditions depending on the quality of the input material and the desired properties of the output products; that is, a given destination (mill) may have several operating alternatives (operate the mill with high silica in the input material or with low silica content). Some additives may be considered at each destination, e. g., cyanide in a hydrometallurgical plant. Each operating alternative at a given destination may have its corresponding costs, metallurgical recoveries, blending requirements, and additive demand (Fig. 1). Considering the example of the mill operated with higher silica content in the input material, given the hardness of the silica, the demand of energy (seen as an additive of the mill) to decrease the average particle size to 80 µm is higher when compared to a low silica content operation. This may affect the costs and recoveries there-after; that is, having different blending requirements in the input material for a particular process originates different consumption of additives, costs and metallurgical recoveries. However, these differences may also be originated from different conditions of operation in the process, or different product specifications; e.g., 120 µm instead of 80 µm as average particle size desired in the output material of the mill.The production plan of a mining complex must state for any mining block when to mine it, but also where to send it; and, for each destination, which operating alternative to use.

Optimization model
Two main decision variables are used to evaluate tonnages and cash flows in the stochastic formulation of the mine production scheduling problem: X itd is a binary variable representing whether or not a block i is mined in period t and sent to destination d, and Y tdo is a continuous variable representing the proportion of material sent to destination d in period t that is processed using alternative o. The amount of material mined at a given period t in simulation s can be evaluated as follows: where m is is the mass of block i in simulation s, N is the number of mining blocks, and D is the set of available processing destinations. For modelling purposes, destination d = 0 represents sending the block to the waste dump, and d = D + 1 represents sending the block to the stockpiles. The amount of material that will be processed at a given destination accounts for the material that comes directly from the mine and the material that comes from the stockpiles: where tonne_rehandle(s, t, d) is the material sent from the stockpiles to destination d in period t and simulation s. This amount of material rehandle depends on the tonnage available in the stockpiles and the idle capacity associated with a given destination, and follows the mass conservation expression: At a given processing destination, all the material must be processed using any of the available operating alternatives: The revenue in each period and destination is obtained by accounting for the different metals recovered from the processing destinations. For each operating alternative at any given destination, there is a recovery expression associated to each material type and metal. Similarly, processing costs and requirement of additives depend on the type of material input to a destination and the operating alternative implemented. The costs associated to the different activities of the mining complex are evaluated as follow: where mine_cost (s, t) is the total cost of mining in period t and simulation s, m c is the per-unit mining cost, process_cost (s, t) is the total processing cost in period t and simulation s, p c (d, o) is the per-unit processing cost, stockpile_cost (s, t) is the total stockpiling cost in period t and simulation s, sp c , rehanlde_cost (s, t) is the total rehandle cost in period t and simulation s, and rh c is the per-unit rehandle cost. The objective function of the proposed formulation is given by the sum of the discounted cash flows in the different periods and simulations of the mineral deposit considered: where S is the set of orebody simulations, T is the number of years considered for the project and d is the discount rate. Given the time value of money and the geological uncertainty associated with the deposit, the blocks with higher and more certain profit must be mined in early periods and sent to their optimal destinations, whereas the blocks which are more certain to be non-profitable must be delayed for latest periods and sent to the waste dump. The amount of metal v sent to a given destination can be evaluated as: where g isv is the grade of metal v in block i and simulation s. The average grade of metal v in processing destination d in period t and simulation s is given by: Grade requirements can be impossed as blending constraints as follows: where low and high_range are the operational ranges of property v in destination d using operating alternative o. The requirement of additive a at each destination can be evaluated as: where k(d, o, a) is the per-unit demand of additive a in destination d using operating o. Availability contraints ensure that the consumption of additive in a given destination and period does not exceed the available amount: The feasible domain is also constrained by mine and processity capacity constraints: A mining block can only be mined once and sent to a unique destination: Precedence constraints ensure that a block can only be extracted if the set of predecessor blocks has been extracted. The set of predecessor blocks is given by the slope angles required based on a geotechnical study: All variables must be greater or equal to zero. Y tdo ≤ 1 and X itd ∈ {0, 1}. Given the complexity of the problem derived from the flexibility considered at the different stages of the mining complex, the use of an exact method incorporated in any conventional optimization software, such as CPLEX, will not be able to generate an optimal solution in a feasible amount of time. A heuristic methodology is proposed to solve the problem. The methodology is presented in the next section. Capacity, availability and blending constraints are called target constraints in the method for simplicity.

Solution
The method proposed in this paper uses iterative improvement over an initial solution until it converges to a final one. The procedure that the method uses can be divided in three stages: (i) assign periods and destinations to the mining blocks based on the initial solution; (ii) calculate the overall profitability per block per destination based on the orebody simulations; and (iii) perturb the solution until a stopping criteria is reached to generate the final solution (Fig. 2).

Stage 1
In this stage, the method assigns periods and destinations to the mining blocks from the initial solution. It also assigns a geotechnical zone for each block based on the geotechnical zone file. Different zones can have different slope angles and therefore different set of predecessor and successor blocks. If there are some slope constraint violations in the initial solution, block mining period corrections are performed based on the different slope angles. For doing so, when a slope constraint violation is found, the mining period of the block is moved to a feasible period based on the set of successor and predecessor blocks; that is, the range between the latest period of the predecessor blocks and the earliest period of the successor blocks (Fig. 3).

Stage 2
At this stage, the profits and costs for each simulation and period are evaluated. From the material types and grades at each simulation, the proposed method calculates, for Overall profitability per block per destination each block, the overall profitability per available destination; that is, it evaluates the profit (or loss) obtained by sending a particular block to a given destination (considering the best operating alternative for that block at that destination) and accumulates it through the set of simulations (Fig. 4). From there, the method evaluates the optimal destination for a particular block, but, this optimal destination may not be the final destination due to capacity, availability and blending constraints (target constraints).

Stage 3
This is the perturbation stage. A block is selected randomly and the available destinations for that particular block are sorted based on its overall profitability. If the best destination has a positive overall profitability, i.e., it increases the value of the objective function, the block is pushed to early periods, otherwise it may be pushed to later periods.
For positive overall profitable blocks, the method defines four possible options for periods and destinations (Fig. 5). The first option is to send the block to its best destination in the previous period (current period-1). If there are no slope and target constraint violations this option is chosen. The second option is to mine the block in the previous period and send it to a profitable destination different from the optimal without violating slope and target constraints; that is, it considers the destinations with positive overall profit. The third option is to randomly select another block mined in the previous period from which a double swapping that increase the objective function can be performed without violating slope and target constraints. The double swapping consists of two different blocks switching mining periods. If the double swapping is non-feasible or non-profitable, the block is sent to the stockpile, which is the last option.
If the block has a negative overall profit for all the different destinations, it is sent to the waste dump. To decide the period when the block is going to be mined, the method evaluates the overall profitability of the set of closest successor blocks. If the sum of the overall profitability of the closest successor blocks is positive, the period of the block does not change to allow the successor blocks to move to early periods. This permits the schedule to access profitable areas early even when waste blocks are overlying them. If the sum is negative, the block and the predecessors belonging to the same period are sent to the next period without violating slope and mine capacity constraints. It should be noted that the method uses an overall revenue cut-off instead of a grade cut-off that conventional methods use to discriminate between ore and waste. The material sent to stockpiles is a profitable material that cannot be processed immediately due to capacity, availability and blending constraints. After a certain number of iterations, the method re-evaluates the destination policy of the stockpiled material based on the tonnage available, its profitability and the respect of the target constraints at the different destinations.The method stops after the maximum number of iterations, swaps, or iterations without substantial improvement of the objective value, are reached. The maximization of the objective value is driven by means of the swapping mechanism: sending the most profitable blocks to early possible periods and the best available destinations and sending to the waste dump the blocks with negative overall profit in later possible periods without violating slope and target constraints. The constraints are respected by means of the checking mechanism throughout the iterations of Stage 3. Figure 6 shows the flowchart of the method in the Stage 3.

Case study: a copper deposit
This section describes the implementation of the proposed method at an actual copper open pit mining operation containing one deposit. The mineral deposit within the pit contains 176,138 mining blocks scheduled to be extracted over a time horizon of 18  years, which results in an optimization problem with more than 3 million integer variables. Additional examples of the application of the method can be found in Montiel (2014). Fifty simulations of the copper orebody are available to quantify the geological uncertainty of the mineral deposit. Note that 25 of these simulations are used in the application of the proposed approach and the remaining to perform a risk analysis of solutions generated. The sensitivity of the proposed method to the number of geological scenarios is also assessed. Figure 7 shows the different material types of the deposit and the available processing destinations in the mining complex. Figure 8 shows a solution generated using a conventional optimizer that could be used as an initial solution for the proposed method. The method seeks to generate a mine and destination schedule that maximizes the NPV and respect capacity, availability and blending constraints. For doing so, an initial solution is iteratively perturbed to improve the objective value. To identify the number of perturbations required in the perturbation stage, different numbers of perturbations are tested and the deviations from target constrains evaluated. Figure 9 shows the evolution of the 50th percentile (P50) of the tonnage sent to process 1 and the total tonnage mined with the number of perturbations. It is observed a large deviation from capacities at small number of perturbations and a substantial reduction in the deviation driven by the increment of the number of perturbations; i.e., the reduction in deviation from capacity of process 1 decreases from 9% on average to 0.2% when increasing the number of perturbations from 100 thousand to 1 million. Regarding the total mine production, the average deviation in the first 16 years remains in the same level (around 4% from the mine capacity). An analysis based on the value of the expected NPV is performed. Figure 10 shows the sensitivity of the expected NPV with the number of perturbations and the number of simulations. It is observed that the increment of the expected NPV is marginal after 1 million perturbations. There is no substantial benefit in increasing the number of perturbations there-after. The objective value is increased by 30% when compared to the initial solution. The same analysis is done for the number of simulations required. It can be observed that after 10-15 simulations, no particular improvement in the expected NPV is presented. This is consistent with the results observed when implementing other stochastic approaches for mine production scheduling (Albor and Dimitrakopoulos 2009;Dimitrakopoulos and Jewbali 2013) and it is due to the so-called support scale effect: a schedule groups thousand of mining blocks in just one period, which means that the sensitivity of a schedule is not the same observed at the scale of individual mining blocks and explains the insensitivity of the stochastic schedule to additional scenarios after a point. Figure 11 shows the tonnage sent to Process 1, 2, 3 and the total tonnage mined. Given that the solution states the destination of the blocks, the differences in tonnage of the material sent to a given destination through the different simulations are negligible. These minor differences are generated from different tonnages of blocks among simulations derived from simulated densities. If the tonnage of the blocks were similar along the different simulations, no differences were presented in terms of tonnage among simulations. It can be observed that the Process 1 and the total tonnage mined are controlled by their corresponding capacities, whereas the Processes 2 and 3 are controlled by the amount of profitable reserves for those destinations. Although the material sent to the different destinations does not vary significantly between simulations, the amount of metal that can be recovered has significant fluctuations (Fig. 12). This is originated from the grade and material type uncertainties; that is, the amount of metal sent to a process change in the simulations due to copper grade uncertainty, and the metallurgical recovery at a given destination vary in the simulations due to material type uncertainty. Figure 13 shows the amount of the different material types sent to the stockpiles. Sending waste material to the stockpile may be seen as a misclassification error. Although there are some risks of misclassifying material by following the stochastic solution generated, the algorithm seeks for minimizing the misclassification errors; that is, sending waste material to the stockpiles or sending a given material type to a non-profitable destination. A given material type sent to a wrong destination may produce a very high cost with low or negligible recovery. The way the algorithm controls misclassification errors is by maximizing the objective function, given that misclassification errors are very costly. The amounts of material type 3, 4 and 5 sent to stockpile are marginal and may be generated from misclassification; however, the Fig. 13 Material types sent to stockpile algorithm will send that material to a particular destination if there is some profit associated. By having a look at the output material from the stockpiles, it is observed that only ore types 1 and 2 are rehandled and sent to process 1. Figure 14 shows the P10, P50 and P90 of the cumulative discounted cash flows. It is observed that during the first two decades of the project, the expected NPV is around $7.8 billion. Although, no blending constraints were considered in this case study, the method attempts to maximize net present value expectations while maintaining target constraints within acceptable tolerable limits. It discriminates blocks between ore and waste based on the overall profitability; that is, the profit (or loss) obtained by sending a block to a given destination accumulated through the set of simulations. When profitable material cannot be processed due to target constraints, it is sent to the stockpiles for being rehandled in future periods.
The results show the ability of the proposed method to improve the possible deviation of the production targets by incorporating and managing uncertainty in the perturbation mechanism. To assess the results, this solution is compared with a solution generated using a conventional approach that is a standard mining industry practice. This practice aims to maximize net present value and includes the use of mine planning  and design software packages that are typically based on a parametric implementation of the Lerchs-Grossman algorithm (L-G) (Lerchs and Grossman 1965), combined with mixed integer linear programming. A solution generated following the standard conventional approach is compared with the risk-based solution from the proposed method. For comparison purposes, both solutions are evaluated using the same set of 25 simulations (geological scenarios); note that these simulated scenarios are different from the ones used to generate the risk-based solution. Figure 15 shows the tonnage sent to process 1 with the conventionally generated solution. Large impractical deviations from the capacity of Process 1 can be observed. Figure 16 shows the risk profile of the cumulative discounted cash flow of the conventional schedule. It is observed that during the first two decades of the project, the expected net present value of the risk-based schedule is 30% greater than the conventional initial schedule. This shows the ability of the method to handle two conflicting objectives: maximize expected net present value while approaching target constraints. Furthermore, the risk-based solution shows similar risk profiles when considering different orebody simulations, which highlights its robustness in the presence of different geological scenarios. The method can handle processes with multiple operating alternatives, additives and blending constraints. It shows substantial gains when compared with solutions generated using standard industry practices; further work includes testing the ability of the method to handle complex blending operations and assessing the possibility for further improvements based on other heuristic strategies. Figure 17 shows a cross section of the final schedule generated using the method. It can be observed that slope constraints are controlled by means of the correcting and checking mechanism utilized throughout the stages of the method.

Conclusions
An iterative improvement heuristic method is presented for generating mine production schedules in single-pit mining complexes that can contain multiple metals or attributes, multiple material types, stockpiles and processing options. The method considers relaxed capacity, availability and blending constraints. The implementation of the method in a copper deposit shows its ability to control target constraints by reducing the deviations from the capacity of Process 1 from 9 to 0.2% while increasing the expected net present value 30% when compared with a solution generated using a standard industry practice.
Regarding the expected NPV, there were no additional benefits from increasing the number of simulations after 15. This is originated from the so-called scale support effect. Although the method allows for improving an initial solution in terms of meeting target constraints and improving expected net present value, different heuristic strategies with diversification should be implemented to explore better the solution domain. Another possibility is to implement the method iteratively by considering several initial solutions simultaneously.
The possibility of adapting the method to multi-pit mining complexes is a future research avenue. Although the method requires practical amount of time for solving single-pit mining complexes (no more than 3 hours for dozens of millions of perturbations in actual-size deposits), its requirement in terms of computational time for multi-pit mining complexes needs to be addressed given the large size of the multi-pit problems. stochastic optimisation and high-order stochastic models for optimising mining complexes with uncertainty'; industry members of COSMO Laboratory: Anglogold Ashanti, Barrick Gold, BHP Billiton, De Beers, Newmont Mining, Vale and Kinross, and NSERC Discovery Grant 239019.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.