Accepted Manuscript Title: Global Optimization of Open Pit Mining Complexes with Uncertainty Author: Ryan C. Goodfellow Roussos Dimitrakopoulos PII: DOI: Reference: S1568-4946(15)00756-5 http://dx.doi.org/doi:10.1016/j.asoc.2015.11.038 ASOC 3349 To appear in: Applied Soft Computing Received date: Revised date: Accepted date: 5-3-2014 28-10-2015 21-11-2015 Please cite this article as: Ryan C. Goodfellow, Roussos Dimitrakopoulos, Global Optimization of Open Pit Mining Complexes with Uncertainty, Applied Soft Computing Journal (2015), http://dx.doi.org/10.1016/j.asoc.2015.11.038 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. Highlights:   Ac ce pt ed M an  t  A stochastic global optimization framework for open pit mining complexes is developed. The method simultaneously optimizes production schedules and downstream processes. The modeling is flexible and may be applied to numerous types of mining complexes. Three combinations of metaheuristics are tested using simulated annealing, particle swarm optimization and differential evolution. An example demonstrates that the stochastic design is 6.6% higher than the deterministic and 22.6% higher than commercial software. us cr ip  Page 1 of 53 us cr ip t Global Optimization of Open Pit Mining Complexes with Uncertainty Ryan C. Goodfellowa,∗, Roussos Dimitrakopoulosa a COSMO Laboratory, 3450 University Street, FDA Bldg., Room 123A, Montreal, QC, Canada, H3A 0E8 an Abstract Global optimization for mining complexes aims to generate a production schedule for the various mines and processing streams that maximizes the M economic value of the enterprise as a whole. Aside from the large scale of the optimization models, one of the major challenges associated with op- ed timizing mining complexes is related to the blending and non-linear geometallurgical interactions in the processing streams as materials are trans- pt formed from bulk material to refined products. This work proposes a new two-stage stochastic global optimization model for the production schedul- ce ing of open pit mining complexes with uncertainty. Three combinations of metaheuristics, including simulated annealing, particle swarm optimization Ac and differential evolution, are tested to assess the performance of the solver. Experimental results for a copper-gold mining complex demonstrate that the optimizer is capable of generating designs that reduce the risk of not meeting production targets, have 6.6% higher expected net present value than the Corresponding author Email addresses: ryan.goodfellow@mail.mcgill.ca (Ryan C. Goodfellow), roussos.dimitrakopoulos@mcgill.ca (Roussos Dimitrakopoulos) ∗ Preprint submitted to Applied Soft Computing October 28, 2015 Page 2 of 53 deterministic-equivalent design and 22.6% higher net present value than an scheduling, metaheuristics, destination policy 1. Introduction us cr ip Keywords: open pit mine design, global optimization, production t industry-standard deterministic mine planning software. Global optimization for mining complexes addresses the issue of inte- an grated mining and processing operations with multiple pits or underground mines, multiple metals or minerals, stockpiles, blending options and alternative processing streams to yield distinct products [1, 2]. The primary M objective of a mining enterprise is to maximize the net present value (NPV) of its cash flows [3, 4], which requires optimizing the long-term mine extrac- ed tion sequences and the use of the materials that have been mined. Extraction sequences define the inventories of raw materials produced from the mines. pt Downstream optimization defines how to use the mining complex’s processing streams to maximize the utility of the available materials, and addresses ce both the destination policies (where to send material from the mines) and processing stream decisions (where to send stockpiled or processed material). Ac Historically, these components have been optimized independently, leading to sub-optimal solutions for the mining complex as a whole [5]. Many of the existing attempts at global optimization ignore the compounded effect that uncertainty (i.e., geological or economic) has on the value and operational feasibility of the supply chain [6, 7]. As the complexity of the supply chain increases, with respect to the number of mines, processing stream options and methods of distribution, it becomes increasingly necessary to integrate all 2 Page 3 of 53 elements simultaneously while considering the uncertainty that arises within t the mining complex’s various components. us cr ip Recent work has focused on integrating geological, or supply, uncertainty into open pit mine production scheduling optimization models. Ramazan and Dimitrakopoulos [8] propose a two-stage stochastic integer programming (SIP) formulation [9] that seeks to maximize the NPV of a production schedule while minimizing the risk of not meeting production targets. an By using a risk discounting parameter, the optimizer aims to strike a balance between extracting high-value and low-risk material at the beginning of a mine’s life, and deferring riskier material to later periods when more M information is available. The basic SIP model has been tested and improved upon [10, 11, 12, 13, 14], and results consistently demonstrate that the NPV ed of the production schedule that considers geological uncertainty can be substantially higher than that of a conventional schedule; this is a direct result pt of managing the impact of risk to ensure production targets are attained throughout the life of the business. ce The previous formulations assume an a-priori decision of what is ore (valuable) and waste material, which is commonly referred to as a cut-off grade Ac policy [15, 16]. Some past research has attempted to integrate dynamic destination decisions with long-term deterministic production scheduling by exploiting the structure of the linear optimization model [17, 18, 19, 20]. Boland et al. [21] propose a multistage stochastic optimization model that decides the destinations for each scenario; this is overly-optimistic because it assumes perfect knowledge of the material that will be extracted and allocated to the available processing streams. Other research efforts have investigated the use 3 Page 4 of 53 of scenario-independent cut-off grade destination policies [22] and scenario- t independent block destinations [23, 24]. Each of these methods has a severe us cr ip limitation. Cut-off grade policies are primarily useful for mining complexes that treat a single metal and do not have material quality constraints in the processing streams. Using scenario-independent block destinations leads to misrepresenting the inherent operational flexibility to change block destinations based on short-term grade control information. This results in material an misclassification, where waste is sent for processing. In order to provide a useful destination policy for mining complexes in general, it is necessary to develop new approaches to define destination policies under uncertainty. M One of the underlying challenges for globally optimizing mining complexes is the inherent non-linearity related to the blending and stockpiling ed of materials [25], as well as the transformations that occur when refining the bulk input material into a set of output products. There have been several pt attempts to optimize mining complexes [1, 2, 26, 27, 28, 29, 30, 31, 32, 33], however all models ignore geological uncertainty, and are often limited in ce the degree of flexibility in modeling the non-linear transformations in the processing streams. Existing optimization approaches avoid non-linear and Ac non-additive geo-metallurgical interactions in the processing streams by preprocessing economic values, recoveries and throughput rates (among other attributes of interest) for each individual block. This assumes that each block is processed independently of others, which is often an unrealistic assumption for mining complexes that require blending and homogenization prior to processing. To provide an alternative assessment of the economic viability of the mining operation, it is necessary to explore new models that 4 Page 5 of 53 shift the focus from the economic value of the blocks that are mined to the t economic value of the products that are sold to customers, which ultimately us cr ip permits integrating non-linear processing transformations directly in the optimization model. Metaheuristics provide a useful platform for global optimization because of their ability to handle large-scale non-linear optimization models. There has been an increased effort to adapt existing metaheuristics to the open pit an mine design and production scheduling problem [14, 24, 34, 35, 36, 37, 38, 39, 40, 41]. These methods have consistently demonstrated the ability to optimize large-scale mine production schedules in a reasonable amount of time, M however, most of these works have been limited by defining the destination policy a priori, and only optimize single-tier processing stream configurations ed (i.e., from the single mine to a final processor, without stockpiles, multiple refinery options, ports or customers). pt This work focuses on the global optimization of open pit mining complexes with supply uncertainty, and addresses many limitations of previous work. ce A new type of destination policy is proposed, which is useful for decisionmaking under uncertainty and mining complexes with blending constraints Ac on secondary elements. Moreover, a new global optimization model is developed, which can be used to model multi-tier processing stream configurations with non-linear transformations. In the following section, an overview of the modeling approach is given. Subsequently, a novel two-stage SIP formulation is proposed, where the first-stage decisions optimize the extraction sequence and destination policies, and second-stage recourse decisions are used to optimize the various processing streams of the supply chain. Following this, three 5 Page 6 of 53 metaheuristic solvers are discussed, which use a combination of simulated an- t nealing, particle swarm optimization and differential evolution. The solvers conclusions and future work are presented. us cr ip are then tested and compared for a copper-gold mining complex. Finally, 2. Modeling Mining Complexes with Uncertainty This section discusses a generalized modeling framework for mining com- an plexes with uncertainty. First, definitions for materials and related attributes in the context of mining complexes with uncertainty are provided. Following this, the decision variables that govern the extraction sequence, destination M policies and processing streams are discussed. ed 2.1. Materials, flow and attributes In a mining complex, materials are products that are mined or generated through blending or processing. Materials are considered to have unique pt mineralogy or metallurgical attributes and, as a result, may only be sent to ce a set of certain locations in a mining complex for further treatment. In order to define the flow of materials from the sources (mines) to the final products (refined metals), it is useful to describe a mining complex as a directed graph. Ac Let the graph G (N , A) represent the flow of materials through the mining complex. The set of nodes, N , is comprised of three disjoint subsets: 1. C: clusters of materials at the mines that have similar attributes (e.g., metal content). See Section 2.3 for a detailed description. 2. S: destinations that are able to stockpile material over time. The input material is not treated or transformed at these nodes. 6 Page 7 of 53 3. P: destinations that process (transform) and send all output material us cr ip t to the subsequent nodes, if available. The set of directed arcs, A, defines the ability to send material from i ∈ N to a subsequent destination j ∈ S ∪ P. Let O (i) represent the set of nodes that receive material from node i, which is defined by the outgoing arcs in A from node i. Additionally, let I (j) represent the set of nodes that send materials to j, which is defined by the incoming arcs to node j in A. Let an T = {1, ..., T } describe the set of periods of time in which the mining complex operates (e.g., months, years), where T represents the end of the life for the M mining complex. In order to simplify future notations, the general case where destinations j are able to receive or produce multiple distinct materials (e.g., multiple concentrates to be sent to various smelters), or nodes i operating ed in non-contiguous periods is not discussed. The model, however, may be generalized to include these complexities. pt Attributes are used to quantify information of interest in the optimization model, such as metal quantities and costs. Uncertainty in the attributes may ce be quantified using a set of joint scenarios S = {1, ..., S}, where a scenario defines a realization, or sampling, from all sources of uncertainty. Attributes Ac are categorized into two classes: 1. Primary attributes (p ∈ P): fundamental variables of interest to the entire model (e.g., metal and total tonnages) that are sent from node i ∈ N to a node j ∈ O (i). The quantity of the attribute is denoted by vp,i,t,s ∀i ∈ S ∪ P ∪ M, t ∈ T, s ∈ S. These often originate at the mines (M), and may flow through the mining complex to the final products. The amount of attribute recovered after treatment is denoted by rp,i,t,s . 7 Page 8 of 53 2. Hereditary attributes (h ∈ H): information that is relevant to the t optimization model, but is not necessarily passed from node i to j us cr ip (e.g., feed material chemistry, processing costs, throughput rates, energy consumption and revenues from sales). These attributes may be expressed as (non-) linear functions, fh,i (vp,i,t,s), of the primary attributes p ∈ P. The quantity of a hereditary attribute is denoted by vh,i,t,s ∀i ∈ S ∪ P ∪ M, t ∈ T, s ∈ S. an 2.2. Mine extraction sequencing Mines are the suppliers of bulk materials to the mining complex, and are M represented by the set M. Each mine m ∈ M is discretized into volumes of material called blocks, Bm . To quantify the geological uncertainty for both the ed materials and attributes, it is assumed that each block b ∈ Bm has a simulated material classification and simulated attributes, βp,b,s ∀p ∈ P, s ∈ S. Fig. 1 shows a cross-section of two stochastic simulations of a real-world copper- pt gold deposit used in the case study (Section 4). It is noted that there are large differences in the material classifications and copper grades for the same ce block between the two simulations. The extraction sequence is determined by the decision variables xb,t ∈ Ac {0, 1}, which define whether (1) or not (0) block b ∈ Bm is extracted in period t ∈ T. In order to safely extract a block b ∈ Bm , it is necessary to uncover b by extracting a set of overlying blocks, Ob , in their entirety. The overlying blocks may be identified for each block b in a preprocessing step by creating an inverted cone from the center of b and verifying which blocks lie inside the cone. In more complex cases, there may be multiple geotechnical zones in a deposit, each requiring variable slope angles for the North, East, 8 Page 9 of 53 Simulated Material Types Simulated Copper Grades Simulation #2 0 Material Code 7 us cr ip t Simulation #1 0% Copper Grade 1% an Figure 1: Example of a cross-section showing a comparison of simulated material types and copper grades for the copper-gold mine used in the case study. M South and West walls. Fig. 2 gives a 2D example of how the overlying blocks are defined with variable slope angles. For a more detailed description of the ed 3D preprocessing algorithm, the reader is referred to Khalokakaie et al. [42]. ce pt + Elevation α1 + East Overlying blocks b Ac Block b α2 Figure 2: 2D example of blocks that must be uncovered (Ob ) prior to extracting block b. Note that variable slope angles are defined for the East-facing (α1 ) and West-facing (α2 ) walls. 9 Page 10 of 53 2.3. Destination policies t In this work, each material from the mines is decomposed into sub-groups us cr ip based on attributes with similar quantities (e.g., valuable or deleterious metal content). A destination policy outlines where each sub-group of material is sent for all scenarios. A similar concept is introduced by Menabde et al. [22], who separate the univariate distribution of the metal content (grades) into “bins” (i.e., categories based on ranges of metal content), and create a time- an varied cut-off grade policy based on the bins (Fig. 3A). Rather than looking at the individual blocks in the mine [21, 23, 24], the optimizer requires substantially fewer decision variables because it focuses on the distribution of M grades. In the more general case proposed herein, the sub-groupings of materials, called clusters, may be created on multivariate distributions, which ed permit a higher degree of flexibility when defining the policies (Fig. 3B). In both cases, the destination of a single block may change between simula- pt tions, depending on the simulated attributes. The general method, however, addresses many of the limitations of cut-off grade destination policies be- ce cause it is able to consider the impacts of deleterious elements, blending and stockpiling on the performance of the value chain. Ac In order to define destination policies, it is necessary to classify the simulated blocks into clusters, C ⊂ N . The k-means++ clustering algorithm [43, 44] is a useful method for grouping information with similar attributes. The algorithm first creates a predetermined number of cluster centroids for each mine and material based on the simulated attributes. For a given scenario, a block’s cluster membership is determined by the closest centroid, measured using a Euclidean distance from the block’s simulated attributes 10 Page 11 of 53 B5 B6 B7 <=B4 Waste >=B5 Mill Copper (%) Copper (%) t B4 Simulation #2 us cr ip B3 Simulation #1 Gold (ppm) B2 Probability B1 (B) All Simulations Gold (ppm) (A) Waste Stock1 Stock2 Mill Copper (%) Figure 3: (A) Simulated univariate copper distributions with a cut-off grade policy defined an using bins. (B) Destination policies using clusters, where the points represent a block’s simulated copper and gold attributes. M to the centroid’s multivariate attributes. This cluster membership is subject to change for a block, given that the material classification and attributes may vary between simulations. Let θb,c,s represent the preprocessed param- ed eter that defines whether (1) or not (0) block b ∈ Bm is a member of cluster c ∈ C in scenario s ∈ S. The destination policies are determined using the pt variable zc,j,t ∈ {0, 1}, which represents the decision of whether (1) or not (0) cluster c ∈ C is sent to destination j ∈ O (c) in period t. It is noted that the ce set of candidate destinations, O (c), is determined by the type of material that the cluster belongs to. Ac 2.4. Processing and stockpiling decisions The destination policy variables described in Section 2.3 define where to send material after it is mined. Depending on the configuration of the processing streams for a given mining complex, it may be necessary to model the transfer of materials between two locations. Processing stream decision variables, yi,j,t,s ∈ [0, 1], define the proportion of an output material sent 11 Page 12 of 53 from i ∈ S ∪ P to destination j ∈ O (i) in period t ∈ T and scenario t s ∈ S. It is noted that, unlike the extraction sequence and destination us cr ip policies, these decisions are designed to be adaptive to uncertainty; after the material is received at the initial destination and the uncertainty is revealed, it is assumed that the mining complex can adapt appropriately. Given that the primary attributes of interest are assumed to be linear and additive, the recovered quantity of the primary attribute sent from destination i ∈ S ∪ P 3. Optimizing Mining Complexes an to j ∈ O (i) is calculated using yi,j,t,s · vp,i,t,s · rp,i,t,s. M 3.1. Generalized two-stage stochastic optimization model Similar to existing deterministic models in long-term production plan- ed ning, the primary objective is to maximize the NPV [1, 3, 4, 26]. To introduce risk management as a critical objective, the proposed mathematical formula- pt tion is based on existing two-stage stochastic integer programming [9] models developed for mine planning [8, 10, 11, 13]. The proposed model discussed ce herein, however, differs substantially from existing SIP formulations because it is generalized to optimize all aspects of the mining complex holistically. Ac Both the long-term extraction sequence and destination policies are defined as first-stage decision variables, which are designed to be robust to the fluctuations that arise from the uncertainty of geological attributes. Recourse (second-stage) variables are used to adapt to the first-stage decisions, which include the processing stream variables and penalties for excessive risk or inability to meet specified targets. To determine discounted revenues and expenses, any attribute in the 12 Page 13 of 53 mining complex may be directly included in the objective function with an t associated discounted value, ph,i,t = ph,i,1 / (1 + d)t , where d is a discount rate. us cr ip − Deviation variables, d+ h,i,t and dh,i,t , are penalized using their related penalty − costs, c+ h,i,t and ch,i,t , respectively. Similar to a discount rate used to calculate the NPV, these penalty costs are be defined to be monotonically decreasing t + with respect to time, i.e., c+ h,i,t = ch,i,1 / (1 + rd) . This phenomenon, referred herein as risk discounting, attempts to defer riskier material to later periods an in the mine life. The risk discount rate, rd, is a parameter that may be used to describe the modeler’s desire to balance the ability to meet production targets in the short- and long-terms. The penalty costs relate to the will- M ingness to pay for a unit of deviation from a capacity constraint, or may be determined experimentally by running the optimization model multiple ed times to obtain a desirable risk profile for all of the constraints of interest [12]. The general global optimization formulation for open pit mining com- ce Objective: pt plexes is defined as follows: Ac max X XXX 1 ph,i,t · vh,i,t,s |S| i∈S∪P∪M t∈T h∈H s∈S | {z } Discounted revenues and costs X XXX
1 + − − c+ − h,i,t · dh,i,t,s + ch,i,t · dh,i,t,s |S| i∈S∪P∪M t∈T h∈H s∈S | {z } (1) Risk-discounted penalties for deviations Subject to: I. Capacity constraints are used to calculate the value of the hereditary attributes and evaluate the surplus or shortage from an upper- (Uh,i,t ) or 13 Page 14 of 53 lower-bound (Lh,i,t ) at the various locations in the mining complex. Examples t of typical constraints may include, but are not limited to, mine production us cr ip capacity, stockpile capacity, processing hours and grade blending constraints. ∀h ∈ H, i ∈ S ∪ P ∪ M, t ∈ T, s ∈ S vh,i,t,s = fh,i (vp,i,t,s) (2) vh,i,t,s − d+ h,i,t,s ≤ Uh,i,t ∀h ∈ H, i ∈ S ∪ P ∪ M, t ∈ T, s ∈ S (3) vh,i,t,s + d− h,i,t,s ≥ Lh,i,t ∀h ∈ H, i ∈ S ∪ P ∪ M, t ∈ T, s ∈ S (4) an II. Reserve and block access constraints ensure a block is extracted at (5) ∀b ∈ Bm , u ∈ Ob , t ∈ T (6) most once, and that the overlying blocks have also been extracted. xb,t ≤ 1 ∀b ∈ Bm t∈T t X xu,t′ ed xb,t ≤ M X t′ =1 pt III. Destination policy constraints calculate the quantities of the primary attributes for each cluster from the blocks, and ensure that a cluster is sent ce to only one destination in the mining complex. γp,c,t,s = X θb,c,s · βp,b,s · xb,t ∀m ∈ M, p ∈ P, c ∈ C, s ∈ S (7) Ac b∈Bm X zc,j,t = 1 ∀c ∈ C, t ∈ T (8) j∈O(c) IV. Mine extraction constraints are used to determine the quantities of the primary attributes that are extracted from each mine (e.g., annual tonnage). vp,m,t,s = X βp,b,s · xb,t ∀m ∈ M, p ∈ P, t ∈ T, s ∈ S (9) b∈Bm 14 Page 15 of 53 V. Processing stream flow constraints calculate the quantity of the pri- t mary attributes that are retained or received at each location, and ensure sequent destinations in the mining complex. vp,j,(t+1),s = X rp,i,t,s · vp,i,t,s · yi,j,t,s + i∈(I(j)\C) us cr ip mass balancing for materials sent from the stockpiles and processors to sub- X γp,c,(t+1),s · zc,j,(t+1) c∈(I(j)∩C)  yj,k,t,s ∀p ∈ P, j ∈ S ∪ P, t ∈ T, s ∈ S (10) an k∈O(j)  (11) M +vp,j,t,s · 1 − X ∀p ∈ P, i ∈ P, t ∈ T, s ∈ S (12) yi,j,t,s ≤ 1 ∀i ∈ S, t ∈ T, s ∈ S (13) yi,j,t,s = 1 ∀i ∈ P, t ∈ T, s ∈ S (14) rp,i,t,s = 1 ∀p ∈ P, i ∈ S, t ∈ T, s ∈ S X ed rp,i,t,s = fh,i (vp,i,t,s) j∈O(i) pt X j∈O(i) ce VI. End-of-year stockpile quantities are used to calculate the quantities Ac of materials that remain in the stockpile at the end of the production period.   X vh,i,t,s = vp,i,t,s · 1 − yi,j,t,s ∀i ∈ S, t ∈ T, s ∈ S (15) j∈O(i) VII. Binary constraints for the extraction sequence (xb,t ) enure that blocks are not partially mined in any given period. Binary constraints for the destination policies variables (zc,j,t) ensure that a cluster c is only sent to a single destination j ∈ O (c) in period t. xb,t ∈ {0, 1} ∀b ∈ Bm , t ∈ T (16) 15 Page 16 of 53 zc,j,t ∈ {0, 1} ∀c ∈ C, j ∈ O (c) , t ∈ T ∀p ∈ P, c ∈ C, t ∈ T, s ∈ S (18) ∀i ∈ S ∪ P, j ∈ O (i) , t ∈ T, s ∈ S (19) ∀p ∈ Bm , i ∈ S ∪ P, t ∈ T, s ∈ S (20) vp,i,t,s ∈ R ∀p ∈ P, i ∈ S ∪ P ∪ M, t ∈ T, s ∈ S (21) vh,i,t,s ∈ R ∀h ∈ H, i ∈ S ∪ P ∪ M, t ∈ T, s ∈ S (22) rp,i,t,s ∈ [0, 1] ∀h ∈ H, i ∈ S ∪ P, t ∈ T, s ∈ S (23) M − d+ h,i,t,s , dh,i,t,s ≥ 0 an yi,j,t,s ∈ [0, 1] us cr ip γp,c,t,s ∈ R t VIII. Other variable definitions (17) 3.2. Optimization using metaheuristics ed The generalized global optimization formulation can be challenging to solve using conventional mathematical programming methods, particularly when the models include non-linear functions. Metaheuristics are algorithms pt that are useful for such cases because they do not require linear formulations or a special structure in the optimization problem. Metaheuristics do not ce guarantee a mathematically optimal solution, but have provided useful solutions for mining-related problems [14, 24, 34, 36, 37, 38, 39, 40, 41, 45, 46]. Ac These models and solution methods, however, have not been extended to the global optimization of mining complexes. Simulated annealing [47, 48, 49] is selected as a base algorithm for com- parison purposes because of previous success using the method for extraction sequencing [24, 34, 37, 38]. While simulated annealing is often used for discrete optimization problems, the algorithm has been adapted to accommodate the continuous processing stream decision variables (yi,j,t,s ) [50, 51]. 16 Page 17 of 53 One of the primary concerns with applying a single metaheuristic is the pos- t sibility of getting trapped in local optima, particularly in the context of a us cr ip mining complex with three distinct sets of decision variables that are inter- related. The ability to optimize both discrete and continuous variables is an important feature when selecting metaheuristics for the model at hand. Both particle swarm optimization (PSO) [52, 53] and differential evolution (DE) [54, 55] are known for their ability to optimize mixed integer non-linear an optimization models, without the need for complex encoding and decoding schemes. PSO has been used for mine extraction sequencing [39] by encoding the M extraction sequence as a two-dimensional surface. The decision variables relate to the depth mined at an (x, y) coordinate pair (column) in a given period ed [36, 41]. While this encoding scheme is particularly suitable for populationbased metaheuristics, initial tests for the proposed global optimization model pt using purely PSO or DE indicate that the method is sensitive to the initial sequences and destination policies generated for the population. Addition- ce ally, these methods require a substantial computational effort in order to normalize the extraction sequence and enforce the slope constraints (Eq. Ac (6)). These methods, however, remain appealing, particularly for optimizing the destination policies and processing streams (i.e., downstream optimization), given that they simultaneously modify both sets of interrelated decision variables. To assess the performance of the basic simulated annealing (SA) algorithm, two additional variants are tested: simulated annealing with downstream particle swarm optimization (SA-D-PSO), and simulated annealing with downstream differential evolution (SA-D-DE). Algorithm 1 17 Page 18 of 53 (see Appendix) provides an overview of the global optimization metaheuris- t tic developed. Computational results for a copper-gold mining complex are us cr ip discussed in Section 4.2. A solution vector, Φ = {x, z, y}, is used to store all decision variables. The extraction sequence vector (x) stores an encoded version of all xb,t variables, where a discrete-valued element xb ∈ x represents the extraction period of block b ∈ Bm , and may take on any value in T ∪ {T + 1}, where an extrac- an tion period of T + 1 inidcates that the block is not mined. The destination policy vector (z) stores an encoded version of the zc,j,t variables, where each element zc,t ∈ z represents the encoded destination for cluster c ∈ C in period M t ∈ T, and may be in the range [1, ..., |O (c) |]. A decoding scheme is used to convert the value to the appropriate destination for the cluster. Finally, ed the processing stream vector (y) stores all yi,j,t,s variables, where each element y ∈ y maps directly to a yi,j,t,s variable, and may lie in the range [0, 1]. pt To initialize the algorithm, an initial sequence (x) may be used, or may be initialized with all blocks un-mined; the destination policies and processing ce stream variables are randomly generated within their bounds. It is noted that Equations (3) and (4) are inherently modeled as soft con- Ac straints in the original mathematical formulation; in these cases, deviations from the capacity constraints are penalized in the objective function. While other deterministic models solved with metaheuristics [36, 39, 41] treat penalties as a tool to eliminate constraint violations, penalties are used as a tool to manage the distribution of risk through time in the stochastic approach [8, 10, 11, 12, 13, 14]. The slope constraints in Eq. (6) are guaranteed during the perturbation mechanism of the algorithm [37], and the reserve constraints 18 Page 19 of 53 in Eq. (5) are guaranteed through the encoding scheme (x). Similarly, the t single-destination constraints in Eq. (8) are inherent in the encoding scheme us cr ip (z). Finally, the processing stream proportion constraints from Eq. (13) and Eq. (14) are guaranteed by normalizing the values after the optimizer makes any changes. It follows that all other constraints are automatically feasible because they are simply calculations derived from other variables. The basic algorithm uses a modified SA algorithm (Algorithm 2) to im- an prove the global best solution vector, Φg . Three classes of perturbations (solution changes) are considered during annealing: M 1. Extraction sequence perturbations (x ∈ Φ): a block is randomly selected, and its mining period is changed (possibly to not being mined at all). Blocks that would violate the slope constraints as a result of ed this change are also considered [37]. 2. Destination policy perturbations (z ∈ Φ): a cluster destination decision possible. pt variable is randomly selected and sent to a different destination, if ce 3. Processing stream perturbations (y ∈ Φ): a processing stream variable is randomly selected and the value is modified using a random normal Ac number, N (yi,j,t,s, 0.1). Note that the variance of the normal distribution is sufficiently small to permit both local and global exploration. All processing stream variables are then normalized to respect Eq. (13) and Eq. (14). The perturbation mechanism is outlined in Algorithm 3. Two parameters, probseq and probdest , are used to determine the probability of selecting a 19 Page 20 of 53 sequencing, destination policy perturbation or processing stream perturba- t tion. The probability of accepting a perturbation of the solution vector for us cr ip a maximization problem is based on the following distribution:  1 if g (Φ′ ) ≥ g (Φ) ′ (24) P (g (Φ) , g (Φ ) , δ) = exp (− |g (Φ′ ) − g (Φ)| /δ) otherwise where g (Φ) and g (Φ′ ) are the objective function values before and after the perturbation, respectively, and δ is the annealing temperature. As the an algorithm progresses, the temperature is gradually reduced until only minor changes are accepted. This is often controlled by the initial temperature M at the start of the algorithm, δ (0), the cooling schedule, which is defined by a reduction factor, k ∈ [0, 1), and a number of iterations before the reduction factor is applied, niter . Fig. 4 shows a comparison of the cumulative ed distributions between the changes in objective function values (g (Φ)−g (Φ′ )) for non-improving perturbations for the three neighborhoods. When using a pt single temperature in the classic SA algorithm, for a very large temperature, the optimizer may limit the number of sub-optimal changes in the destination ce policy, but will likely accept all extraction sequence and processing stream changes. As the temperature decreases, it becomes more likely that only Ac non-improving processing stream changes are accepted. This behavior is not desirable, given that the neighborhoods are strongly related. In the proposed modified SA algorithm, the cumulative probability dis- tributions (cdfseq , cdfdest and cdfproc), shown in Fig. 4, are constructed by proposing random perturbations to the solution. Rather than using a single temperature, δ, for the neighborhoods in Eq. (24), the optimizer uses a single parameter, ρ ∈ [0, 1], which represents a probability of accepting a 20 Page 21 of 53 1 0.9 t 0.7 0.6 0.5 0.4 0.3 Extraction sequence 0.2 Destination policy 0.1 Processing stream 0 0.E+00 4.E+06 an Sub-optimal Objective Function Change g(Φ)-g(Φ') us cr ip Acceptance Probability 0.8 Figure 4: Example of cumulative probability distributions of objective function changes M for perturbations that do not improve the solution quality. non-improving perturbation. For a fixed ρ, the correct temperature variable ed δ is retrieved from the appropriate cumulative probability distribution using −1 −1 −1 a look-up function (cdfseq (ρ), cdfdest (ρ) or cdfproc (ρ)). The cooling schedule pt (k, niter ) is then applied to ρ, rather than δ. As the algorithm progresses, the information garnered from any new proposed non-improving perturba- ce tions is used as feedback to update the cumulative distributions; this better reflects the current search space, rather than the search space when the SA Ac algorithm commenced. In the global optimization algorithm (Algorithm 1), once SA is complete, the acceptance probability ρ is reset to its initial value, and the algorithm is re-started from the global best solution Φg . This diversification strategy is repeated until no improvement is found. In order to assess the performance of the basic SA algorithm, two variations (PSO and DE) are proposed to improve the solution of the downstream variables (z, y ∈ Φg ) prior to di21 Page 22 of 53 versifying with annealing. An outline of the downstream optimization is t provided in Algorithm 4. Both PSO and DE are population-based meta- us cr ip heuristics that employ unique approaches to modify existing solutions. Both approaches require defining the size of the population, NP , as a parameter. In PSO, a member of the population (particle) is comprised of three equally-sized vectors: its solution vector, Φq , its best solution vector, Φbest q , and its velocity vector, vq . Initially, the vectors Φq and vq are randomly an generated. To maintain a high-quality solution in the swarm and reduce the chance of premature convergence, the last particle’s best vector, Φbest N P , is initialized to the global best solution previously found from annealing (Φg ). M Algorithm 5 provides an outline of how the particle’s vectors are updated. The key parameters that need to be calibrated for PSO are c1 , c2 and c3 , ed which represent the weights for the particle’s inertia, its best solution and the local-best solution, respectively. It is noted that in this algorithm, a localbest ring topology [56] using the nearest NP local particles on either side of ce pt particle q is used. The local-best solution, Φbest lbest , represents the best solution  of the particles with adjacent indices, i.e., q − NP local , q + NP local , and is determined experimentally to delay swarm convergence. Ac The DE algorithm operates similar to the PSO algorithm. The population is represented by a set of agents. Algorithm 6 provides an outline of how an agent’s solution vector is modified. In order to update an agent q, three other agents (a, b and c) are randomly selected from the population. An agent’s best solution vector, Φbest q , is modified by randomly selecting variables to be crossed over with agent a and a weighted difference between agents b and c. While there are several types of variations that may be used to modify 22 Page 23 of 53 the solution [57, 58], it has been determined experimentally that a simple t modification structure provides adequate results. This algorithm requires us cr ip calibrating two parameters: CR, a cross-over ratio that defines a probability for modifying a decision variable, and F , a parameter that defines the weight applied to a difference between the two solution vectors, b and c. Unlike PSO, DE only modifies a proportion of the decision variables for each agent in each iteration; it has been found through testing that this helps to prevent an the population from converging prematurely. 4. Case Study – Application at a Copper-Gold Mining Complex M The proposed integrated mine planning optimization framework is demonstrated on a real-world copper-gold mining complex. ed 4.1. Overview of the mining complex In the given case study, a single mine supplies materials to a mining pt complex that produces gold and copper. Fig. 5 summarizes the definition of the mine’s materials and the processing options. The mine contains three ce main material groups: sulfides, transition and oxides. In order to respect the chemistry requirements at the sulfide heap leach (processor), the sulfide Ac and transition material groups are both separated into two different material types based on being above or below 0.2% copper. The oxide materials are classified as ore or waste based on chemistry. The deposit’s uncertainty is represented by a set of 25 equally probable geological simulations with variable copper, gold, tonnages and material types. With the exception of the oxide waste dump, all destinations (processors) have variable grade-recovery curves that are based on the average grade of 23 Page 24 of 53 (1 Mtpa) Sulfide Mill (3 Mtpa) Copper-Gold Mine (8 Mtpa) Sulfides us cr ip Sulfide Heap Leach (25 Mtpa) <0.2% Cu Sulfide Waste Dump ≥ 0.2% Cu (Unlimited) Transition Heap Leach Transition <0.2% Cu (Unlimited) ≥ 0.2% Cu Oxide Heap Leach Oxides (Unlimited) Potential Cu Au t Stockpile Cu Cu Au Au Oxide Waste an Waste (Unlimited) M Figure 5: Definition of material types at the copper-gold mine, along with the various destinations. ed the incoming material at a process in a given period (Fig. 6). The non-linear grade-recovery relationships have interesting implications when considering pt the transition materials: for a given block or cluster that has (hypothetically) similar economic values for two processing options, the selected destination ce is the one that profits the most from an increase in recovery. As a result, one cannot assume that the destinations can be specified a-priori in a greedy Ac manner because it is the recovery of the aggregated material sent to a given processor that determines the potential value. For this reason, the model never considers the economic value of a block; the economic values are evaluated using only the recovered copper and gold at the processors (mill and leach pads). The objective for the optimizer is to maximize the NPV of the mining complex, which considers the sale of copper and gold from each of the desti- 24 Page 25 of 53 100% 90% 90% 70% 60% 50% 40% 30% Sulfide Heap Leach 20% 80% t 80% 70% 60% 50% 40% 30% Transition Heap Leach 20% Sulfide Dump Leach Oxide Heap Leach 10% 10% Sulfide Mill Sulfide Mill 0% 0% 0 0.5 1 us cr ip Au Recovery (% Relative to Max) Cu Recovery (% Relative to Max) 100% 1.5 0 Copper Head Grade (% Cu) 0.5 1 1.5 Gold Head Grade (gpt) an Figure 6: Grade-recovery curves for copper (left) and gold (right) at each of the processors. nations, along with the processing and mining costs. All cost-related param- M eters in Table 1 are expressed relative to the mining cost for confidentiality purposes. Table 2 summarizes the constraints and penalty costs used for the ed deterministic and stochastic optimization models. A risk discount rate [8] of 10% is used to penalize the deviations from the production capacities, and pt ensures that riskier material is deferred to later periods when more geological information is available. The mine model contains 34 057 blocks that may ce be scheduled over 22 years, and a slope angle of 45◦ is used [42]. 4.2. Deterministic optimization and risk analysis Ac A deterministic orebody model is created by averaging the grades in the given simulations and re-classifying the material types in the same manner that the simulations are classified. This model is herein referred to as the “E-type orebody model” [59]. A base-case design is first obtained using the E-type orebody model and an industry-standard commercial mine planning software package. The results from optimizing the deterministic model with the proposed optimizer will be referred to as the “deterministic-equivalent” 25 Page 26 of 53 Value $1.00/t Sulfide mill cost* $11.30/t Sulfide heap leach cost* $2.98/t Sulfide dump leach cost* $1.87/t Transition heap leach cost* $2.15/t Oxide heap leach cost* $2.06/t Copper price $2.88/lb Gold price $1480/oz 7% 10% ed Risk discount rate (rd) M Economic discount rate (d) an Mining cost* us cr ip Parameter t Table 1: Economic parameters used in the model. * For confidentiality reasons, all costs are ex- pt pressed relative to the mining cost. ce Table 2: Lower and upper bounds (Lh,i,1 and Uh,i,1 ) and respective penalty costs (c− h,i,1 and c+ h,i,1 ) used in the deterministic and stochastic optimization models. Ac Description of constraint (Mtpa) Lower, Upper Bounds (*106) Penalty ($/tonne) Mine capacity -, 25.0 -, 10 Stockpile capacity -, 1.0 -, 20 Sulfide mill capacity 2.8 a,b , 3.0 50, 50 Sulfide heap leach capacity 7.8 a,b , 8.0 10, 25 a Lower-bounds are not enforced for the deterministic model. b Lower-bounds are removed from periods 11-22. 26 Page 27 of 53 design. For the deterministic-equivalent model, each material is clustered t into 15 groups, with the exception of the oxide waste material, which only us cr ip considers one. As a result, there are 1 672 encoded destination policy vari- ables (z), along with 34 057 encoded block extraction decisions (x) and 22 stockpile processing stream variables (y). Table 3 summarizes the parameters used for the simulated annealing, particle swarm optimization and differential evolution algorithms. It is noted an that in the first iteration of SA (igopt = 0, Algorithm 1), a different set of parameters is used to compensate for the large number of new blocks that are introduced into the extraction sequence. Table 4 shows a comparison of M the SA, SA-D-PSO and SA-D-DE algorithms, where 10 trials are performed for each method. All trials are run on Amazon’s EC2 c4.8xlarge virtual ed machines with 36 virtual CPUs (Intel Xeon E5-2666 v3) and 60 GB RAM. Each trial is initialized without the use of an initial extraction sequence pt (i.e., all xb,t = 0), and the downstream solution vectors (z and y) have been randomly generated. It is apparent that the SA-D-PSO and SA-D-DE ce algorithms consistently outperform the algorithm that only uses simulated annealing, with an average increase in net present value (NPV) of 1.91% Ac and 2.57%, respectively. It is noted that the final solutions do not contain capacity constraint violations, thus the NPV is analogous to the objective function value. These slight improvements, however, come with a drastically longer computing time. On average, it takes over 2.4 and 2.9 times longer to run the SA-D-PSO and SA-D-DE algorithms, respectively, than the basic SA algorithm. The basic SA algorithm converges after an average of 14.7 diversifications. When SA is combined with another metaheuristic to optimize 27 Page 28 of 53 the downstream variables, the number of diversification iterations increases. t Despite the SA-D-DE algorithm having a smaller initial population than SA- us cr ip D-PSO (16 agents versus 30 particles), the average run time is higher due to the fact that more diversifications are performed, and the SA algorithm requires more iterations to catch up and improve upon the global best solution after optimizing with DE. Fig. 7 (left) shows a comparison for the sulfide mill, the sulfide heap an leach and the cumulative discounted cash flows for the commercial and deterministic-equivalent extraction sequences. It is noted that the deterministicequivalent solution generates a 13.8% higher NPV than the solution gen- M erated by the commercial mine planning tool, and the life of the mine is extended by a year. This drastic increase is largely a result of the global ed optimization approach, whereby the production schedule, destination policy and stockpiles are optimized simultaneously, rather than the commercial pt approach, which relies on an iterative approximation (ultimate pit limit and phase design definition prior to scheduling) and heuristics in order to provide ce a production schedule. A risk analysis is performed by taking the schedule and destination poli- Ac cies generated from the deterministic optimizer and testing how the 25 geological simulations perform. It is noted that the extraction sequence and destination policies remain fixed, however, the stockpile variables (y) are optimized for each scenario using the DE algorithm. Fig. 7 (left) also provides the risk profiles, which are defined by the exceedance probabilities for 10%, 50% and 90% (P90, P50, P10, respectively) of the simulations’ responses to the deterministic-equivalent design. The simulations indicate that the deter- 28 Page 29 of 53 us cr ip t Table 3: Simulated annealing, particle swarm optimization and differential evolution parameters used for both the deterministic and stochastic optimization trials. Parameter Value Initial annealing acceptance probability (ρ) Cooling factor (k) 0.4∗ , 0.3 0.99 Cooling iterations (niter ) 600 an Perturbation probability - extraction sequence (probseq ) Perturbation probability - destination policy (probdest ) Annealing global best updates before diversification (igbu ) Number of particles (NP ) M Total annealing iterations before diversification (itotal ) 0.3∗ , 0.5 0.6∗ , 0.4 2 000 500 000 30 ±1 Particle inertia (c1 ) 0.9 ed Number of particles in neighborhood (NP local ) pt Particle best inertia (c2 ) Local best inertia (c3 ) 0.8 0.6 ce Population convergence criteria (pcc) Maximum number of iterations (itotal ) 2 500 16 Ac Number of agents (NA) 0.0025 Cross-over ratio (CR) 0.2 Weighting factor (F ) 0.8 Population convergence criteria (pcc) 0.0025 Maximum number of iterations (itotal ) 2 500 * Parameters used for the first SA iteration to compensate for a large number of blocks introduced into the extraction sequence. 29 Page 30 of 53 t us cr ip Table 4: Performance comparisons for deterministic optimization using the simulated annealing (SA), simulated annealing with downstream particle swarm optimization (SA- an D-PSO) and simulated annealing with downstream differential evolution (SA-D-DE) algorithms. SA Number of trials M Parameter SA-D-PSO SA-D-DE 10 10 10 1.085 1.105 1.104 Average relative NPV* 1.094 1.115 1.122 Maximum relative NPV* 1.104 1.123 1.138 Relative NPV standard deviation* 0.007 0.006 0.009 Average computing time (h) 5.345 13.033 15.499 Computing time standard deviation (h) 1.57 1.396 3.113 Average number of diversifications 14.7 17.8 18.3 ce pt ed Minimum relative NPV* * For confidentiality purposes, NPVs are expressed relative to the base-case solution Ac generated using an industry-standard commercial mine planning tool. 30 Page 31 of 53 ministic design performs noticeably worse than the E-type orebody model t indicates, which is attributed to the differences in the univariate grade dis- us cr ip tributions between the estimated model and the simulations. The spikes in sulfide mill tonnage above the capacity are in the order of 5-10%, and are related to the fact that the simulations contain excessive quantities of high-grade material in the production schedule; these tonnage spikes are unrealistic and inflate the NPV of the risk analysis. Moreover, the risk analysis an indicates that the deterministic design is consistently unable to meet the production targets at the sulfide heap leach. Interestingly, the NPV of the risk analysis is unaffected by its inability to fill the sulfide heap leach to capacity. M This is a result of sending excessive and higher valued materials to the sulfide mill, which leads to an increased recovery and profits from the sulfide mill, ed despite the underutilization of the sulfide heap leach. 4.3. Stochastic optimization pt The stochastic optimizer aims to optimize the extraction sequence, destination policies and the use of the stockpile while considering the geological ce uncertainty. Table 2 outlines the penalty costs and constraint bounds that are used to guide the risk profiles. A total of 10 trials are run for the stochas- Ac tic example using the SA-D-DE optimizer. The average computing time for the stochastic optimizer is 75.6 hours, with a standard deviation of 8.6 hours. When comparing the objective function value of the best stochastic design to the simulations through the best deterministic design (calculated using the same penalty costs and constraint bounds), the best stochastic design has a 16.6% higher value, with an average of 15% across all 10 trials. This is a direct result of being able to manage the impact of risk on the sulfide mill 31 Page 32 of 53 and sulfide heap leach tonnages over time. t To elucidate these results, Fig. 7 (right) shows the risk profiles of the us cr ip stochastic solution for the sulfide mill tonnages, sulfide heap leach tonnages and the cumulative discounted cash flows. It is apparent that the stochastic design is capable of meeting the target capacity of 3 Mtpa at the sulfide mill on average, with substantially less risk than the deterministic design (Fig. 7, left). The stochastic design also meets the target sulfide heap leach tonnage an over the first 12 periods, after which, the quantities begin to decline. This implies that the riskier sulfide heap leach material is being deferred until the end of the mine life, which is a result of incorporating risk discounting in M the stochastic optimization model. Finally, it is noted that the NPV of the stochastic design is 6.6% higher than the deterministic-equivalent design, ed when comparing the P-50 values, and 22.6% higher than the commercial Ac ce pt (deterministic) mine planning software. 32 Page 33 of 53 Sulfide Mill Tonnage Deterministic Design Sulfide Mill Tonnage Stochastic Design 5.0 4.0 Deterministic-Equiv. 3.5 Commercial 3.0 2.5 2.0 1.5 1.0 0.5 0.0 P-10/90 4.5 t P-50 P-50 4.0 3.5 3.0 2.5 2.0 1.5 1.0 0.5 0.0 0 5 10 15 20 25 0 Period 6.5 6.0 P-10/90 P-50 Deterministic-Equiv. 4.0 0 5 10 15 20 8.0 7.5 7.0 6.5 6.0 5.5 P-10/90 5.0 P-50 4.5 4.0 0 25 5 10 80% 60% 40% 20% 0% 0 5 P-10/90 P-50 Deterministic-Equiv. Commercial 10 15 20 25 15 20 25 Period Relative Cumulative DCF (Compared to Commercial Software) pt ce 100% Ac Relative Cumulative DCF (Compared to Commercial Software) Cumulative Discounted Cash Flows Deterministic Design 120% 25 8.5 Period 140% 20 9.0 M 7.0 Commercial 15 an 7.5 Sulphide Heap Leach Tonnage (Mt) 8.0 ed Sulphide Heap Leach Tonnage (Mt) 8.5 4.5 10 Sulfide Heap Leach Tonnage Stochastic Deisgn 9.0 5.0 5 Period Sulfide Heap Leach Tonnage Deterministic Design 5.5 us cr ip P-10/90 4.5 Sulphide Mill Tonnage (Mt) Sulphide Mill Tonnage (Mt) 5.0 Cumulative Discounted Cash Flows Stochastic Design 140% 120% 100% 80% P-10/90 60% P-50 40% 20% 0% 0 5 10 15 20 25 Period Period Figure 7: (Left) Comparison of the deterministic design generated using a commercial mine planning package, the best deterministic solution generated and its associated risk profiles. (Right) Risk profiles of the best stochastic design. 33 Page 34 of 53 5. Conclusions us cr ip t This work presents a framework for global asset optimization of mining complexes under uncertainty, whereby the solutions provide robust long-term open-pit mine extraction sequences and destination policies. The proposed framework permits a high degree of flexibility and detail in modeling a mining complex, including the opportunity to integrate non-linear relationships that are generally ignored in existing models because of the challenges associated an with non-linear optimization. As a result, the proposed method overcomes a major limitation of existing global optimization models: modeling the eco- M nomic value of the products sold, rather than the value of the blocks mined. A new form of destination policies is also developed that overcomes the severe limitations of other methods, such as quality constraints on secondary ed elements in the processing streams. Three solvers, based on a combination of simulated annealing, particle swarm optimization and differential evolution, pt are developed and compared. The method is tested on a copper-gold mining complex. When optimiz- ce ing the deterministic example, the global optimizer achieves a design with 13.8% higher net present value than an industry-standard mine planning Ac tool. While optimizing with both simulated annealing and differential evolution demonstrated improved results, this comes at the expense of a substantially longer computing time. For the stochastic optimization, results indicate that the stochastic design is able to satisfy the target tonnage capacities, thus ensuring that the mine is able treat profitable material over the life of mine. Additionally, the stochastic solution indicates a 6.6% higher net present value when compared to the deterministic-equivalent design and 34 Page 35 of 53 22.6% higher value than the industry-standard software. Given that the pro- t posed method seeks to generate a single set of destination policies, future us cr ip research will investigate the use of multistage stochastic optimization in order to permit adaptive policies under both supply (geological) and demand (metal price) uncertainty, which will likely lead to higher economic value. Acknowledgements an The work in this paper was funded from NSERC CDR Grant 335696 and BHP Billiton, NSERC Discovery grant 239019 and the members of the COSMO Stochastic Mine Planning Laboratory: AngloGold Ashanti, Barrick M Gold, BHP Billiton, De Beers, Newmont Mining, and Vale. ed Appendix A. Metaheuristic Algorithms for the Global Optimiza- Ac ce pt tion of Mining Complexes with Uncertainty 35 Page 36 of 53 t us cr ip Algorithm 1 Algorithm for the global optimization of mining complexes. Require: Φ = {x, y, z} ⊲ Initial solution vector. ⊲ T rue if particle swarm optimization will be used. an useP SO useDE ⊲ T rue if differential evolution will be used. M procedure GlobalOptimization(Φ, useP SO, useDE) Φg ← Φ ⊲ Initialize the global best solution vector. igopt ← 0 while true do ed ⊲ Global optimization iteration counter. igopt ← igopt + 1 ⊲ Copy the solution vector for later comparison. pt Φ ← Φg ⊲ Update the iteration counter. ⊲ See Algorithm 2. if useP SO = true or useDE = true then ⊲ See Algorithm 4. ce Φg ←SimulatedAnnealing(Φg ) Φg ← DownstreamOptimization(Φg , useP SO, useDE) Ac if g (Φg ) = g (Φ) then break ⊲ No update in objective function. 36 Page 37 of 53 Algorithm 2 Simulated annealing for open pit mining complexes. Require: ρ, k, niter probseq , probdest igbu us cr ip t ⊲ Acceptance probability, cooling factor, cooling schedule. ⊲ Probability of a schedule or destination perturbation. ⊲ Total number of global best updates before terminating. itotal ⊲ Total number of iterations before terminating. cdfseq , cdfdest , cdfproc ⊲ Distributions of non-improving scheduling, destination policy and processing stream changes. an function SimulatedAnnealing(Φg ) Φ, Φ′ ← Φg ⊲ The current and perturbed solution vectors. i, iu ← 0 M ⊲ Annealing iteration and global update counters. δ←0 ⊲ A temporary annealing temperature. i←i+1 ed while true do ⊲ Update the iteration counter. if i mod niter = 0 then pt ρ← ρ∗k ⊲ Update the acceptance probability. Φ′ , δ ← PerturbSolution(Φ, ρ) ce r ← U [0, 1] ⊲ See Algorithm 3. ⊲ A random uniform number. if P (g (Φ) , g (Φ′ ) , δ) ≥ r then Φ ← Φ′ ⊲ See Eq. (24). Ac ⊲ Update the current solution. if g (Φ′ ) < g (Φ) then Update cdfseq , cdfdest or cdfproc with |g (Φ′ ) − g (Φ) | if g (Φ′ ) ≥ g (Φg ) then Φg ← Φ′ iu ← iu + 1 ⊲ Update the global best solution. ⊲ Increase the global best update counter. if i = itotal or iu = igbu then break 37 return Φg Page 38 of 53 t Algorithm 3 Perturb the current solution and the retrieve the appropriate us cr ip annealing temperature with respect to the selected perturbation neighborhood. function PerturbSolution(Φ, ρ) r ← U [0, 1] ⊲ A random uniform number. δ←0 ⊲ A temporary annealing temperature. if r < probseq then an Randomly select xb , an encoded variable from x ∈ Φ. Find the set of blocks xb (t′ ) that must be extracted in t′ to satisfy Φ′ ← [x ⊕ xb (t′ ) , z, y] −1 δ ← cdfseq (ρ) M Eq. (6), see [37]. ⊲ Update the solution vector. ed ⊲ Retrieve the sequencing temperature.
else if r < probseq + probdest then Randomly select zc,t , an encoded variable from z ∈ Φ. ⊲ Get a new (encoded) cluster destination. pt ′ zc,t ← U [0, |O (c) |]   ′ Φ′ ← x, z ⊕ zc,t ,y ce −1 δ ← cdfdest (ρ) ⊲ Update the solution vector. ⊲ Retrieve the destination policy temperature. else Ac Randomly select yi,j,t,s, a variable from y ∈ Φ. ′ yi,j,t,s ← yi,j,t,s + N (yi,j,t,s , 0.1)   ′ Φ′ ← x, z, y ⊕ yi,j,t,s ⊲ Add random normal number. ⊲ Update the solution vector. Normalize y ∈ Φa to obey Eq. (13) and Eq. (14). −1 δ ← cdfproc (ρ) return Φ′ , δ ⊲ Retrieve the processing stream temperature. ⊲ Return the new solution vector and temperature. 38 Page 39 of 53 t Algorithm 4 Downstream optimization using particle swarm optimization us cr ip or differential evolution. Require: NP ⊲ Size of the population. NP local ⊲ Number of particles in the local neighborhood (PSO). c1 , c2 , c3 ⊲ Inertia, personal best and local-best coefficients (PSO). CR, F ⊲ Crossover ratio and weighting factor (DE). ⊲ Population convergence criteria. an pcc itotal ⊲ Total number of iterations before terminating. for all q ∈ {1, ..., NP } do M function DownstreamOptimization(Φg , useP SO, useDE) ed Randomize Φq (PSO, DE) and velocity vq (PSO). (x ∈ Φq ) ← (x ∈ Φg ) ⊲ Extraction sequence remains constant. Normalize y ∈ Φq to obey Eq. (13) and Eq. (14). ⊲ Initialize the member’s best solution vector. pt Φbest ← Φq q ⊲ Copy the global best solution to the last member. ce g Φbest NP ← Φ while true do i←i+1 ⊲ Update the iteration counter. Ac for all q ∈ {1, ..., NP } do if useP SO = true then local Get Φbest with the best lbest , the member within q ± NP objective function. best Φq ← PSOUpdate(vq , Φq , Φbest q , Φlbest ) 39 Page 40 of 53 else if useDE = true then ⊲ See Algorithm 6. Correct z ∈ Φq to obey Eq. (8). us cr ip best best best Φq ← DECrossover(Φbest q , Φa , Φb , Φc ) t Randomly select members a, b and c, where (a, b, c) 6= q. Normalize y ∈ Φq to obey Eq. (13) and Eq. (14).
if g (Φq ) ≥ g Φbest then q Φbest ← Φq q if g (Φq ) ≥ g (Φg ) then an ⊲ Update member’s best solution. M Φg ← Φq ⊲ Update the global best solution.
P P best g Φ ⊲ Get the average objective value. avg ← N1P · N q q=1

if i = itotal or g Φbest − avg /avg < pcc ∀q = {1, ..., NP } then q return Φg ed break pt Algorithm 5 Particle swarm optimization update for particle q. best function PSOUpdate(vq , Φq , Φbest q , Φlbest ) ce Let vqz , vqy ∈ vq represent the velocities of the downstream variables. Let Φzq , Φyq ∈ Φq represent the values of the downstream variables. ⊲ Vector of random uniform numbers of size |z|. Ac r1 , r2 ← U [0, 1] r3 , r4 ← U [0, 1] ⊲ Vector of random uniform numbers of size |y|.

best vqz ← c1 · vqz + c2 · r1 · zbest − z + c · r · z − z q 3 2 q q lbest

y y best best vq ← c1 · vq + c2 · r3 · yq − yq + c3 · r4 · ylbest − yq zq ← zq + vqz ⊲ Update zq ∈ Φq . yq ← yq + vqy ⊲ Update yq ∈ Φq . return Φq 40 Page 41 of 53 t us cr ip an Algorithm 6 Differential evolution (“DE/rand/1/bin”) for agent q. best best best function DECrossover(Φbest q , Φa , Φb , Φc ) Φq ← Φbest q ⊲ Copy the best solution back. r ← U [0, 1] if r ≤ CR then ← a,best zε,t +F ·  b,best zε,t ed q zε,t M q for all zε,t ∈ zq do ⊲ Random uniform number. ⊲ Mutate the destination policies.  ⊲ Modify Φq . − c,best zε,t q for all yi,j,t,s ∈ yq do pt r ← U [0, 1] ⊲ Random uniform number. if r ≤ CR then ce q yi,j,t,s ⊲ Mutate the processing streams.   a,best b,best c,best ← yi,j,t,s + F · yi,j,t,s − yi,j,t,s ⊲ Modify Φq Ac return Φq 41 Page 42 of 53 References us cr ip t [1] G. Whittle, Global asset optimisation, in: R. Dimitrakopoulos (Ed.), Orebody modelling and strategic mine planning: Uncertainty and risk management models, AusIMM Spectrum Series, 331–336, 2007. [2] J. Whittle, The global optimiser works - what next?, in: R. Dimitrakopoulos (Ed.), Advances in Orebody Modelling and Strategic Mine an Planning I, AusIMM Spectrum Series, 3–5, 2009. [3] T. Johnson, Optimum open pit mine production scheduling, Ph.D. the- M sis, University of California: Berkeley, Berkeley, 1968. [4] K. Dagdelen, Optimum multi period open pit mine production schedul- ed ing, Ph.D. thesis, Colorado School of Mines, Golden, 1985. [5] M. E. Gershon, Optimal mine production scheduling: evaluation of large scale mathematical programming approaches, Int. J. Mining Eng. 1 (4) pt (1983) 315–329, doi:10.1007/BF00881548. ce [6] P. H. L. Monkhouse, G. Yeates, Beyond naive optimisation, in: R. Dimitrakopoulos (Ed.), Orebody modelling and strategic mine planning: Un- Ac certainty and risk management models, AusIMM Spectrum Series, 3–8, 2007. [7] B. S. Pimentel, G. R. Mateus, F. A. Almeida, Mathematical mdels for optimizing the global mining supply chain, in: B. Nag (Ed.), Intelligent Systems in Operations: Methods, Models and Applications in the Supply Chain, IGI Global, 133–163, doi:10.4018/978-1-61520-605-6.ch008, 2010. 42 Page 43 of 53 [8] S. Ramazan, R. Dimitrakopoulos, Production scheduling with uncertain us cr ip 14 (2) (2013) 361–380, doi:10.1007/s11081-012-9186-2. t supply: a new solution to the open pit mining problem, Optim. Eng. [9] J. Birge, F. Louveaux, Introduction to Stochastic Programming, Springer series in operations research and financial engineering, Springer, New York, 2 edn., 2011. an [10] A. Leite, R. Dimitrakopoulos, Stochastic optimisation model for open pit mine planning: application and risk analysis at copper deposit, Min. M Technol. 116 (3) (2007) 109–118, doi:10.1179/174328607X228848. [11] A. Jewbali, Modelling geological uncertainty for stochastic short-term production scheduling in open pit metal mines, Ph.D. thesis, University ed of Queensland, Brisbane, 2006. [12] J. Benndorf, R. Dimitrakopoulos, Stochastic long-term production pt scheduling of iron ore deposits: integrating joint multi-element ge- ce ological uncertainty, J. of Min. Sci. 49 (1) (2013) 68–81, doi: 10.1134/S1062739149010097. Ac [13] F. R. Albor Consuegra, R. Dimitrakopoulos, Algorithmic approach to pushback design based on stochastic programming: method, application and comparisons, Min. Technol. 119 (2) (2010) 88–101, doi: 10.1179/037178410X12780655704761. [14] A. Lamghari, R. Dimitrakopoulos, A diversified Tabu search approach for the open-pit mine production scheduling problem with 43 Page 44 of 53 metal uncertainty, Eur. J. Oper. Res. 222 (3) (2012) 642 – 652, doi: us cr ip t 10.1016/j.ejor.2012.05.029. [15] K. F. Lane, The Economic Definition of Ore: Cut-off Grades in Theory and Practice, Mining Journal Books Limited, London, 1988. [16] J. Rendu, An Introduction to Cut-off Grade Estimation, Society for Mining, Metallurgy, and Exploration, Englewood, 2 edn., 2011. an [17] A. Gleixner, Solving large-scale open pit mining production scheduling problems by integer programming, Master’s thesis, Technische Univer- M sität Berlin, Berlin, 2008. [18] N. Boland, I. Dumitrescu, G. Froyland, A. M. Gleixner, LP-based disag- ed gregation approaches to solving the open pit mining production scheduling problem with block processing selectivity, Comput. Oper. Res. 36 (4) pt (2009) 1064–1089, doi:10.1016/j.cor.2007.12.006. [19] D. Bienstock, M. Zuckerberg, Solving LP relaxations of large-scale prece- ce dence constrained problems, in: F. Eisenbrand, F. B. Shepherd (Eds.), IPCO’10 Proceedings of the 14th international conference on integer Ac programming and combinatorial optimization, Lecture Notes in Computer Science, Springer-Verlag, 1–14, doi:10.1007/978-3-642-13036-6 1, 2010. [20] W. B. Lambert, A. M. Newman, Tailored Lagrangian relaxation for the open pit block sequencing problem, Ann. of Oper. Res. 202 (1) (2013) 1–20, doi:10.1007/s10479-012-1287-y. 44 Page 45 of 53 [21] N. Boland, I. Dumitrescu, G. Froyland, A multistage stochas- scheduling with approach uncertain to open geology, pit mine Optimization production t programming Online URL us cr ip tic http://www.optimization-online.org/DB FILE/2008/10/2123.pdf. [22] M. Menabde, G. Froyland, P. Stone, G. Yeates, Mining schedule optimization for conditionally simulated orebodies, in: R. Dimitrakopoulos (Ed.), Orebody modelling and strategic mine planning: Uncertainty and an risk management models, AusIMM Spectrum Series, 379–384, 2007. [23] M. Kumral, Incorporating geo-metallurgical information into mine pro- 10.1057/jors.2009.174. M duction scheduling, J. Oper. Res. Soc. 62 (1) (2011) 60–68, doi: ed [24] M. Kumral, Optimizing ore–waste discrimination and block sequencing through simulated annealing, Appl. Soft Comput. 13 (8) (2013) 3737 – pt 3744, doi:10.1016/j.asoc.2013.03.005. ce [25] A. Bley, A. Gleixner, T. Koch, S. Vigerske, Comparing MIQCP solvers to a specialised algorithm for mine production Scheduling, in: H. G. Bock, X. P. Hoang, R. Rannacher, J. P. Schlöder (Eds.), Modeling, Ac Simulation and Optimization of Complex Processes, Springer, 25–39, doi:10.1007/978-3-642-25707-0 3, 2012. [26] S. Hoerger, J. Bachmann, K. Criss, E. Shortridge, Long term mine and process scheduling at Newmont’s Nevada operations, in: K. Dagdelen (Ed.), Proceedings of 28th APCOM Symposium, Society for Mining, Metallurgy, and Exploration, 739–748, 1999. 45 Page 46 of 53 [27] S. Hoerger, L. Hoffmann, F. Seymour, Mine planning at Newmont’s us cr ip t Nevada operations, Min. Eng.-Littleton 51 (10) (1999) 26–30. [28] E. K. Chanda, Network linear programming optimization of an integrated mining and metallurgical complex, in: R. Dimitrakopoulos (Ed.), Orebody modelling and strategic mine planning: Uncertainty and risk management models, AusIMM Spectrum Series, 149–155, 2007. an [29] P. Stone, G. Froyland, M. Menabde, B. Law, R. Pasyar, P. H. L. Monkhouse, Blasor – blended iron ore mine planning optimization at Yandi, Western Australia, in: R. Dimitrakopoulos (Ed.), Orebody mod- M elling and strategic mine planning: Uncertainty and risk management models, AusIMM Spectrum Series, 133–136, 2007. ed [30] C. Wharton, The use of extracted blending optimization for improved profitability, in: R. Dimitrakopoulos (Ed.), Orebody modelling and pt strategic mine planning: Uncertainty and risk management models, ce AusIMM Spectrum Series, 293–300, 2007. [31] T. Sandeman, C. Fricke, P. Bodon, C. Stanford, Integrating optimization and simulation - a comparison of two case studies in mine planning, in: Ac B. Johansson, S. Jain, J. Montoya-Torres, J. Hugan, E. Yucesan (Eds.), Proceedings of the 2010 Winter Simulation Conference (WSC), ISSN 0891-7736, 1898–1910, doi:10.1109/WSC.2010.5678882, 2010. [32] R. Epstein, M. Goic, A. Weintraub, J. Catalán, P. Santibáñez, R. Urrutia, R. Cancino, S. Gaete, A. Aguayo, F. Caro, Optimizing long-term 46 Page 47 of 53 production plans in underground and open-pit copper mines, Oper. Res. us cr ip t 60 (1) (2012) 4–17, doi:10.1287/opre.1110.1003. [33] G. Singh, R. Garcı́a-Flores, A. Ernst, P. Welgama, M. Zhang, K. Munday, Medium-term rail scheduling for an iron ore mining company, Interfaces 44 (2) (2013) 222–240, doi:10.1287/inte.1120.0669. [34] M. Godoy, The effective management of geological risk in long-term an production scheduling of open pit mines, Ph.D. thesis, University of Queensland, Brisbane, 2004. M [35] J. Ferland, J. Amaya, M. Djuimo, Application of a particle swarm algorithm to the capacitated open pit mining, in: S. Mukhopadhyay, G. Sen Gupta (Eds.), Autonomous Robots and Agents, Studies in Com- ed putational Intelligence, Springer-Verlag, 127–134, 2007. [36] J. Sattarvand, Long-term open-pit planning by ant colony optimization, pt Ph.D. thesis, RWTH Aachen University, Aachen, 2009. ce [37] R. Goodfellow, R. Dimitrakopoulos, Algorithmic integration of geological uncertainty in pushback designs for complex multipro- Ac cess open pit mines, Min. Technol. 122 (2) (2013) 67–77, doi: 10.1179/147490013X13639459465736. [38] L. Montiel, R. Dimitrakopoulos, Stochastic mine production scheduling with multiple processes: application at Escondida Norte, Chile, J. Min. Sci. 49 (4) (2013) 583–597, doi:10.1134/S1062739149040096. [39] A. Khan, C. Niemann-Delius, Production scheduling of open pit mines 47 Page 48 of 53 using particle swarm optimization algorithm, Adv. in Oper. Res. 2014 us cr ip t (2014) 1–9, doi:10.1155/2014/208502. [40] A. Lamghari, R. Dimitrakopoulos, J. A. Ferland, A variable neighbourhood descent algorithm for the open-pit mine production scheduling problem with metal uncertainty, J. Oper. Res. Soc. 65 (9) (2014) 1305– 1314, doi:10.1057/jors.2013.81. an [41] M. S. Shishvan, J. Sattarvand, Long term production planning of open pit mines by ant colony optimization, Eur. J. Oper. Res. 240 (3) (2015) M 825–836, doi:10.1016/j.ejor.2014.07.040. [42] R. Khalokakaie, P. A. Dowd, P. A. Dowd, R. J. Fowell, R. J. Fowell, Lerchs - Grossmann algorithm with variable slope angles, Min. Technol. ed 109 (2) (2000) 77–85. [43] S. Lloyd, Least squares quantization in PCM, IEEE T. Inform. Theory pt 28 (2) (1982) 129–137, doi:10.1109/TIT.1982.1056489. ce [44] D. Arthur, S. Vassilvitskii, K-means++: The advantages of careful seeding, in: G. Habow (Ed.), Proceedings of the Eighteenth Annual ACM- Ac SIAM Symposium on Discrete Algorithms, SODA ’07, Society for Industrial and Applied Mathematics, 1027–1035, 2007. [45] J. Sattarvand, C. Niemann-Delius, Past, present and future of metaheuristic optimization methods in long-term production planning of open pits, BHM 158 (4) (2013) 146–154, ISSN 0005-8912, doi: 10.1007/s00501-013-0127-y. 48 Page 49 of 53 [46] B. Denby, D. Schofield, Open-pit design and scheduling by use of genetic t algorithms, Trans. Inst. Mining Met. Sect. A: Mining Ind. 103 (1994) us cr ip 21–26. [47] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, E. Teller, Equation of state calculations by fast computing machines, J. Chem. Phys. 21 (6) (1953) 1087–1092, doi:10.1063/1.1699114. simulated annealing, an [48] S. Kirkpatrick, C. D. Gelatt, M. P. Vecchi, Optimization by Science 220 (4598) (1983) 671–680, M 10.1126/science.220.4598.671. doi: [49] S. Geman, D. Geman, Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images, IEEE T. on Pattern. Anal. PAMI- ed 6 (6) (1984) 721–741, doi:10.1109/TPAMI.1984.4767596. [50] I. O. Bohachevsky, M. E. Johnson, M. L. Stein, Generalized simulated pt annealing for function optimization, Technometrics 28 (3) (1986) 209– ce 217, doi:10.2307/1269076. [51] A. Corana, M. Marchesi, C. Martini, S. Ridella, Minimizing multi- Ac modal functions of continuous variables with the simulated annealing algorithm, ACM Trans. Math. Softw. 13 (3) (1987) 262–280, doi: 10.1145/29380.29864. [52] J. Kennedy, R. Eberhart, Particle swarm optimization, in: IEEE International Conference on Neural Networks, IEEE Press, 1942–1948, doi:10.1109/ICNN.1995.488968, 1995. 49 Page 50 of 53 [53] R. Poli, J. Kennedy, T. Blackwell, Particle swarm optimization, Swarm us cr ip t Intell. 1 (1) (2007) 33–57, doi:10.1007/s11721-007-0002-0. [54] R. Storn, K. Price, Differential Evolution – a simple and efficient heuristic for global optimization over continuous spaces, J. Global Optim. 11 (4) (1997) 341–359, doi:10.1023/A:1008202821328. [55] K. Price, R. M. Storn, J. A. Lampinen, Differential Evolution: A Springer-Verlag, Secaucus, 2005. an Practical Approach to Global Optimization, Natural Computing Series, M [56] J. Kennedy, R. Mendes, Population structure and particle swarm performance, in: D. B. Fogel, M. A. El-Sharkawi, X. Yao, G. Greenwood, H. Iba, P. Marrow, M. Shackleton (Eds.), Proceedings of the 2002 ed Congress on Evolutionary Computation, vol. 2, IEEE Press, 1671–1676, doi:10.1109/CEC.2002.1004493, 2002. pt [57] J. Brest, S. Greiner, B. Boskovic, M. Mernik, V. Zumer, Self-adapting ce control parameters in differential evolution: a comparative study on numerical benchmark problems, IEEE Trans. on Evol. Comp. 10 (6) Ac (2006) 646–657, doi:10.1109/TEVC.2006.872133. [58] D. Zaharie, Influence of crossover on the behavior of differential evolution algorithms, Appl. Soft Comput. 9 (3) (2009) 1126–1138, doi: 10.1016/j.asoc.2009.02.012. [59] P. Goovaerts, Geostatistics for Natural Resources Evaluation, Applied geostatistics series, Oxford University Press, Oxford, 1997. 50 Page 51 of 53 us cr ip t Biographies Ryan Goodfellow is a post-doctoral research fellow at McGill Univer- an sity’s COSMO Laboratory. He received a Bachelor’s of Mining Engineering pt ed McGill University in 2014. M in 2009 from McGill and a Ph.D. in Mining and Materials Engineering from ce Roussos Dimitrakopoulos is a Professor at McGill University, Montreal, and holds the Canada Research Chair (Tier I) in Sustainable Mineral Ac Resource Development and Optimization under Uncertainty. He leads the COSMO Research Consortium of major global mining companies, namely, AngloGold Ashanti, Barrick Gold, BHP Billiton, DeBeers, Newmont Mining, and Vale. 51 Page 52 of 53 Geological Simulations Ac ce pt ed M an us cr ip t *Graphical abstract (for review) Clustering for Destination Policies Optimization Material Simulations Cross-section of production schedule All Simulations - Material #1 Simulation #1 Analysis of Stochastic Design Destination Policy Downstream Optimization with PSO or DE Robust destination policies Sulfide Waste Dump Adaptive processing streams Sulfide Heap Leach Stockpile Sulfide Mill Simulation #2 Global best solution Simulation #1 Mining Complex Model Stockpile Sulfide Mill Simulation #2 Copper-Gold Mine Global Optimization with Simulated Annealing Robust production schedule Sulfide Heap Leach Sulfide Waste Dump Transition Heap Leach Oxide Heap Leach Robust destination policies Adaptive processing streams Cumulative Discounted Cash Flows (Compared as % to Commercial) Grade Simulations Cumulative Discounted Cash Flows 140% 120% 100% 80% 60% Stochastic P-10/90 Stochastic P-50 Deterministic-Equivalent Commercial (Deterministic) 40% Page 53 of 53 20% 0% 0 5 10 15 Period Oxide Waste 20 25 25