Implementing a parametric maximum flow ... - Semantic Scholar

6 downloads 142278 Views 293KB Size Report
time-dependent discounted values of mining blocks when designing optimal production phases and ultimate pit .... application of the proposed stochastic optimization app- roach at a ... interest, λ, usually metal price, identifies successive pit.
Journal of the Operational Research Society (2012), 1–13

© 2012 Operational Research Society Ltd. All rights reserved. 0160-5682/12 www.palgrave-journals.com/jors/

Implementing a parametric maximum flow algorithm for optimal open pit mine design under uncertain supply and demand MWA Asad and R Dimitrakopoulos COSMO Stochastic Mine Planning Laboratory, McGill University, Canada Conventional open pit mine optimization models for designing mining phases and ultimate pit limit do not consider expected variations and uncertainty in metal content available in a mineral deposit (supply) and commodity prices (market demand). Unlike the conventional approach, a stochastic framework relies on multiple realizations of the input data so as to account for uncertainty in metal content and financial parameters, reflecting potential supply and demand. This paper presents a new method that jointly considers uncertainty in metal content and commodity prices, and incorporates time-dependent discounted values of mining blocks when designing optimal production phases and ultimate pit limit, while honouring production capacity constraints. The structure of a graph representing the stochastic framework is proposed, and it is solved with a parametric maximum flow algorithm. Lagragnian relaxation and the subgradient method are integrated in the proposed approach to facilitate producing practical designs. An application at a copper deposit in Canada demonstrates the practical aspects of the approach and quality of solutions over conventional methods, as well as the effectiveness of the proposed stochastic approach in solving mine planning and design problems. Journal of the Operational Research Society advance online publication, 25 April 2012 doi:10.1057/jors.2012.26 Keywords: open pit mine optimization; maximum flow algorithm; Lagrangian relaxation; subgradient method

1. Introduction An open pit mine refers to a large-scale surface mining operation dedicated to the extraction of raw materials, including metallic ores. The existence of economic ore reserves, that is, the availability of the sufficient amount of valuable raw material, is a prerequisite to commencing a mining operation. Exploratory drilling and subsequent geological interpretation establish the size, shape, depth, and geology of the economic ore reserves. Spatial or geostatitical modelling of this data is based on the division of the volume of the reserve into a number of fixed size mining units, while generating the quantitative data used subsequently in optimization, such as the tonnage and available metal content in each mining unit. A mining unit is traditionally referred to as a mining block, whose value is derived from an estimated grade or metal content and economic parameters such as the commodity price and mining and processing costs. A mining block having 

Correspondence: MWA Asad, COSMO Stochastic Mine Planning Laboratory, Department of Mining and Materials Engineering, McGill University, FDA Building, 3450 University Street, Montreal, Quebec, Canada. E-mail: [email protected]

positive value becomes an ore block destined for the processing stream, while a block with negative value is identified as a waste block destined for the waste dump. A typical three-dimensional economic orebody model consists of millions of blocks with their respective values and becomes the basic input to forthcoming production planning and optimization (Hustrulid and Kuchta, 2006; Ramazan, 2007). The objective of long-term open pit mine planning is to maximize net present value over the life of the mining operation. This is accomplished with an optimal production schedule identifying blocks to be mined and processed over a number of periods subject to the mining and processing capacity constraints. However, the size of this mixed integer programming-based production scheduling problem limits the solution strategies to heuristics that may lead to possibly suboptimal solutions (Denby et al, 1998; Akaike and Dagdelen, 1999; Caccetta and Hill, 2003). Alternatively, it is an established practice to develop separate and sequential models to arrive at a complete long-term open pit mine plan (Lerchs and Grossmann, 1965; Tachefine and Soumis, 1997; Whittle, 1988; Hochbaum and Chen, 2000; Hustrulid and Kuchta, 2006). Thus, long-term planning constitutes a two-tier

2

Journal of the Operational Research Society

process which identifies: (i) the optimal production phase design that refers to the phase-wise sequence of extraction reaching to the ultimate pit limit or overall extent of extraction (that is, mining beyond ultimate pit limit would be uneconomical), as illustrated in Figure 1; and (ii) the optimal production schedule of the blocks within each phase (that is, a subset of the blocks included in an orebody model), such that the discounted cash flows are maximized, subject to the production capacities of a variety of processes comprising an open pit mining operation. Hence, a production phase within the optimal ultimate pit limit constitutes a manageable unit of reserves that may be exhausted over a number of time periods, using dedicated human and physical resources (Hustrulid and Kuchta, 2006). The work presented herein addresses the first segment of the long-term plan, that is, the development of an optimal phase design and an optimal ultimate pit limit. Optimization of long-term open pit mine plans may be accomplished through conventional (deterministic) or stochastic approaches. The conventional approaches do not consider input uncertainty and are based on an estimated orebody model and constant commodity prices and operating costs. The underlying assumption in conventional approaches is that the inputs derived from geological data (metal content) and economic data (commodity price, mining, and processing costs) are known and constant. As a result, the economic values of mining blocks calculated are also treated as constant and are processed as if they are the actual economic values of the mining blocks. The assumption of constant economic values of mining blocks is not necessarily realistic, particularly given the inherent uncertainty in the input data, which is well-known for both the metal content of mining blocks (Dowd, 1994) as well as commodity price (Chen, 2010). The assumption of constant block values may result in technically unrealistic mine design analysis and significant economic losses and a possible failure of a project requiring substantial capital investment (Dimitrakopoulos, 2011; Godoy and Dimitrakopoulos, 2011). Although metal supply uncertainty has been considered in the past through stochastic integer programming formulations (Ramazan and Dimitrakopoulos, 2012) the existing formulations are

Phase A

Ultimate Pit Limit

relatively restricted to a few tens of thousands of mining blocks. Optimization of mine design jointly integrating metal and commodity price uncertainty has not been addressed in the past and is contributed herein. To further elucidate as appropriate, it is noted that the uncertainty in metal content of a particular block relates to the uncertainty in supply of raw material to the upstream ore processing facilities. Similarly, the uncertainty in commodity price corresponds to the demand of that particular commodity because apart from other factors, the uncertainty in metal demand is the major factor contributing to the metal price fluctuations (Bower et al, 2007; Chen, 2010). As such, a unique aspect of the proposed approach is that it not only accounts for the uncertainty over the supply of raw materials, but also looks into the consequences of uncertainty on the demand side of the open pit mine planning problem. The uncertainty of these inputs suggests that there is a possibility of a block being wrongly identified as an ore block when it may actually be a waste block at the time of extraction, and vice versa. The proposed stochastic approach utilizes a set of simulated orebody realizations or equally probable models as an input for developing optimal long-term open pit mine plans. The stochastic optimization model addresses the ultimate pit limit and phase design problem by representing it as a graph (Johnson, 1968; Picard, 1976). More specifically, the paper (i) proposes the structure of the graph in a stochastic framework that incorporates uncertainty in metal content and commodity prices, that is, uncertainty in both potential supply and the demand; (ii) establishes that the structure of the graph in a stochastic framework is different from the graph representing a conventional or deterministic framework; (iii) applies a widely accepted maximum flow algorithm to this stochastic framework; (iv) incorporates production capacity constraints and Lagrangian relaxation of these constraints to exploit the classical structure of the maximum closure problem, while utilizing the subgradient method for systematic selection of the Lagrangian parameters to develop successive phases to the ultimate pit limit; and, finally, (v) demonstrates the application of the proposed stochastic optimization approach at a copper deposit in Canada, which suggests an improvement in optimal ultimate pit limit and phase design as compared to the existing conventional solutions. In the following sections the paper presents: (i) a review of previous research; (ii) the structure of the graph for a stochastic framework utilizing discounted values of mining blocks; (iii) an implementation of the parametric maximum flow (minimum cut) algorithm; (iv) application of the method and comparisons; and then, (v) conclusions follow.

Phase B

1.1. Previous research Phase C

Figure 1 Production phases within the ultimate pit limit.

The ultimate pit limit and phase design problem may be represented as a graph defined on a directed graph

MWA Asad and R Dimitrakopoulos—Implementing a parametric maximum flow algorithm

G ¼ (V, A), in which each block of value ci becomes a node in V and a directed arc, in A, extends from node i to node j where block j overlies block i. This sets an extraction precedence of block j over block i. The problem lies in finding a set of nodes in V0 DV, which includes maximum P value nodes and all successor nodes such that iAV0 ci maximized. This generates a maximum closure on the graph (Hochbaum and Chen, 2000; Chandran and Hochbaum, 2009). The most common solution strategy in the mining industry is the graph theory-based (Lerchs and Grossmann, 1965) algorithm, or LG algorithm. The LG algorithm solves optimally only the ultimate pit limit problem, and it is implemented in several commercial mine planning software (Whittle, 1988). The LG algorithm represents the established conventional pit planning models. A separate or parametric implementation of the LG algorithm solves for phase design in a conventional mine planning framework, in which increase in a parameter of interest, l, usually metal price, identifies successive pit shells that are used to develop phases within the ultimate pit limit (Whittle, 1988). This separates the problem into two parts and requires repeat computation and evaluation of phase design with a range of l values, selected through trial and error, resulting in a series of nested pits, a process referred to as parameterization. Several variations of pit parameterization algorithms include Francois-Bongarcon and Guibal (1984), Coleou (1989), and Wang and Sevim (1993). A common drawback of pit parameterization methodologies is the so-termed ‘gap’ (Meagher et al, 2009). The Gap refers to operationally unmanageable phase design within pit limit primarily due to the drastic difference in size and quantity of material to be mined between successive phases. In some cases, a phase may not even be connected to the previous phase, thereby requiring frequent equipment relocations that render the phase design impractical. As an alternative to the LG algorithm, Johnson (1968) established a relationship between the ultimate pit limit problem and the maximum flow algorithm. A mathematical formulation of the problem is given as (Tachefine and Soumis, 1997; Hochbaum and Chen, 2000; Hochbaum, 2003) X ci xi ð1Þ Maximize: Z ¼ i2V

3

ci(l) ¼ hi þ lqi, where hi is a parameter such as mining or processing cost and qi is parameter such as metal price. A particular value of l may result in an increase or a decrease in metal price, which may result in an ultimate increase or decrease in the block value ci. The totally unimodular constraint matrix in the mathematical formulation reflects that this problem is the dual of the maximum flow problem, that is, a minimum cut problem (Hochbaum and Chen, 2000; Hochbaum, 2004). Picard (1976) demonstrates that the minimum cut algorithm applied to a related graph including a source node, S, and a sink node, R, solves the maximum closure problem, that is, the minimum cut algorithm generates the optimal ultimate pit limit. Numerous algorithms such as Goldberg and Tarjan (1988), Ahuja and Orlin (1989), Gallo et al (1989), Faaland et al (1990), Tachefine and Soumis (1997), and Hochbaum (2001, 2008) are available for solving maximum flow and parametric maximum flow problems. Most studies present comparative analyses between algorithms and evaluate algorithmic performance through an application on the challenging open pit mining problem. Hochbaum and Chen (2000) conclude that the push-relabel algorithm (one of the efficient algorithms to compute maximum flow, see Goldberg and Tarjan (1988)) performs better than the LG algorithm in solving for ultimate pit limit and phase design. The study also demonstrates the lack of robustness of the LG algorithm due to its sensitivity to the distribution of ore and waste blocks in the input data. Hochbaum (2001) provides detailed analysis of the LG algorithm. The LG algorithm is more complex than the push-relabel algorithm in cases where an increase occurs in the number of ore blocks, the mine size, or the number of successor blocks included in the set xi related to variants of the precedence constraint pattern. Exploiting the structure of the maximum closure problem proposed in Picard (1976), Tachefine and Soumis (1997) include resource or open pit production capacity constraints and relax the problem through Lagrangian relaxation. If K is the set of production capacity constraints, aki is the consumption of capacity (resource) k at node i, that is, the quantity of material in block i and bk is the availability of resource k, then the Equations (1)–(3) may be augmented with the constraint X aki xi pbk ; k 2 K: ð4Þ i2V

Subject to: xi  xi 0 p0; i 0 2 xi ; i 2 V

ð2Þ

xi 2 f0; 1g; i 2 V:

ð3Þ

The relaxed formulation then becomes Maximize: ZðlÞ ¼

XX

½½cik  lk aki xi  þ

i2V k2K

where xi is a binary variable equal to 1 if node i is in the closure, and 0 otherwise, and xi is the set of successors of node i. In developing phase design through parameterization, the block value may be modified using the equation

X

lk bk ð5Þ

k2K

Subject to: xi  xi 0 p0; i 0 2 xi ; i 2 V

ð6Þ

xi 2 f0; 1g; i 2 V:

ð7Þ

4

Journal of the Operational Research Society

where lkX0 denotes the Lagrangian multiplier associated with the production capacity constraint in Equation (4). The production capacity constraint is multiplied by lk and introduced in the objective function of Equation (1) as a penalty, which produces the objective function of Equation (5). Equation (4) shows that only the upper bound on the production capacity constraint is considered, here. This allows using the ultimate pit limit as a guide in the optimization model. Note that, in case the lower bound is present, some of the blocks included in the solution could be outside the optimal ultimate pit limit. Tachefine and Soumis (1997) solve the relaxed problem through the subgradient, the bundle (Kelley’s cutting plane) method, and heuristic approaches implemented through a pushrelabel algorithm. However, these studies do not incorporate uncertainty, assuming that block values are constant with negligible effect on ultimate pit limit and phase design solutions. Meagher et al (2009) addressed the ultimate pit limit and phase design problem by proposing a theoretical background utilizing a parametric minimum cut algorithm in a stochastic framework that incorporates multiple simulated orebody models and economic parameter uncertainty. Unlike the conventional framework, which ignores the available information on quantifying risk (such as the probability of being an ore block or a waste block) by essentially averaging multiple simulated orebody models, this approach blends ore block with high metal content and low risk inside the pit with blocks carrying high risk. For example, if a mining block is worth 3000, 3000, and þ 9000 in three different orebody realizations, the average value of this block is þ 1000. Another block may be worth þ 1000 in all three realizations. The average value is again þ 1000. The conventional framework fails to distinguish between high or low risk blocks. The stochastic approach will include in the pit the second block because its probability of being an ore block is equal to 1, but will exclude the first block because of the higher risk. Meagher et al (2009) create single directed graphs incorporating multiple simulated realizations of the orebody and expected economic parameters. Therefore, the modified graph honouring this stochastic framework is different from the one created for conventional framework in that the latter utilizes estimated (expected) block values. Previous implementations of the maximum flow algorithm do not jointly account for geological and market uncertainties. This paper proposes an extended stochastic approach to solve for ultimate pit limit and phase design, which utilizes both metal price and grade simulations in calculating an orebody model carrying a set of simulated realizations of the discounted block values (Dimitrakopoulos and Abdel Sabour, 2007). We propose an extension of the concepts presented by Meagher et al (2009). We honour production capacity constraints and incorporate Lagrangian relaxation of these constraints to

exploit the classical structure of the maximum flow algorithm. Unlike experimental procedures implemented in previous studies (Whittle, 1999), the subgradient method systematically selects the Lagrangian parameter, l, over a number of iterations, thereby decreasing ore block values to generate ultimate pit limit and successive phases with the ultimate pit limit. We address the gap problem by introducing a modification in the subgradient method, such that the successive phases contain uniform amounts of material and are thus operationally manageable. We also demonstrate the value of proposed stochastic approach in an application that promises an improvement in optimal ultimate pit limit and phase design.

2. Phase design and ultimate pit limit Developing a stochastic framework to solve for ultimate pit limit and phase design is a complex process. It includes determining multiple realizations of timedependent discounted block values, creation of a single graph for multiple realizations, and implementation of the parametric maximum flow (minimum cut) algorithm. A detailed outline of this sequence of steps is presented in the following subsections.

2.1. Calculating economic block values The procedure for calculating the discounted values of mining blocks incorporates the uncertainty in commodity price and metal content jointly. Commodity price modelling initially captures the market or demand uncertainty by employing a stochastic model describing the evolution of commodity prices over time. Then, the value of each mining block is evaluated by considering simulated grade realizations and mining, processing and refining/smelting costs in combination with commodity prices. The following model parameters are defined: i ¼ block indicator, where i ¼ 1, . . . , I; t ¼ period (year) indicator, where t ¼ 1, . . . , T: g ¼ grade realization (or simulation) indicator, where g ¼ 1, . . . , G; p ¼ price realizations indicator, where p ¼ 1, . . . , P; ubVpgi ¼ undiscounted value of block i in grade realization g and price realization p; bVgi ¼ average undiscounted value of block i in grade realization g; dbVgti ¼ discounted value of block i in period t and grade realization g; Qi ¼ available quantity (tons) of material in block i; ggi ¼ grade of block i in realization g; P ¼ current or initial metal price ($/lb. or oz. of metal); dP ¼ change in metal price over an interval between two periods dt;

MWA Asad and R Dimitrakopoulos—Implementing a parametric maximum flow algorithm

Z ¼ speed at which the metal price returns back to the equilibrium level after a shock in the market; P ¼ long-term equilibrium metal price; s ¼ standard deviation of the metal price; o ¼ N(0, 1) random number; y ¼ metallurgical recovery (%); m ¼ mining cost ($/ton of material); p ¼ processing cost ($/ton of ore); r ¼ refining cost ($/lb. or oz. of metal); d ¼ discount rate (%). The steps for determining discounted block values are as follows: 1. For a time interval dt, simulate price realizations using the mean reverting price model (Schwartz, 1997): pffiffiffiffi dP ¼ ZðP  PÞdt þ so dt

ð8Þ

Equation (8) consists of two terms, the deterministic and the stochastic terms, respectively. Given the price volatility s, the stochastic term perturbs the commodity price away from the equilibrium level. The deterministic term acts to restore the commodity price to the equilibrium level (Schwartz, 1997). 2. For each price realization, p, and grade realization, g, calculate the undiscounted value of Block i using the relation ubVpgi ¼ ðP þ dP  rÞQi ggi y  mQi  pQi :

ð9Þ

3. Calculate the average undiscounted value of Block i using the relation PP bVgi ¼ 4. Adjust the average that 8 bVgi > > > < bVgi ¼ > > > mQi :

ubVpgi : P

p¼1

ð10Þ

undiscounted value of Block i such if bVgi 40; ie, i is an ore block in grade realization g; if bVgi p0; ie, i is waste block in grade realization g:

The cutoff grade is the established economic criterion for identifying a particular block i as an ore block or a waste block. The cutoff grade is derived from the price of the metal, metallurgical recovery, mining, processing, and refining costs. Hence, a block of grade higher than the cutoff grade becomes an ore block with positive value, ensuring that it pays off the cost of mining, processing, and refining (Equation (9)). A block with grade lower than the cutoff grade becomes a waste block with a negative value that pays the cost of mining and transport to the waste dump (Lane, 1964, 1988;

5

King, 2011). The cutoff grade is dynamic in the sense that it declines during the life of the operation due to the exhaustion of the reserves (Asad, 2007; Bascetin and Nieto, 2007). However, it is the established practice to assume a fixed cutoff grade schedule at least for determining the ultimate pit limit, and the work herein also assumes a fixed cutoff grade (Lerchs and Grossmann, 1965; Whittle, 1988, 1999; Hochbaum and Chen, 2000; Hustrulid and Kuchta, 2006; Hochbaum, 2008). In addition, a stochastic model that considers metal content or price uncertainty does not guarantee that the ore blocks or waste blocks are correctly identified as a result of assuming a dynamic cutoff grade. The available metal content in the Block i then dictates the block as an ore block or a waste block. The model parameter Qi is related to the available metal content. If the available metal content generates value after paying the costs of mining, processing, and refining, then Qi equals the amount (tons) of ore. Otherwise, the available material in Block i is equal to the amount of waste. Uncertainty in metal content and commodity price allows a block to be an ore block in one grade realization and a waste block in another grade realization. 5. Calculate the discounted block value as dbVgti ¼

bVgi : ð1 þ dÞt

ð11Þ

The steps outlined above generate an orebody model consisting of G realizations of time-dependent, discounted mine block values. In a stochastic framework, the amount of the input data (number of block values) increased to ‘GTI ’ from ‘I ’ in a conventional framework.

2.2. The graph structure Having described the orebody model, the next step is to present the structure of the graph including multiple realizations of time-dependent discounted block values. For simplicity, we describe the process of generating a single graph in three steps. First, the relationship between the ultimate pit limit and phase design problem and the maximum closure on the graph is established, considering a single orebody realization. Second, the relationship is extended to multiple realizations with undiscounted block values. Third, the proposed method is demonstrated on multiple realizations with time-dependent discounted block values.

2.2.1. Single realization. The ultimate pit limit and phase design problem is equivalent to a maximum closure problem (Picard, 1976). A source, S, and a sink, R, are augmented to the graph G ¼ (V, A) to create the related graph G˜ ¼ (V˜, A˜) such that V˜ ¼ V [ {S, R}

6

Journal of the Operational Research Society

(Hochbaum and Chen, 2000). The related graph consists of sets of nodes V þ ¼ {iAV|ci40} and V  ¼ {iAV|cip0} representing ore blocks and waste blocks, respectively, and a set of arcs {A} [ {(S, v)|vAV þ } [ {(v, R)|vAV }. The capacity of all arcs in A is set to N. The capacity of each arc connected to the source or sink is |cv|, such that c(S, v) ¼ cv for vAV þ and c(v, R) ¼ cv for vAV . The source set of minimum cut separating S and R is the maximum closure in the related graph, since minimizing the cut capacity is equivalent to maximizing the sum of the node values in the source set of the cut. The structure of the related graph G˜ ¼ (V˜, A˜) is explained through an example in Figure 2. Figure 2 shows a related graph consisting of 15 blocks. The top-left block is labelled Block 1. The bottom-right block is Block 15. The economic values of the blocks are indicated inside each block. The related graph consists of the following elements: 1. Seventeen nodes including the source and sink nodes (S and R), seven ore block nodes, and eight waste block nodes. 2. Seven arcs from the source node to ore block Nodes 1, 2, 4, 7, 9, 12, and 13 of capacity 4, 8, 4, 14, 7, 6, and 6, respectively.

s

4 8 4 7

4

8

-4

4

-4

-6

14

-2

7

-3

-4

6

6

-9

-2

6

2 9

4

t

Figure 2 Structure of the related graph in a single realization.

3. Eight arcs from waste block Nodes 3, 5, 6, 8, 10, 11, 14, and 15 to the sink node of capacity (absolute values) 4, 4, 6, 2, 3, 4, 9, and 2, respectively. 4. Twenty six arcs of capacity N, representing block’s extraction precedence or slope requirements. For example, Blocks 2, 3, and 4 must be extracted before Block 8 is considered for extraction. In addition, Figure 2 outlines the ultimate pit limit obtained by the minimum cut algorithm. The value of this cut is equal to the sum of the flow capacities of the arcs within the cut (mined portion) originating at the source and terminating at the sink. The minimum cut algorithm for an open pit mining problem minimizes the number of waste blocks on the source node side of the cut and the number of ore blocks on the sink node side of the cut. Alternatively, it simultaneously maximizes the ore and minimizes waste inside the valid pit while respecting slope constraints (Meagher et al, 2009). The value of a cut is given by the sum of the capacities of waste blocks on the source side and the capacities of ore blocks on the sink side; that is, the value of a cut corresponds to the relative cumulative value of the waste and ore blocks in the source and sink side, respectively. In this example, the value of the cut is 16.

2.2.2. Multiple realizations—undiscounted block values. A stochastic framework that considers geological uncertainty consists of multiple realizations of undiscounted block values. The related graph concept can be extended so that the source node is connected to the ore blocks and all waste blocks are connected to the sink node. However, a particular block may be simulated as an ore block in one realization and a waste block in another realization. Maintaining a difference in arc capacity, the source node will be connected to this block in one realization, and will be connected as a waste block to the sink node in another realization. The slope constraints are maintained through infinite-capacity arcs in each realization. Figure 3 shows the arcs from the source node to ore blocks and from waste blocks to the sink node in a stochastic framework consisting of three realizations. Block 1 is an ore block in Realizations 1 and 2, but a waste block in Realization 3. A source to Block 1 arc is maintained in Realizations 1 and 2, with capacities of 4 and 2, respectively. A Block 1 to sink arc of capacity 10 is maintained in Realization 3. Once the related graph is constructed, the minimum cut algorithm is applied to solve for the ultimate pit limit. The approach keeps a given block in different realizations on the same side of the cut when generating a valid pit. Regardless of the number of realizations, the decision to mine a block is always binary, and a particular block is either inside or outside the pit. Hence, unlike the conventional approach, the stochastic

MWA Asad and R Dimitrakopoulos—Implementing a parametric maximum flow algorithm

s

s

6 2

4

5

-4

-6 14 -2

7

-3

-4

-9 -2

6

6

2

6

-2 -4 -4

-3

8

4

9

-3

-2

4

7

5

-1



-10 3

-6

6

-4

-4

3

2

5

6

-1

2

2

6

5

Realization 3

4

Realization 2

-4

Realization 1

8

12

17 -10

0



0

-12 6

-4 21

0 12

5 10

0

-13

11

0

25

0

9

17

6

6 4

7

6

-2 15

-12

0 11

-6 5

10

-7

0 7

t

0 4

10

-9 6

-3 3

Figure 3 Graph structure in multiple realizations with undiscounted block values. t

framework requires an additional constraint that ensures that a given block falls on the same side of the cut in all realizations (Meagher et al, 2009). This constraint may be satisfied by maintaining bidirectional infinite-capacity arcs among the same blocks in different realizations. Figure 3 indicates two bidirectional infinite-capacity arcs for Block 11 from Realization 1 to Realization 2 and from Realization 2 to Realization 3. A total of 30 bidirectional arcs are introduced in this example case. In a stochastic framework, the number of nodes and arcs is greater than the number of nodes and arcs in a conventional framework. An implementation of the maximum flow algorithm for a stochastic framework requires processing a significant amount of data, which introduces computational complexity as compared to the conventional framework. By maintaining bidirectional arcs across realizations, a block will be on the same side of the cut in different realizations. Therefore, bidirectional arcs allow merging of the nodes representing a particular block in different realizations into a single node. This newly merged (single) node eliminates the need for bidirectional infinitecapacity arcs and reduces the number of nodes and arcs, thus leading to reduction in computational complexity. Finally, a single arc is introduced from the source to this new node if the new node is an ore block or from the new node to the sink if the new node is a waste block. The capacity of the new arc is the sum of the capacities of the arcs from the source to the merged nodes or the sum of capacities of the arcs from the merged waste nodes to the sink over all realizations. Figure 4 depicts this concept. Three realizations as in Figure 3 have been merged and the bidirectional infinite-capacity arcs were eliminated to simplify the related graph.

Figure 4 Structure of the merged related graph in multiple realizations with undiscounted block values.

In the merged graph, an arc from S to Block 1 with a capacity of 6 is maintained. This capacity is the sum of the capacities of the arcs from S to Block 1 in Realizations 1 and 2 (Figure 3). Similarly, an arc from Block 1 to R with a capacity of 10 is maintained, and this capacity is equal to the capacity of the same block in Realization 3 (Figure 3). Thus, merging Block 1 eliminates one arc from S to this block and two bidirectional infinite-capacity arcs. In this example, an overall reduction of 54 arcs including 30 bidirectional infinite-capacity arcs is achieved. For an actual orebody model, merging the nodes of a particular block over a number of realizations results in a significant decrease in the size of the input matrix, which in turn simplifies the computation. The proposed structure of the related graph has the following features: 1. The number of nodes remains constant regardless of the number of realizations. 2. The capacities of arcs from source to ore blocks and from waste blocks to sink differ. K

K

If a block is an ore block in all realizations; then a single arc from the source to the node is maintained, and the capacity of this arc is the sum of the capacities of the arcs from the source to this block across all realizations. If a block is a waste block in all realizations; then a single arc is maintained from the node to the

8

Journal of the Operational Research Society

K

2.3. The parametric maximum flow (minimum cut) algorithm Extending Equations (1)–(4) to develop a related graph G˜, which considers discounted block values (as described in Section 2.2.3) in combination with pit production capacity

-4

4

-4

-6 14 -4

4

8

7

-3

-4

-9 -4

6

6

3-2

3

-3

3

-3

-4 11 -3

4

-2

-3

-6 -3

| - 9 + 6| = 3

6

2

4

4

| - 6 + 4| = 2

2

3

-2

2

-2

-3

9

-2

2

-1

-2

3

3

-4 -2

Period 3

2.2.3. Multiple realizations—discounted block values. A stochastic framework takes multiple realizations of time-dependent discounted mining block values as an input; as such, it incorporates metal price and geological uncertainties jointly. While bidirectional arcs exist between nodes for block i across G realizations, if the nodes are merged as explained in the previous section, then the graph G˜ contains one node per period for a given block i. If toT and v(i, t) is the ore Node i in period t, then the capacity of the arc (S, v(i, t)) is equal to {dbV(i, t)dbV(i, t þ 1)}, which is positive because dbV(i, t)4dbV(i, t þ 1). Similarly, if toT and v(i, t) is waste Node i in period t, then the capacity of the arc (v(i, t), R) is equal to |dbV(i, t)dbV(i, t þ 1)|. The capacities of the arcs from the source to ore nodes and from waste nodes to the sink in the final period t ¼ T are equal to |dbV(i, t)|. In order to ensure that Block i removed in period t stays removed in the subsequent periods, the graph G˜ also contains infinite-capacity arcs from Node i in period t to Node i in period t þ 1. Finally, the infinitecapacity arcs representing slope constraints for a particular block i in Period t are also included in this graph. Figure 5 presents a two-dimensional example. Figure 5 shows the arcs of Block 1 in three periods. If Block 1 is removed (mined) in Period 1, no cost is added to the minimum cut. However, if it is removed in Period 2, the arc with capacity 4–3 crosses the cut, which represents a loss of $1. Similarly, if the block is removed in Period 3, the arcs with capacity 4–3 and 3–2 will cross the cut, representing a loss of $2. If this block is left in the ground, then all three arcs will cross the cut, which represents a loss of $4.

4-3

Period 2

This clarifies the structural difference between the related graphs for the stochastic and conventional frameworks mentioned in the literature review.

s

Period 1

sink, and its capacity is the sum of the capacities of the arcs from this block to the sink across all realizations. If a block is identified as an ore block in some realizations and a waste block in other realizations, then source to node and node to sink arcs are maintained with capacities equal to the sum of the capacities of the source to node arcs in some realizations, and equal to the sum of capacities of node to sink arcs in other realizations, respectively.

| - 4| = 4

t

Figure 5 Graph structure in multiple realizations and discounted block values.

constraint, gives the following formulation Maximize: Z ¼

I X T X

cit xit

ð12Þ

i¼1 t¼1

Subject to: xit  xi 0 t p0; i 0 2 xi for tpT; and i 2 V 8t: I X

Qi xit pbt ; 8t:

ð13Þ

ð14Þ

i¼1

xit 2 f0; 1g; i 2 V 8t; and

T X

xit p1; 8i:

ð15Þ

t¼1

where cit ¼ |dbV(i, t)dbV(i, t þ 1)|, bt is the available pit production capacity in Period t in tons of material, and xit is the binary variable representing node i in period t. Equations (12)–(15) present the classic open pit mine production scheduling formulation. However, this formulation solves for the optimal ultimate pit limit and phase design, which is accomplished by exploiting the network structure of the constraint matrix after relaxing the pit production capacity constraint. As such, it facilitates systematic selection of the Lagrangian parameter l to decrease ore block values (Equation (16)), thus helps to solve for ultimate pit limit and to obtain successively smaller phases within that ultimate pit limit. Similarly, the variable xit in the formulation (12)–(15) is strictly used in the context of discounting block values as explained in Section 2.2.3 and Figure 5. Lagrangian relaxation of pit production capacity constraints produces the classical maximum closure problem

MWA Asad and R Dimitrakopoulos—Implementing a parametric maximum flow algorithm

of Equations (1)–(3). If ltX0 are the multipliers associated with Equation (14), then a modification of the relaxed formulation given in Equations (5)–(7) is given by

Maximize: ZðlÞ ¼

I X T X

½½cit  lt Qi xit  þ

i¼1 t¼1

T X

Subject to: xit  xi 0 t p0; i 0 2 xi for tpT; and i 2 V 8t:

xit 2 f0; 1g; i 2 V 8t; and

T X

xit p1; 8i:

ð17Þ

ð18Þ

The solution to the relaxed problem is obtained by applying the parametric maximum flow (minimum cut) algorithm over a number of iterations in which lt values are selected systematically so that the optimum pit limit and phase design are identified.

2.3.1. The subgradient method. The relaxed formulation describes a maximum closure problem on graph G˜, where the value of ore Node i in Period t is updated over a number of iterations by changing the value of lt in [citltQi]. When lt ¼ 0 the largest possible pit Yl, that is, the ultimate pit limit is obtained and an increase in lt over a number of iterations generates successive pits/phases satisfying Y1DY2D?DYl. From Equation (16) Z(l) is a piecewise linear convex function. At each point lP t, there exists a subgradient st with components stt ¼ bt Ii Qixit(ltt), where xit(lt) is the ith component of the solution vector for all values of t in Equations (16)–(18). The natural approach for systematically updating lt and the node values in Period t is the subgradient method. At the first iteration, the subgradient method initializes the multipliers lt0 to 0. Then, the multipliers are updated in successive iterations t as follows (Tachefine and Soumis, 1997):

such that; dt ¼

   dt stt t ¼ max 0; lt  m r½Zðlt Þ  Z jjst jj2

; with 0oro2:

PI m¼

t¼1

lttþ1

function corresponding to a feasible solution is used. In addition, a procedure to address the gap problem is proposed herein, which employs scaling lt by m, a possible number of phases, defined as:

lt bt ð16Þ

t¼1

ð19Þ

ð20Þ

P take negative values over Since stt ¼ bt Ii Qixit(ltt), stt may P I t successive iterations until Then, i Qixit(lt)pbt. t lt þ 1 in Equation (19) increases and cit ¼ |dbV(i, t) dbV(i, t þ 1)| decreases. This corresponds to ore blocks in Equation (16) and ultimately to the smallest possible pit Y1. The optimal solution Z in Equation (20) is not known in advance. The best known value of the objective

9

i¼1

PT

t¼1 Qi xit ðl0 Þ : bt Tn

ð21Þ

Equation (21) ensures that m remains constant throughout the iterative process as it is a function of the total amount of material (tons) in the ultimate pit (at the first iteration), the user-defined pit production capacity and the number of years, n, in which a particular phase may be mined. The parameter m is used in Equation (19), which is modified form of the original equation (l0 t þ 1 ¼ max[0, [lttdtstt]]) to determine l0 t þ 1 and it controls the step size, dt, so that ltt þ 1 is restricted to the values defined by bt and n, ensuring convergence of the algorithm toward consistent phase sizes within the ultimate pit limit. However, depending upon the orebody model, some instances may require far more iterations than others to develop consistent size phases, which may be mitigated by relaxing the tolerance on phase size using m to further increase or decrease step size and facilitate convergence. The stopping criterion in this iterative process is based on the size of the smallest pit (the first phase, Y1) depending on the pit production capacity and the number of years n in which this particular phase may be mined.

2.3.2. Steps for implementing the algorithm. The following steps facilitate implementation of the parametric maximum flow push-relabel algorithm for developing phase design and ultimate pit limit: 1. Construct the graph so that all arcs are of infinitecapacity as described in Section 2.2.3. 2. Set initial l values to zero at iteration t. 3. Update the discounted value of ore Block i to dbVgti ¼ [dbVgtiltQi]. 4. Merge block values over G grade realizations to compute dbV(i, t), as described in Section 2.2.2. 5. Determine node value, so that cit(l) ¼ |dbV(i, t) dbV(i, t þ 1)| if toT and cit(l) ¼ |dbV(i, t)| if t ¼ T. 6. Update the graph by appending the arcs from source to nodes if cit(l)40 and from nodes to sink if cit(l)p0. 7. Solve this graph using the maximum flow push-relabel algorithm. If t ¼ 1, then the solution vector, xit(ltt), is the optimum ultimate pit limit. Otherwise, the solution is a phase. 8. If the amount of material (tons) inside the pit or phase is less than or equal to the pit production capacity, stop iterating. Otherwise go to Step 9.

10

Journal of the Operational Research Society

P 9. Compute the subgradient stt ¼ bt Ii Qixit(ltt). 10. Evaluate dt and m. 11. Set t ¼ T þ 1 and evaluate new l values. Go to Step 3.

Level 3

Level 2

Level 1

a x

b

c e

d

1

2

3

4

5

6

7

8

9

3. Case study: application to a copper deposit The application to a copper deposit demonstrates the practical aspects of the proposed approach. This application presents two streams of results of stochastic approach: (i) ‘Case A’ in which l values were scaled using Equation (19); and (ii) ‘Case B’ in which l was not scaled but modified instead using the equation ltt þ 1 ¼ max[0, [lttdtstt]]. Apart from evaluation of Cases A and B in a stochastic approach, this application also presents a comparison of the proposed stochastic and conventional approaches. For the proposed stochastic approach, the procedure commences with the development of the copper orebody model that constitutes time-dependent discounted values incorporating the uncertainty in metal content and commodity price jointly. The ultimate pit limit and phase design follows the construction of a graph representing this stochastic framework and its solution implementing the steps presented in Section 2.3.2. The geological database of the deposit includes 185 drill holes on a pseudo regular grid of 50 m by 50 m over an approximately rectangular area of 1600  900 m2. Spatial modelling of the geological information generates 20 grade realizations using the technique of conditional simulations (Leite, 2008; Boucher, 2009; Horta and Amilcar, 2010). Similarly, the uncertainty in commodity price was modelled by generating 20 000 realizations of copper price using Equation (8). Fixing the block tonnage (Q) at 10 800 tons, the initial metal price (P) at $2.00 per lb. of copper, mining cost (m) at $2.00 per ton of material, processing cost (p) at $10.00 per ton of ore, refining or selling cost (r) at $0.1 per lb. of copper, metallurgical recovery (y) at 90% and discount rate at 8%, the procedure outlined in Section 2.1 produces an orebody model consisting of time-dependent discounted block values spanning over 20 simulations of the orebody. Generally, orebody models consist of several hundred thousand blocks. The orebody model in this case study consists of 39 244 mining blocks and block values have been discounted over 12 periods. As described in Section 2.2.3, the graph is created considering the discounted block values, where ore blocks are connected to the source node and waste blocks are connected to the sink node. With a slope angle at 451, we follow the block precedence pattern presented in Figure 6, which requires removal of 14 overlying blocks, that is, Blocks a, b, c, d, e at Level 2 and Blocks 1, 2, 3, 4, 5, 6, 7, 8, 9 at Level 1 before an underlying block x at Level 3 is considered for extraction. Infinite capacity arcs maintaining the extraction precedence are

Figure 6 Precedence constraint pattern.

included in the graph. To ensure that a block is removed only once, infinite capacity arcs between the nodes in various periods are included in the graph. Therefore, in this case study, the algorithm solves a graph constituting 470 928 (number of mining blocks multiplied by the number of periods, ie, 39 244  12 ¼ 470 928) nodes, 552 476 arcs from source to (ore blocks) nodes and (waste blocks) nodes to sink, and 6 716 308 infinite capacity arcs. The infinite capacity arcs representing precedence and block mining constraints remain the same for a particular orebody model, but the arcs from source to ore blocks and from waste blocks to sink are updated in successive iterations with a decrease in the value of ore blocks as a function of l. A parametric maximum flow push-relabel algorithm is applied according to the steps outlined in Section 2.3.2. Maintaining a pit production capacity of 27 million tons per period, Equation (21) controls the size of phases leading to the development of six optimal phases within the ultimate pit limit. As explained in Section 2.3.1, the algorithm utilizes m and repeats until a consistent size phase is generated. For the instance presented here, the algorithm converges in single iteration for Phases 1, 4, 5, and 6; while it converges in 7 and 16 iterations for Phases 2 and 3, respectively. This difference in number of iterations corresponds to the production capacity, which is controlling the consistency in phase sizes by varying m and l values. Cases A and B present a comparative analysis to evaluate the effectiveness of the l scaling in the subgradient method. Table 1 compares the phases and the ultimate pit limit in both cases. Figure 7 presents section maps of the ultimate pit limit and the sequence of phases. Figure 8 presents the minimum, maximum, and average stripping ratios of phases in both cases. The stripping ratio is an economic indicator that defines the amount of overlying waste to be mined per ton of ore or it is the ratio of the number of waste blocks and the number of ore blocks in a particular mining phase or pit. The life of the operation is approximately 11 years (Table 1) in both cases. In Case A the phase sizes are relatively uniform. The scaling of l values has thus effectively addressed the gap problem. The optimal solution in Case B is mining 3019 blocks in Phase 1 over a period of approximately 1 year followed by mining the remaining

MWA Asad and R Dimitrakopoulos—Implementing a parametric maximum flow algorithm

11

Table 1 Phase design and ultimate pit limit solutions in Cases A and B Case

Phase Number

A: With scaling of the l values

B: Without scaling of the l values

Blocks

Pit

Tons

Life (years)

Number

1

2413

26 060 400

0.97

1

2413

2

5773

62 348 400

2.31

2

8186

3

4055

43 794 000

1.62

3

12 241

4

5720

61 776 000

2.29

4

17 961

5

4771

51 526 800

1.91

5

22 732

6 Total

4887 27 619

52 779 600 298 285 200

1.95 11.05

Ultimate pit

27 619

1 2 Total

3019 24 600 27 619

32 605 200 265 680 000 298 285 200

1.21 9.84 11.05

1 Ultimate pit

3019 27 619

Case A

Case B

5

5

10

10

15

15

20

20

25

1

30

10

2

20

Blocks

3

30

4

40

5

50

25

6

60

70

1

30

10

2

20

30

40

50

60

70

Figure 7 A representative section showing phase and ultimate pit solutions in Cases A and B.

Minimum

9 8 7 6 5 4 3 2 1

Case B

Maximum Average

5 Stripping Ratio

Stripping Ratio

Case A Minimum

Maximum

Average

4 3 2 1

1

2

3 4 5 Phase (Pushback)

6

1

2 Phase (Pushback)

Figure 8 Minimum, maximum and average phase stripping ratios in Cases A and B.

24 600 blocks from Phase 2 over the remainder of the operation. This solution presents a gap problem between Phase 1 and Phase 2, which is a large discrepancy in phase size. Figures 7 and 8 support the results listed in Table 1. Figure 7 presents an even distribution of material production among six phases in Case A and highly variable phase

material amounts in Case B. In Case A the incremental stripping ratio increases smoothly (Figure 8) from Phase 1 to Phase 6. In Case B, the increase is abrupt. Smooth production spanning over a number of phases as in Case A may potentially bring higher discounted cash flows during the early years of the open pit mining operation.

12

Journal of the Operational Research Society

Apart from addressing the gap problem, this application also demonstrates a comparison of the stochastic approach presented in this paper and the corresponding conventional equivalent shown in Figure 2. For the conventional approach, an ‘E-type’ orebody model is generated by averaging 20 grade realizations. Again, fixing the block tonnage (Q) at 10 800 tons, the metal price (P) at $2.00 per lb. of copper, the mining cost (m) at $2.00 per ton of material, the processing cost (p) at $10.00 per ton of ore, the refining or selling cost (r) at $0.1 per lb. of copper, and the metallurgical recovery (y) at 90%, the block values are calculated using a modified Equation (9) as follows: bVi ¼ ðP  rÞQi gi y  mQi  pQi

ð22Þ

where, bVi corresponds to the value of block i, which is adjusted such that 8 if bVi 40; ie; i is an ore block bVi > > < in the E-type model; bVi ¼ mQi if bVi p0; ie; i is an waste block > > : in the E-type model: A graph is created as described in Section 2.2.1, which consists of 39 244 nodes, 39 244 arcs from source to (ore block) nodes and (waste block) nodes to sink, and 770 389 infinite capacity arcs. Using the scaling option, a maximum flow push-relabel algorithm solves this graph according to the steps outlined in Section 2.3.2, which develops three optimal phases to the ultimate pit limit. Phases 1, 2, and 3 constitute 4512, 4181, and 6111 blocks, respectively; that is, a total of 14 804 blocks within the optimal ultimate pit limit. Therefore, the size of the ultimate pit limit is reduced by 46.40% as compared to the stochastic approach proposed herein. A reduction in the size of ultimate pit or the scale of operation has major impact on the overall monetary value, amount of metal to be produced, and life of operation. These results show that jointly incorporating uncertainty in metal content and commodity price in a stochastic framework promises substantial improvement as compared to existing conventional solutions.

4. Conclusions This paper proposes a method for developing ultimate pit limit and phase design in long-term open pit mine planning that jointly incorporates metal content and commodity price uncertainties. The three-stage optimization procedure utilizes a set of price and grade simulations to determine the time-dependent discounted block values, creates a graph based on this stochastic framework, and implements the parametric maximum flow algorithm. Lagrangian relaxation of pit production capacity constraints and the subgradient method facilitate selection of the l values for developing successive phases.

The case study presented an example in which the stochastic approach proposed develops an optimal design with an approximately 45% larger pit, leading to a higher economic value and metal production compared to the conventional approach. These differences are consistent with past findings for mine design based on stochastic optimization approaches. The proposed stochastic network flow approach has no particular size limits and can accommodate millions of nodes and arcs. Ongoing research aims to both further test the method presented above, as well as include variability and stochasticity in other parameters such as mining and processing costs, and discount rate. Acknowledgements —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, BHP Billiton, De Beers, Newmont, and Vale. Thanks are in order to Brian Baird, Peter Stone, Darren Dyck, and Gavin Yates of BHP Billiton for their support, collaboration, and technical comments.

References Ahuja RK and Orlin JB (1989). A fast and simple algorithm for the maximum flow problem. Operations Research 37(5): 748–759. Akaike A and Dagdelen K (1999). A strategic production scheduling method for an open pit mine. Proceedings of the 28th International Symposium on Application of Computers and Operations Research in the Mineral Industry, Colorado School of Mines: Golden, Colorado, pp 729–738. Asad MWA (2007). Optimum cutoff grade policy for open pit mining operations through net present value algorithm considering metal price and cost escalation. Engineering Computations 24(7): 723–736. Bascetin A and Nieto A (2007). Determination of optimal cut-off grade policy to optimize NPV using a new approach with optimization factor. Journal of the South African Institute of Mining and Metallurgy 107(2): 87–94. Boucher A (2009). Sub-pixel mapping of coarse satellite remote sensing images with stochastic simulations from training images. Mathematical Geosciences 41(3): 265–290. Bower U, Geis A and Winkler A (2007). Commodity price fluctuations and their impact on monetary and fiscal policies in Western and Central Africa. European Central Bank, Occasional Paper Series, No. 60. Caccetta L and Hill S (2003). An application of branch and cut to open pit mine scheduling. Journal of Global Optimization 27(2–3): 349–365. Chandran BG and Hochbaum DS (2009). A computational study of the pseudoflow and push-relabel algorithms for the maximum flow problem. Operations Research 57(2): 358–376. Chen MH (2010). Understanding world metals prices—Returns, volatility and diversification. Resources Policy 35(3): 127–140. Coleou T (1989). Technical parameterization of reserves for open pit design and mine planning. Proceedings of the 21st International Symposium on Application of Computers and Operations Research in the Mineral Industry, Las Vegas, Nevada, pp 485–494. Denby B, Schofield D and Surme T (1998). Genetic algorithms for flexible scheduling of open pit operations. Proceedings of the 27th Symposium on Application of Computers and Operations

MWA Asad and R Dimitrakopoulos—Implementing a parametric maximum flow algorithm

Research in the Mineral Industry, Royal School of Mines: London, UK, pp 605–616. Dimitrakopoulos R (2011). Stochastic optimization for strategic mine planning: A decade of developments. Journal of Mining Science 47(2): 138–150. Dimitrakopoulos R and Abdel Sabour SA (2007). Evaluating mine plans under uncertainty: Can the real options make a difference? Resources Policy 32(3): 116–125. Dowd PA (1994). Risk assessment in reserve estimation and openpit planning. Transactions of the Institute of Mining and Metallurgy 103: 148–154. Faaland B, Kim K and Schmitt T (1990). A new algorithm for computing the maximal closure of a graph. Management Science 36(3): 315–331. Francois-Bongarcon D and Guibal D (1984). Parameterization of optimal designs of an open pit—Beginning of a new phase of research. Transactions of Society of Mining Engineers of AIME 274: 1801–1805. Gallo G, Grigoriadis MD and Tarjan RE (1989). A fast parametric maximum flow algorithm and applications. SIAM Journal of Computing 18(1): 30–55. Godoy M and Dimitrakopoulos R (2011). A risk quantification framework for strategic mine planning: Method and application. Journal of Mining Science 47(2): 235–246. Goldberg AV and Tarjan RE (1988). A new approach to the maximum flow problem. Journal of the Association for Computing Machinery 35(4): 921–940. Hochbaum DS (2001). A new-old algorithm for minimum-cut and maximum-flow in closure graphs. Networks 37(4): 171–193. Hochbaum DS (2003). Efficient algorithms for the inverse spanning tree problem. Operations Research 51(5): 785–797. Hochbaum DS (2004). Selection, provisioning, shared fixed costs, maximum closure, and implications on algorithmic methods today. Management Science 50(6): 709–723. Hochbaum DS (2008). The pseudoflow algorithm: A new algorithm for the maximum-flow problem. Operations Research 56(4): 992–1009. Hochbaum DS and Chen A (2000). Performance analysis and best implementations of old and new algorithms for the open-pit mining problem. Operations Research 48(6): 894–914. Horta A and Amilcar S (2010). Direct sequential co-simulation with joint probability distributions. Mathematical Geosciences 42(3): 269–292. Hustrulid W and Kuchta M (2006). Open Pit Mine Planning and Design. 2nd edn. Taylor and Francis/Balkema: The Netherlands. Johnson TB (1968). Optimum open pit mine production scheduling. PhD Thesis, Department of IEOR, University of California, Berkeley, CA. King B (2011). Optimal mining practice in strategic planning. Journal of Mining Science 47(2): 247–253.

13

Lane K (1964). Choosing the optimum cutoff grade. Colorado School of Mines Quarterly 59(4): 811–829. Lane K (1988). The Economic Definition of Ore, Cutoff Grade in Theory and Practice. Mining Journal Books: London. Leite A (2008). Stochastic optimization approaches to open pit mine planning: Applications for and the value of stochastic approaches. MEng Thesis, Mining and Materials Engineering, McGill University, Canada. Lerchs H and Grossmann IF (1965). Optimum design of open pit mines. Transactions CIM 58: 17–24. Meagher C, Abdel Sabour SA and Dimitrakopoulos R (2009). Pushback design of open pit mines under geological and market uncertainties. Proceedings of the International Symposium on Orebody Modeling and Strategic Mine Planning: Old and New Dimensions in Changing World, pp 297–304. Picard JC (1976). Maximal closure of a graph and applications to combinatorial problems. Management Science 22(11): 1268–1272. Ramazan S (2007). The new fundamental tree algorithm for production scheduling of open pit mines. European Journal of Operational Research 177(2): 1153–1166. Ramazan S and Dimitrakopoulos R (2012). Production scheduling with uncertain supply: A new solution to the open pit mining problem. Optimization and Engineering, DOI: 10.1007/s11081012-9186-2. Schwartz ES (1997). The stochastic behavior of commodity prices— Implications for valuation and hedging. Journal of Finance 52(3): 923–973. Tachefine B and Soumis F (1997). Maximal closure on a graph with resource constraints. Computers and Operations Research 24(10): 981–990. Wang Q and Sevim H (1993). An alternative to parameterization in finding a series of maximum-metal pits for production planning. Proceedings of the 24th International Symposium on Application of Computers and Operations Research in the Mineral Industry, Canadian Institute of Mining and Metallurgy and Ecole Polytechnique, Montreal: Quebec, Canada, pp 168–175. Whittle J (1988). Beyond optimization in open pit design. Proceedings of the Canadian Conference on Computer Applications in the Mineral Industries, Quebec: Canada (Balkema, Rotterdam), pp 331–337. Whittle J (1999). A decade of open pit mine planning and optimization—The craft of turning algorithms into packages. Proceedings of the 28th International Symposium on Application of Computers and Operations Research in the Mineral Industry, Colorado School of Mines: Golden, Colorado, pp 15–24.

Received August 2010; accepted February 2012 after two revisions