A local branching heuristic for the open pit mine production scheduling problem

Accepted Manuscript

A local branching heuristic for the open pit mine production scheduling problem Mehran Samavati , Daryl Essam , Micah Nehring , Ruhul Sarker PII: DOI: Reference:

S0377-2217(16)30541-0 10.1016/j.ejor.2016.07.004 EOR 13830

To appear in:

European Journal of Operational Research

Received date: Revised date: Accepted date:

14 July 2015 4 July 2016 4 July 2016

Please cite this article as: Mehran Samavati , Daryl Essam , Micah Nehring , Ruhul Sarker , A local branching heuristic for the open pit mine production scheduling problem, European Journal of Operational Research (2016), doi: 10.1016/j.ejor.2016.07.004

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.

ACCEPTED MANUSCRIPT Highlights A local branching heuristic is used to tackle the open pit mine scheduling. An adaptive branching scheme is developed and combined with local branching. A heuristic is developed to generate starting feasible solutions for local branching. Minimum resource requirements are considered. A comparison is drawn between the proposed methodology and other techniques.

AC

CE

PT

ED

M

AN US

CR IP T

    

1

ACCEPTED MANUSCRIPT A local branching heuristic for the open pit mine production scheduling problem Mehran Samavati a*, Daryl Essam a, Micah Nehring b, and Ruhul Sarker a a

School of Engineering and Information Technology, University of New South Wales, Canberra, Australia

b

School of Mechanical and Mining Engineering, University of Queensland, Brisbane, Australia *Corresponding Email: [email protected]

CR IP T

Abstract

M

AN US

This paper considers the well-known open pit mine production scheduling problem (OPMPSP). Given a discretisation of an orebody as a block model, this problem seeks a block extraction sequence that maximizes the net present value (NPV) over a horizon of several periods. In practical applications, the number of blocks can be large, and therefore, the problem can be difficult to solve. It is even more challenging when it incorporates minimum resources requirements that are represented as lower bounds on resource constraints. In this study, we propose to tackle OPMPSP by using a novel metaheuristic technique known as local branching. To accelerate the search process, we combine local branching with a new adaptive branching scheme, and also develop a heuristic to quickly generate a starting feasible solution. Despite consideration of minimum requirements being seldom taken into account in the literature, this method yields near-optimal solutions for a series of data sets we have conceptually generated. To judge the performance of our methodology, the results are compared to those of two techniques from the literature, as well as to those obtained by a mixed integer linear programming (MILP) solver.

1. Introduction

ED

Keywords: Open pit mine production scheduling; Minimum resource requirements; Local branching; Heuristics

AC

CE

PT

Open pit mining is the process of first removing valuable material from the earth‟s surface and then selling it in the market to derive profit. The economic viability of a mining project significantly relies on careful management. Operations research has been extensively used to make decisions on how the extraction process should be accomplished in order to obtain the best profit. To make the process of planning for material extraction more tractable, mine planners categorise the orebody into individual small cubic (but not always cubic) parts called blocks. Thus, a block model is an orebody that has been discretised into blocks. While a number of models are available, the 3D fixed-block model has attracted the most attention, as it best suits the application of computerized optimization techniques (Caccetta & Hill, 2003). A cut-off grade is the level below which material within an orebody does not contain sufficient value to economically justify excavation. Blocks that contain valuable mineral above the cut-off grade will be sent to a processing plant. All other blocks will be sent to a waste dump. Both types of blocks are scheduled across the life of the mining operation. The objective of this strategic mine plan is to maximize net present value (NPV) while complying with a variety of operational restrictions. In the planning of mining operations, one of the major issues is to determine the optimal schedule for block extractions in a certain block model. This gives rise to the open pit mine production scheduling problem (OPMPSP), where the objective is to find which block should be extracted, if at all, in which time period, in order to gain the highest NPV. Several operational restrictions must be taken into 2

ACCEPTED MANUSCRIPT account while solving OPMPSP using operations research techniques. For example, precedence constraints, known as geometric sequencing constraints, or wall slope constraints, ensure that blocks are excavated such that the deposit does not collapse inwards. This means that a given block cannot be excavated before its overlying blocks. Resource capacities enforce limitations on the amount of material removed from the orebody and the amount of ore (valuable material) processed at the mill. The time horizon of the project is subdivided into time periods, each having a certain amount of resource-processing and resource-production, known as the maximum processing capacity and the maximum production capacity, respectively.

CR IP T

Minimum mining and ore processing requirements are other constraints in OPMPSP. These are important operational considerations that are principally associated with consistent acceptable equipment utilisation between scheduling periods. Minimum mining and processing rates in each period across set scheduling horizons also allows other resources associated with operating heavy equipment, such as labour and maintenance, to be appropriately allocated. This ensures overall operational feasibility.

PT

ED

M

AN US

Since the 1960s, numerous researchers have been addressing OPMPSP. Early studies are only focused on the ultimate pit limit problem, a simplified version of the problem in which only precedence constraints are considered (Underwood & Tolwinski, 1998). Different variants of OPMPSP are subsequently formulated as integer or mixed-integer linear programs (MILP). The corresponding mathematical models involve binary variables that determine whether or not a given block is removed in a certain time period. In MILPs, additional continuous variables denote the amount of each block that is sent to a particular destination in a certain time period. The constraints of the formulations vary depending on the assumed operational restrictions, and thereby, form variants of the scheduling problem. In order to make the problem more tractable for solution, OPMPSP is often addressed with certain simplifications, such as fixed cut-off grade and no inventory. This simplified version of the problem is known as the constrained pit limit problem (CPIT), where the corresponding formulation involves binary variables, and the objective is to maximize NPV subject to precedence constraints, maximum resource capacities, and minimum resource requirements. The model will be explained in detail in Section 3.

AC

CE

Although minimum resource requirements are important, they can substantially add complexity to a problem. Cullenbine et al. (2011) conduct an experiment to analyse the influence of minimum capacities on complexity of CPIT. They find that a particular instance can be solved to optimality by CPLEX within 176 seconds, while it requires 3688 seconds when considering minimum requirements. There are a number of reasons that minimum requirements make the problem more complicated to solve. Below are some of them: Firstly, maximum processing and production capacities normally appear as two separate capacity (knapsack) constraints in the mathematical model. The minimum requirements result in lower bounds on these constraints, which may conflict with the knapsack‟s upper bounds. For instance, a lower bound in a processing knapsack constraint may be in direct conflict with an upper bound on extraction capacity. Secondly, when only an upper bound is considered, the problem can be seen as the extraction of high grade blocks (high value blocks) as early as possible, in order to minimise money value effects. However, when a lower bound is integrated, a look-ahead approach is needed to keep ore material for all periods. 3

ACCEPTED MANUSCRIPT Finally, certain heuristics and techniques that are used for the problem cannot be applied when considering minimum requirements. For instance, some heuristics can reduce the size of the model a priori, by first solving the ultimate pit limit problem, and then incorporating only those blocks that are present in the solution of this problem into their solution algorithm. However, such a size reduction is not possible in a model incorporating minimum requirements, as a lower bound on processing can require a non-economical block to be sent to the mill so as to preserve the feasibility of the instance.

AN US

CR IP T

In spite of their importance, minimum resource requirements are addressed in very few studies. One of the earliest methods that takes these requirements into consideration is the branch and cut (B&C) algorithm by Gaupp (2008). This algorithm uses special cutting plans produced by using the earliest and latest possible extraction time available for every block. The author develops a variable reduction technique to reduce the number of variables before applying the B&C algorithm. Gaupp also exploits Lagrangian relaxation in a separate attempt to reduce solution times. Even though their results by both B&C and Lagrangian relaxation are substantially better than those obtained from direct solution of the entire problem, the largest instance addressed consists of only 10,819 blocks and 6 time periods. Other methods are the algorithms of Cullenbine et al. (2011) and Lambert and Newman (2014). The former develops a combination of a sliding time window heuristic with Lagrangian relaxation to solve CPIT, and the latter uses a sliding time window heuristic to generate an initial feasible solution for the use of a tailored Lagrangian relaxation procedure. Both of these studies obtain good results for instances of up to 25,620 blocks and 10 time periods. A relatively small case study with 7,020 blocks and 5 periods is also addressed by Kumral (2012).

PT

ED

M

Therefore, we propose a new efficient technique to address larger instances of OPMPSP while including both minimum requirements and maximum capacities. Our technique uses a new high-level solution strategy, presented by Fischetti and Lodi (2003), called local branching. This strategy is, in fact, a hybrid metaheuristic that combines mixed integer programming (MIP) solution techniques with typical metaheuristic concepts, such as local searches, neighbourhood definitions, and intensification and diversification mechanisms. Local branching uses a generic MIP solver to explore feasible subspaces (neighbourhoods) that are constructed by introducing linear inequalities known as local branching constraints. We now give a brief description of this technique.

+

AC

*

CE

We begin by considering the following general 0-1 linear programming problem:

where , , and . Let be a feasible solution (referred to as a reference * +, solution) to , { }, and be a positive integer. Also, let ̃ be a starting incumbent solution, and bestLB ̃ be a starting lower bound. Having inequality ( ) ∑ ( ) ∑ , the feasible region of can be divided by taking, on ) one branch (referred to as left branch), the values of by which inequality ( is satisfied, ) and on the other branch (referred to as right branch), the values of by which inequality ( is satisfied. The subproblem generated in left branch is then optimised using a generic solver, and a new solution ̃ is obtained. Let ̅ ̃, and set (i.e., the best solution obtained by local branching) to . In the right branch‟s subproblem, constructed by reversing the local branching ) ) constraint ( into ( , if ̃ > best , then ̃ and best ̃. 4

ACCEPTED MANUSCRIPT Subsequently, the process is iterated, that is, the local branching constraint ( model .

̅)

is added to

To our best knowledge, there is no published research that utilizes local branching for addressing OPMPSP. In our approach, we develop a new adaptive branching scheme, and combine it with local branching. This strategy makes local branching extremely fast and efficient. We also develop a new heuristic to find a starting feasible solution to further accelerate the procedure.

CR IP T

To judge the efficiency of the proposed methodology, we develop a series of conceptually generated data sets and compare the results of solving them by our methodology to those obtained by other techniques from the literature: B&C and Lagrangian relaxation. The proposed algorithm shows consistently better performance compared to other techniques. Interestingly, our methodology tackles an instance of 50,383 blocks and 10 time periods and obtains a near-optimal solution while, as aforementioned, the largest instance solved in the literature includes only 25,620 blocks and 10 time periods.

AN US

The remainder of this paper is organised as follows. In the next section, we review relevant literature and summarises our motivation for this work. In Section 3, we present the mathematical programming model considered in this study, and in the following section, explain our solution methodology. Section 5 raises experimental results. Section 6 concludes the work and provides some directions for future studies. 2. Literature review

M

This section briefly discusses the relevant literature for OPMPSP, as well as the local branching procedure. 2.1. OPMPSP

AC

CE

PT

ED

Due to the fact that the size of OPMPSP in practice does not allow MIP solvers to solve it in its entirety, a variety of aggregation, relaxation, and heuristic approaches have been developed for different variants of this problem in the literature. One of the methods is the “Fundamental Tree Algorithm” by Ramazan (2007), which considers fixed cut-off grade, minimum processing requirements, and maximum processing and production capacities. This study reduces the number of variables by aggregating blocks. Boland et al. (2009) address OPMPSP with variable cut-off grade and maximum resource capacities. They use aggregated blocks to schedule the extraction process, while individual blocks are used for processing decisions. Another algorithm is the “aggregation and disaggregation” heuristic of Jélvez et al. (2016). While aggregation techniques substantially reduce solution times, the obtained schedules are not, however, practically applicable from the mining operational point of view when it comes to annual ore and/or waste production, average grade of ore, and so on. A combination of aggregation and Lagrangian relaxation is proposed by Kawahata (2007) to address variable cut-off grade alongside maximum capacities. Some researchers have recently focused their efforts towards heuristics that are based on topological ordering. The most powerful proposed heuristics are the well-known TopoSort heuristics by Chicoisne et al. (2012) and Moreno et al. (2010). The heuristic by Liu and Kozan (2016) is also a good example of topological ordering-based methods. These recent studies have obtained impressive results for CPIT constrained by knapsack constraints with lower bounds of 0.

5

ACCEPTED MANUSCRIPT There are very few studies that employ metaheuristics to address OPMPSP. The most recent and efficient metaheuristic is developed by Shishvan and Sattarvand (2015), which is based on Ant Colony Optimization (ACO). Solving the LP relaxation of CPIT has also been taken into account. The solution procedures of Moreno et al. (2010) and Bienstock and Zuckerberg (2010) can be considered to be the most effective ones for solving the LP relaxation of CPIT. For a detailed literature review on the solution approaches, we refer the reader to Gaupp (2008) and Espinoza et al. (2013). 2.2. Local branching

AN US

CR IP T

The local branching technique was first developed by Fischetti and Lodi (2003) with the aim of addressing complex and large MIPs. The authors prove the efficiency of this technique by solving several difficult MIP problems–network design, crew scheduling, railway crew scheduling, nesting, telecommunication network design, lot-sizing, rolling stock and line planning, and railway line planning. Fischetti et al. (2004) combine local branching with variable neighborhood search (VNS) and apply this methodology to a telecommunication network design problem. Another VNS heuristic for solving MIP problems is described by Hansen et al. (2006). The resource-constrained projectscheduling problem has also been addressed by using local branching. The sole example of this is the combination of local branching with a special variable reduction method by Zhu et al. (2006). The superiority of local branching compared to the best relevant heuristics for solving a network design problem, is demonstrated by Rodríguez-Martín and Salazar-González (2010). 3. Mathematical formulation

ED

M

The mathematical formulation considered in this study is associated with the version of OPMPSP that involves a predefined cut-off grade, minimum mining and processing requirements, and maximum capacity restrictions. Each block only requires one time period to extract, and cannot be partially removed due to practical implementation issues. OPMPSP can be formulated as follows:

AC

CE

PT

Sets and parameters set of blocks * + set of ore blocks * + set of waste blocks * + set of time periods * + set of immediate predecessors for block tonnage of block maximum mining capacity (tons) minimum mining requirement (tons) maximum processing capacity (tons) minimum processing requirement (tons) net present value obtained from block Decision variables if block is removed in time period ; otherwise

Mathematical model

( ) Subject to: ∑ ∑

∑

∑

∑

(1)

∑

(2) (3)

∑

(4) (5)

6

ACCEPTED MANUSCRIPT *

+

(6)

The objective function seeks to maximize NPV. Constraint (2) ensures that no block is extracted before its predecessors. Constraint (3) guarantees that the minimum mining requirement is satisfied in each time period, and that the maximum mining capacity is not exceeded. Constraint (4) performs the same task for lower and upper bounds on processing constraints in each time period. Constraint (5) ensures that each block is removed, if at all, in only one time period. Equation (6) is a binary constraint.

CR IP T

4. Solution methodology

AN US

In this section, we propose our methodology based on the combination of the local branching technique and a new adaptive branching scheme. In the first subsection, we briefly describe local branching and explain the way we employ it. In the second subsection, we develop an adaptive branching scheme to combine with local branching, and in the following subsection, propose our methodology. Finally, we develop a heuristic to generate starting solutions for this methodology in order to increase its efficiency. For convenience of reading, a summary of the notation used throughout this section is given here: Notation

M

the entire successor set of block (including block ) the entire precedence set of block (including block ) earliest possible extraction time period for block latest possible extraction time period for block the set of blocks in the top levels of the block model

ED

4.1. Local branching in open pit mine production scheduling

CE

PT

In model , each integer feasible solution has a maximum of variables that are equal to 1. ̅ *( + indicate the binary support ) Given a feasible reference solution ̅ of model , let ̅ of ̅ . For a given positive integer parameter , we define the left branch‟s subproblem of ̅ as the set of the feasible solutions of model satisfying the additional local branching constraint: (

̅)

∑(

(̅

)

)

∑(

)

̅

(7)

AC

By adding Inequality (7) to the current solution space, the left branch‟s subproblem is generated. This means that in this subproblem at least ̅ blocks are extracted in the same time period as the reference solution ̅ . The right branch‟s subproblem is generated by reversing Inequality (7) to: (

̅)

(8)

The subproblem created by Inequality (7) is very restrictive for OPMPSP due to the precedence constraints. Suppose that is set to 1 in the left branch‟s subproblem according to the support set ̅. This means that every block in set must be extracted by period 3, and none of the blocks in set can be extracted earlier than this period. As a result, many decision variables in the solution of this subproblem have the same value as the solution in the parent node. To cope with + be the this issue, we define a new support set for the use of local branching. Let * ̅ reference solution. We then express the new support set as ̅ = {( ) ̅ }. 7

ACCEPTED MANUSCRIPT 4.2. Adaptive branching scheme In this section, we develop an adaptive branching scheme. In the next section, we combine it with the local branching technique. The earliest possible extraction time (ES) of block is defined as the earliest time period such that this block can be removed. Similarly, the latest possible extraction time (LS) of block is defined as the latest time period such that the block can be extracted. Given ES and LS for each block, bounds on the earliest and latest extraction time ( ) are defined for each to be used in solution techniques.

AN US

CR IP T

Because of Constraint (5) in model , the variables naturally divide into type-1 special ordered sets (Beale & Forrest, 1976), and so the model is suitable for the use of the special ordered set (SOS) branching technique. By fixing some variables, SOS branching reduces the number of variables and constraints. In large instances, however, it may not significantly improve upon the branching-based techniques, e.g., local branching. In response, we develop an adaptive branching scheme that removes variables from the model by not only the SOS technique but also using the time window ( ). We describe this scheme in Algorithm 1.

M

Algorithm 1 is run every time whenever CPLEX calls for SOS branching. This algorithm commences with a selected block, which is the first element of the special ordered set in the context of SOS branching technique. We let CPLEX select the special ordered set, or blocks, based on the priority order it generates. Therefore, Algorithm 1 receives the selected block as part of its input. Another input of Algorithm1 is the fractional LP solution that also comes from CPLEX when SOS branching is called. Let ̅ * ̅ + indicates this LP solution.The third input in this algorithm is associated with the weighted extraction time of the selected block. We choose time period t as the ∑ weight to each variable, so ̅ shows the weighted extraction time of block . The idea

ED

of using the time period as the weight has been considered in Chicoisne et al. (2012).

PT

The algorithm begins by splitting the variables into two subsets. The first subset consists of variables that their value in the LP solution is less than . The second subset consists of those with an LP value of larger than . Two branches 1 and 2 are then created by, respectively, setting the variables of the first and second subsets to zero.

AC

CE

Algorithm 1 then sets more variables to zero, based on the earliest possible extraction times. To illustrate this part of the algorithm, a simple example is provided. Let 6, 1, 3.2, and 6. Then, respectively, in branches 1 and 2, variables * + and * + are all set to zero. Clearly, the earliest possible extraction time in branch 1 will equal 4, and because of this, the earliest possible extraction time for all the successors of block 6 may be restricted. Suppose that block 10 is one of the successors of block 6. If the earliest possible extraction time for block 10 is smaller than that of its immediate predecessors–including block 6–then that time should increase. Suppose that blocks 5, 6, and 7 are the immediate predecessors of block 10, and the maximum earliest time belongs to block 6, which was earlier stated to equal 4. Therefore, the earliest time for block 10 should also increase to 4. However, block 10 may not be removable in period 4 due to mining and processing capacities, and so it might be necessary to postpone removing it until the following time period. If so, the corresponding variable of block 10 in period 4, i.e., will set to zero. Algorithm 1 updates the earliest possible extraction times and sets variables to zero in the same way. In this algorithm, four supplementary variables are used to store certain values: carries the value of the earliest possible extraction time for each block at the 8

ACCEPTED MANUSCRIPT beginning of the algorithm. denotes the immediate predecessors of block (i.e., ), and the largest earliest time among them is selected and stored in . As stated earlier, indicates the entire ∑ * precedence set of block . Now, assume that 5. Then the assignment + transfers the total tonnage of the blocks that are in set , and that possess an earliest possible extraction time of 5, to . As mentioned above, indicates the tonnage of block . represents the same as that of , but for ore blocks. Algorithm 1: Adaptive branching scheme. Inputs:

* ̅ +;

; ̅

selected block

∑ ̅

CR IP T

begin for

do if ̅ then set to zero and add it to branch 1; else set to zero and add it to branch 2 endif; if ⌊ ⌋ then ⌋ *

for

+ do

“for each successor of block

*

∑

+

*

“transfer the total tonnage of the blocks that are in set possess an earliest possible extraction time of to ”

+

else

endif; then *

then

endif; if ⌊

“if insufficient capacity”

} to zero and add them to branch 1

CE

endif; endif; set {

“transfer the total tonnage of the ore blocks that are in set and possess an earliest possible extraction time of to ”

PT

if

+

ED

∑

⌋ then

+ } to zero and add them to branch 2

AC

*

and

“if insufficient capacity”

then

if

”

M

if

set { endif;

AN US

⌊

end

4.3. Combination of local branching and adaptive branching scheme In this section, we describe our proposed methodology based on the combination of local branching and the adaptive branching scheme developed in the previous section. The local branching method and the relevant notation have been described earlier. The overall methodology is illustrated in Algorithm 2. In the first step, given a feasible reference solution ̃ and its objective value ( ̃), we set (i.e., the current lower bound) and (i.e., the 9

ACCEPTED MANUSCRIPT NPV of the best solution obtained by the algorithm) to ( ̃), and set (i.e., the best solution obtained by the algorithm) to ̃. Given a positive integer , the algorithm begins by creating the leftbranch subproblem by adding the local branching constraint ( ̅ ) to model , where ̅ ̃ and . As stated earlier, this subproblem is solved using a generic solver, typically CPLEX, in the local branching procedure. In our methodology, however, this subproblem is solved using our branching scheme within a certain time limit (node_time_limit). After each call to the branching scheme, a new solution ̃ is obtained, and one of the following four cases arises with respect to this solution: The current left branch has been solved to optimality. In this case, a new right branch is constructed by reversing the last local branching constraint into ( ̅ ) . Lower bound , best solution , best NPV, best lower bound , and ̅ are updated and the ( ) process is iterated, that is, the local branching constraint ̅ is added to model .

II.

No feasible solution with objective value strictly better than lower bound is found within the node time limit, but there is no guarantee that such a solution does not exist. In this case, Fischetti and Lodi‟s soft or strong diversification scheme is applied depending on the current value of flag diversity. If diversity=false, then the soft diversification scheme is performed in which the last local branching constraint is removed, , local branching constraint ( ̅ ) is added, and the process is iterated. If diversity=true, the strong diversification scheme is applied in order to escape from the current solution. In strong diversification, is set to , , the last local branching constraint ( ̅ ) is replaced with ( ̅ ) , , the local branching constraint ( ̅) is added, and the process is iterated. Note that is to only count the number of times strong diversification comes to effect.

III.

A solution with objective function strictly better than lower bound LB is found, but the optimality for the current branch is not proven. In this case, the last local branching constraint ( ̅) is replaced with constraint ( ̅ ) , unless flag first true, which means that constraint ( ̅ ) has already been introduced at Step 2, in which case, only the last local branching constraint is removed. Subsequently, lower bound , best solution , best NPV, best lower bound , and ̅ are updated, the local branching constraint ( ̅ ) is added, and the process is iterated.

IV.

The current left branch is proven to have no feasible solution of objective function strictly better than lower bound . In this case, the last local branching constraint is reversed into ( ̅) . If diversity false, the last local branching constraint is removed, rhs rhs+⌈ ⌉, the local branching constraint ( ̅ ) is added, and the process is iterated. If

AC

CE

PT

ED

M

AN US

CR IP T

I.

diversity true, is set to diversity false is followed.

,

, first

true, and the same process as when

We stop local branching and call our adaptive branching scheme when a predefined limit is reached, that is, when . In addition, we set a time limit total_time_limit for the run of the algorithm.

10

ACCEPTED MANUSCRIPT

until (

CE

PT

4.

ED

M

3.

diversify true; endif; if solution ̅ found, optimality not proven then if (first) then delete the last local branching constraint ( ̅ ) ; else replace the last local branching constraint ( ̅ ) with ( ̅) endif; if ( ( ̃) ) then ( ̃) ̃ endif; first diversify false, ̅ ̃ ( ̃) endif; if infeasibility proven then reverse the last local branching constraint ( ̅ ) into ( ̅ ) ; if (diversify) then first true endif; ⌈ ⌉ diversify true; endif;

AN US

2.

reverse the last local branching constraint into ( ̅ ) , diversify first false, ̅ ̃ ( ̃); endif; if no feasible solution found, infeasibility not proven then if (diversify) then replace the last local branching constraint ( ̅ ) with ( ̅) , , ⌈ ⌉ first true; else delete the last local branching constraint ( ̅ ) , ⌈ ⌉ endif;

CR IP T

Algorithm 2: Local branching and adaptive branching scheme. Inputs: feasible solution ̃, objective value ( ̃), dv_max, total_time_limit, ( ̃) ( ̃) node_time_limit, ̃ , diversity false, first true, ̅ ̃. begin repeat Add the local branching constraint ( ̅ ) , and solve the created subregion by the branching scheme within node_time_limit; 1. if optimality proven then ( ̃); endif; if ( ( ̃) ) then ̃

)

AC

solve the problem with added local branching constraints using the branching scheme in the remaining time. return ( ) end

4.4. Starting integer feasible solution Local branching requires a starting feasible solution. In the literature, this solution is generated by either the black-box solver of local branching (e.g., CPLEX) or a heuristic. In this section, we develop a new heuristic to quickly produce a starting feasible solution for local branching. Our heuristic has been built based on topological sort, an idea introduced by Cheng and Gen (1998). Algorithms I and II summarize the procedure in the heuristic. Suppose that ( ) is a directed graph 11

ACCEPTED MANUSCRIPT

ED

M

AN US

Algorithm I Inputs: priority array generated with uniform distribution ( ) directed graph ( ), , , Indegree[i] , top[i] , cut[i] , cut[i] , . begin for to do for to do if (cut[i] ) then * + endif; find the highest priority node from ; top[k] ; cut[ ] ; ; (the number of immediate successors of node ) for to do Indegree[ , -] Indegree [ , -] ; , -; if (Indegree[ , -]= ) then cut[ , -] endif; end

CR IP T

representing the immediate precedence relationships between blocks in a block model, with a set of nodes indicating the blocks and a set of arcs S. For an arc * + , we say that it is linked from node i to node j. Let cut be a array, where cut[i] shows the element. Let and be the set of eligible nodes and a priority array, respectively. top denotes the set of nodes removed from the priority array and added to the topological sort, so top[k] means that node is placed in the element of array top. Indegree[i] is an array that controls the eligible nodes. denotes the set of immediate predecessors of node , so Indegree[i] means that the number of edges linked to node i equals the number of its immediate predecessors. Also, , - indicates the element of array , where is the set of immediate successors of node . For example, {3, 5, 8} shows that nodes 3, 5, and 8 are immediate successors of node i, and so [2] calls the second element of this array, which is node 5. Finally, let be the set of blocks on the surface of the block model.

AC

CE

PT

In Algorithm II, top is an array containing the topological sort created by Algorithm I. For example, top [3]=5 shows that the third element of the topological sort is node 5. Array t[i] represents the time period in which node i has been scheduled. t[2]= 3, for instance, states that node 2 has been scheduled in time period 3. ( ) and ( ) denote the tonnage of node i and its amount of ore, respectively. Note that every constraint except the lower bounds is guaranteed to be satisfied in the first loop of Algorithm II. The feasibility of the lower bounds at the end of each loop is checked. If they are not satisfied in the current schedule, a new topological sort is generated by changing the priority array, and a new schedule is produced accordingly; otherwise Algorithm II returns the schedule as a feasible solution. Algorithm II repeats until the generated schedule meets the lower bounds. Algorithm II repeat until (schedule t satisfies the lower bounds) begin generate a random priority array; generate topological sort “top” by Algorithm 1; t[i] ; top[1]; t[ ] ; ( ); ( );

12

ACCEPTED MANUSCRIPT for

to do top[i]; if ( ) then ; else t[ ] ( is the node with the highest time period in the set endif; repeat until ( ( ) or ( ) ) ; [end of loop] t[ ] ; ( ); ( );

CR IP T

end [end of loop]

)

Therefore, when solving an instance, we first run Algorithm II to generate a starting feasible solution. This solution is then transferred to the methodology developed in Section 4.3, i.e., Algorithm 2. The final solution is returned by this algorithm when its stopping criterion is reached.

AN US

5. Experimental results

For the purpose of testing and validating the newly created solution procedure, we solve a series of conceptually generated data sets. We compare the results to those obtained by the B&C and Lagrangian relaxation techniques proposed in the literature to address CPIT. We also run the same data sets in CPLEX 12.6. The motivation of this is to show that general IP solvers are not able to achieve near-optimal solutions for our data sets.

M

5.1. Data

CE

PT

ED

We create a series of six conceptual block models as outlined in Table 1. While the data for this is conceptual in nature, it is, however, based on real operational scenarios with block size, tonnages, grades, resource limitations and sequencing interactions reflective of real open pit mining operations. The setting of this mine is a typical remote mining region within Australia. The deposit itself has been modelled as a typical porphyry style copper (Cu) deposit. Porphyry deposits are characterised by a steeply dipping (almost vertical) centre core of high grade mineralisation. The core is surrounded by a halo of material which progressively reduces in grade the further away it is located from the centre core. Table 1 Block model data. Block x y model

Max Min processing processing capacity (t) capacity (t) 37 37 11 15,059 30,000,000 6,000,000 4,200,000 3,600,000 1 46 46 12 25,392 37,500,000 6,750,000 4,500,000 3,900,000 2 48 48 13 29,952 56,250,000 7,500,000 5,625,000 4,800,000 3 50 50 14 35,000 90,000,000 8,250,000 5,400,000 4,800,000 4 50 50 16 40,000 105,000,000 9,000,000 8,250,000 5,400,000 5 53 53 18 50,562 127,500,000 9,000,000 9,750,000 5,400,000 6 The first column represents each of the six block models. The next three columns detail the number of blocks in the „x‟, „y‟ and „z‟ direction contained within each block model. Multiplying the number of blocks in each direct provides the total number of blocks as contained in the fifth column. Columns six and seven provide the maximum and minimum tonnage of material to be mined in each year. The final two columns provide the maximum and minimum tonnage of material to be processed in each year.

AC

z

Number of blocks

Max mining capacity (t)

13

Min mining capacity (t)

ACCEPTED MANUSCRIPT

Table 2 Input parameters. Mining cost ($/t):

2.40

Processing cost ($/t):

9.60

Price ($/t Cu):

5000

Recovery (%):

0.90

Discount Rate (%):

0.10

Cut-off grade (% Cu):

0.21

CR IP T

Each block is 30 metres long (x), 30 metres wide (y) and 30 metres high (z), containing a total of 75,000 tonnes of material at a density of 2.78m3/t, regardless of whether it contains copper mineralisation or no mineralisation at all. Associated with this data are maximum and minimum mining and processing capacities appropriate for this size and scale of deposit.

Based on Table 2, the value of blocks containing no grade (waste) is expressed as: 2.4

180,000

(9)

AN US

75,000

And the value of blocks containing grade (potential ore) is expressed as: (

)

)

(10)

(

)

(

For example, for a block containing grade of 0.89% Cu: 0.0089 0.9

5,000

)

(75,000

M

(75,000

2.4 )

(75,000

9.6 )

2,103,750

ED

A final pit wall angle of 45 degrees is assumed in every block model. This requires a total of five overlying blocks to each block to be removed in full (directly above, above north, above east, above south, and above west).

PT

5.2. Data pre-processing

CE

Before implementing our methodology, it is necessary to examine the data sets to remove redundant data. In practice, using data blindly can result in extra computations.

AC

The first data set with 15,059 blocks, for example, is 37 blocks wide in the x-coordinate direction, 37 blocks long in the y-coordinate direction, and 11 blocks deep in the z-coordinate direction (from the bottom up). The cut-off grade is fixed to 0.21%, and any block with a grade of less than that is recognised as waste. The first step of pre-processing is performed by eliminating two groups of blocks. Those two groups are: 1) waste blocks that are located at the bottom level of the block model, and 2) waste blocks that have no ore blocks in their entire set of successors (i.e, set ). Table 3 shows the number of remaining blocks in each block model after the first step of pre-processing is completed.

Table 3 Number of blocks after the first step of pre-processing. Block model Number of blocks 8,245 1

14

ACCEPTED MANUSCRIPT 10,323 15,112 25,160 28,931 34,820

2 3 4 5 6

AN US

CR IP T

Gaupp (2008) develops a technique called Variable Reduction (VR) in which the blocks that cannot be extracted in a particular time period, due to its predecessors and/or successors, are identified, and the corresponding decision variables are removed (see Gaupp (2008) for a full description). In the second step of data pre-processing, we generate the sets of predecessors and successors for each block in order to apply the VR technique. For clarity, we refer to each block model by its number of blocks, as shown in the second column of Table 4. For example, we refer to block model 1 as problem P8,245. The characteristics of each problem are summarized in the same table. Column T shows the number of time periods. Each instance involves 10 time periods (years). The number of variables before and after applying VR appear in columns Variables and Variables (ES+LS), respectively. Also, the total number of constraints before and after VR is represented in columns Constraints and Constraints (ES+LS). Table 4 Problems‟ characteristics. T

Variables

1

P8,245

10

82,450

2

P10,323

10

3

P15,112

4

Constraints 369,452

Variables (ES+LS) 65,905

Constraints (ES+LS) 297,606

103,230

470,603

88,390

401,045

10

151,120

695,811

119,449

548,251

P25,160

10

251,600

643,109

214,202

971,542

5

P28,931

10

289,310

1,321,743

225,603

1,029,067

6

P34,820

10

1,608,172

236,200

1,090,672

M

Problem

ED

Block model

348,200

PT

5.3. Numerical results

AC

CE

We test the performance of our methodology on the problems in Table 4. For all the computations, we use the data generated after applying the VR procedure. The entire code is written in Java and run on a personal computer with Intel Core i7 CPU at 3.40 GHz, 16 Gigabytes of RAM, operating under Microsoft Windows 7. The MIP solver in CPLEX 12.6 is called through ILOG Concert Technology during B&C. The LP relaxation‟s optimal value for all instances is computed in IBM ILOG CPLEX studio 12.6. The solution quality is given with respect to the optimal solution value from the LP relaxation of CPIT, computed as follows: optimality gap = where is the optimal LP relaxation value of the corresponding problem, and solution found by a given algorithm.

represents the

Initial integer feasible solutions are generated by our heuristic and then given to local branching. The time of this is included in the total solution time. 15

ACCEPTED MANUSCRIPT After testing several other options, we set the algorithm‟s parameters to the following values:    

k:=0.8% of the number of variables node_time_limit:=5% of the total time limit :=30% dv_max= unlimited

Thence, the regression equation is: total_time_limit

0.0024 number_of_variables 8.7

(11)

CE

PT

ED

M

AN US

where total_time_limit is produced in minutes.

CR IP T

The second column of Table 6 gives the best total_time_limit for each instance, obtained after conducting some preliminary experiments. We choose the best total_time_limit for each instance in such a way as to reach a good compromise between running time and solution quality. Fig. 1 shows total_time_limit for each instance. To help readers approximate the best total time limits required for different sizes (i.e., different numbers of variables), a regression line is fitted (see Montgomery et al. (2012)).

AC

Fig. 1. Number of variables versus the total time limit. The solid line represents the linear regression fitting of the data.

As stated earlier, our algorithm stops local branching when the number of strong diversifications, i.e., dv, reaches a predefined limit dv_max, and continues with the adaptive branching scheme for the remaining time until the total time limit total_time_limit is reached. In this case, the strong diversification scheme may come into effect several times. However, the concept of strong diversification is associated with finding a new reference solution that is often worse than the current one, which means that no improvement in the solution may be found in spite of the time it uses at this step. Thus, with dv_max=unlimited, the algorithm may waste a significant amount of the total time limit total_time_limit by dedicating it to the implementation of strong diversification. On the other hand, removing this step entirely from the procedure may cause the algorithm to get stuck during its 16

ACCEPTED MANUSCRIPT early stages. Therefore, finding a reasonable value for the maximum number of strong diversifications (i.e., dv_max) is of importance.

dv

Gap % Gap % Gap % Gap % (dv_max=unlimited) (dv_max=0) (dv_max=5) (dv_max=10) P8,245 7 2.3 1.1 1.8 2.3* P10,323 15 4.2 2.4 0.8 2.2 P15,112 9 2.4 2.0 1.6 2.4* P25,160 11 2.7 4.8 1.9 2.6 P28,931 2 1.8 2.3 1.8* 1.8* P34,820 10 2.1 4.1 2.4 2.1* Average ……… 2.58 2.78 1.72 2.23 *The gap here is the same as that of the solution obtained with dv_max=unlimited. This is because the number of times that the strong diversification scheme has come into effect, while solving the instance with dv_max=unlimited, is less than the current value of the parameter dv_max.

PT

ED

M

Problem

AN US

Table 5 Optimality gaps obtained by different values of dv_max.

CR IP T

To analyse the effect of dv_max and to obtain the best value for this parameter, we solve the problems with dv_max= unlimited, and record the number of strong diversifications performed. These numbers are presented in the second column of Table 5. Subsequently, we conduct different experiments by setting dv_max to 0, 5, and 10. All other parameters are set as above. The detailed results are shown in Table 5. In this table, column Gap (dv_max=unlimited) represents the results of solving the problems when parameter dv_max is set to unlimited. Similarly, the fourth, fifth and last coulmns show the results when this parameter is set to 0, 5, and 10, respectively. The algorithm yields optimality gaps of 2.3%, 1.1%, 1.8%, and 2.3%, for P8,245 when parameter dv_max is set to unlimited, 0, 5, and 10, respectively. This means that the algorithm works most effectively with dv_max=0 for this instance. The bottom of the table (Average) shows the average optimality gaps obtained for each value of parameter dv_max. As can be seen, the minimum average optimality gap is obtained with dv_max=5, and thus, we consider this value to be the best for dv_max.

AC

CE

Column Local branching in Table 6 displays the optimality gaps found by our algorithm with dv_max=5 within the given total time limit. We see that our algorithm obtains a solution within less than 2.5% of optimality for every instance. We have also included the results of running model on CPLEX 12.6 with default settings. Column CPLEX in Table 6 shows the optimality gaps obtained by CPLEX with the same time limit as in our algorithm. From the table, CPLEX is only able to find an integer solution for two of our six problems. For P8,245 and P10,323, CPLEX manages to establish optimality gaps of 12% and 5.3%, respectively. This shows that our algorithm obtains strictly better solutions for every instance. Due to the fact that the local branching procedure is based on generating cuts, and also since it uses B&B–in some cases B&C–as the black box, it is an interesting evaluation for our algorithm to be compared to the developed B&C techniques in the literature. To do so, we solve our instances one more time using Gaupp‟s B&C technique (Gaupp, 2008), with CPLEX parameters set in the same way as in this study. Gaupp terminates B&C when the optimality gap drops to 2% or less; however, we terminate it when the total time limit used in our algorithm is reached. Table 6 shows that local branching produces results of better quality than Gaupp‟s B&C. For P10,323, B&C reaches an optimality gap of 2.8% in 4h: 12min, which is the minimum gap among all the obtained results, while 17

ACCEPTED MANUSCRIPT our algorithm reaches 0.8% of optimality. B&C cannot find any integer feasible solution for P34,820 within 10h: 16min. The gap stays above 5% for the other four instances.

Local branching 1.8% 0.8% 1.6% 1.9% 1.8% 2.4%

CPLEX

B&C

12.0% 5.3% * * * *

5.4% 2.8% 6.2% 10.7% 8.6% *

Lagrangian relaxation 1.1% 0.8% ** 4.8% * * *

CR IP T

Table 6 Numerical results. Problem Total time limit P8,245 2h: 43min P10,323 4h: 12min P15,112 4h: 41min P25,160 8h: 06min P28,931 9h: 18 min P34,820 10h: 16 min

*No feasible solution is found within the given time limit. ** The optimality gap of 0.8% is obtained before the given time limit is reached, i.e., within 1h: 11min, but it does not improve in the remaining time.

ED

M

AN US

Lagrangian relaxation is another solution technique that has been employed in the literature to address OPMPSP when considering minimum resource requirements. Therefore, we also compare our results to those obtained by Gaupp‟s Lagrangian relaxation technique (Gaupp, 2008). This technique uses the subgradient method for multiplier updates and employs a feasing routine to transform Lagrangian solutions from infeasibility into feasibility. The method terminates when the optimality gap drops to 2% or less. In our comparison, however, we terminate it once our total time limit is reached. The results appear in column Lagrangian relaxation in Table 6. From the table, Lagrangian relaxation shows a better performance than CPLEX, B&C, and our methodology for small size instances. It yields an optimality gap of 1.1% for P8,245, while our methodology, B&C, and CPLEX give 1.8%, 5.4%, and 12%, respectively, for the same instance. Also, for P10,323, Lagrangian relaxation reaches an optimality gap of 0.8% within only 1h: 11min, which is much faster than our methodology–it could not improve it in the remaining time. However, this method does not yield very good results for larger instances. In fact, it cannot find any feasible solution within the given time limit for instances larger than P15,112.

AC

CE

PT

To further analyse our algorithm as well as the efficiency of the tuned parameters, i.e., dv_max, node_time_limit, , k, and total_time_limit, we aim to solve two larger instances generated in a similar way to previous instances. As in Table 4, Table 7 shows the information about these two instances. The number of blocks in the block models–before pre-processing–are 56,180 and 74,008 for P39,625 and P50,383, respectively. Table 8 indicates the solution results. Using the regression equation, total time limits of 11h: 51min and 13h: 23min are selected for P39,625 and P50,383, respectively. For both of these two large instances, the gap stays below 2.5%, demonstrating that the tuned parameters work properly for other sizes. With identical time limits, CPLEX and Lagrangian relaxation are unable to find a feasible solution for either of the instances, and B&C can only obtain an optimality gap of 11.2% for P39,625. Table 7 Characteristics of two larger instances. Problem Number of T blocks P39,625 39,625 10 P50,383 50,383 10

Variables

Constraints

396,250 503,830

1,825,731 2,283,452

18

Variables (ES+LS) 292,746 331,052

Constraints (ES+LS) 1,349,706 1,526,913

ACCEPTED MANUSCRIPT Table 8 Numerical results of two larger instances. Problem Total time Local limit branching

CPLEX

P39,625 11h: 51min 2.4% * P50,383 13h: 23min 2.2% * *No feasible solution is found within the given time limit.

B&C

Lagrangian relaxation

11.2% *

* *

ED

M

AN US

CR IP T

In this study, we have specified a size-based termination time for our algorithm, which may make it difficult for future studies to compare their solutions to ours. Indeed, due to the lack of appropriate series of benchmarks in this research area, it is difficult, if not impossible, to make a time comparison. To alleviate this situation, we let our algorithm run for 20 hours for every problem and record the behaviour. Fig. 2 shows the optimality gaps found in each hour. The gaps at time zero indicate the initial feasible solution obtained by our heuristic.

Conclusion

CE

6.

PT

Fig. 2. The optimality gaps obtained in each hour over a running time of 20 hours. The gap axis is on a log scale.

AC

The open pit mine production scheduling problem is an NP-hard combinatorial problem that has attracted the attention of many researchers since the 1960s. The major difficulty of this scheduling problem is its natural complexity that occurs due to its operational requirements and large size. The main operational requirements that should be considered when scheduling production are: 1) precedence relationships between blocks rising due to the wall slope restriction; 2) mining and processing capacities; and 3) mining and processing requirements. Due to the fact that real-world problems are extremely large in terms of capacity constraints, destinations, and number of blocks, a simplified version of OPMPSP is generally considered. A variety of simplifying assumptions can be made, such as no-inventory, fixed cut-off grade, and knapsack constraints with lower bounds of zero. This paper has adopted a local branching method to tackle OPMPSP while including lower bounds for knapsack constraints. Local branching is a novel idea whereby mathematical programming and metaheuristic techniques are combined. To use local branching in a more efficient way for our problem, we have developed a new 19

ACCEPTED MANUSCRIPT adaptive branching scheme and combined it with this technique. We have also proposed a heuristic to generate starting feasible solutions to accelerate our methodology. We have tested our methodology on a set of conceptually generated data sets. The outcomes demonstrate that our algorithm outperforms both the B&C and Lagrangian relaxation techniques that are applied to OPMPSP in the literature. Also, in comparison to other similar studies that take into consideration the lower-bounding constraints, we have addressed larger problems.

CR IP T

As future work, we intend to further improve local branching in order to tackle even larger instances. We also aim to apply local branching to OPMPSP with dynamic cut-off grades and stockpiles, which are more complex versions of the problem.

References

AC

CE

PT

ED

M

AN US

Beale, E., & Forrest, J. (1976). Global optimization using special ordered sets. Mathematical programming, 10(1), 52-69. Bienstock, D., & Zuckerberg, M. (2010). Solving LP relaxations of large-scale precedence constrained problems. Integer Programming and Combinatorial Optimization (pp. 1-14): Springer. Boland, N., Dumitrescu, I., Froyland, G., & Gleixner, A. M. (2009). LP-based disaggregation approaches to solving the open pit mining production scheduling problem with block processing selectivity. Computers & Operations Research, 36(4), 1064-1089. Caccetta, L., & Hill, S. P. (2003). An application of branch and cut to open pit mine scheduling. Journal of global optimization, 27(2-3), 349-365. Cheng, R., & Gen, M. (1998). An evolution programme for the resource-constrained project scheduling problem. International Journal of Computer Integrated Manufacturing, 11(3), 274-287. Chicoisne, R., Espinoza, D., Goycoolea, M., Moreno, E., & Rubio, E. (2012). A new algorithm for the open-pit mine production scheduling problem. Operations Research, 60(3), 517-528. Cullenbine, C., Wood, R. K., & Newman, A. (2011). A sliding time window heuristic for open pit mine block sequencing. Optimization letters, 5(3), 365-377. Espinoza, D., Goycoolea, M., Moreno, E., & Newman, A. (2013). MineLib: a library of open pit mining problems. Annals of Operations Research, 206(1), 93-114. Fischetti, M., & Lodi, A. (2003). Local branching. Mathematical programming, 98(1-3), 23-47. Fischetti, M., Polo, C., & Scantamburlo, M. (2004). A local branching heuristic for mixed‐integer programs with 2‐level variables, with an application to a telecommunication network design problem. Networks, 44(2), 61-72. Gaupp, M. P. (2008). Methods for improving the tractability of the block sequencing problem for open pit mining: DTIC Document. Hansen, P., Mladenović, N., & Urošević, D. (2006). Variable neighborhood search and local branching. Computers & Operations Research, 33(10), 3034-3045. Jélvez, E., Morales, N., Nancel-Penard, P., Peypouquet, J., & Reyes, P. (2016). Aggregation heuristic for the open-pit block scheduling problem. European Journal of Operational Research, 249(3), 1169-1177. Kawahata, K. (2007). A new algorithm to solve large scale mine production scheduling problems by using the Lagrangian relaxation method. Colorado School of Mines.

20

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

Kumral, M. (2012). Production planning of mines: Optimisation of block sequencing and destination. International Journal of Mining, Reclamation and Environment, 26(2), 93-103. Lambert, W., & Newman, A. (2014). Tailored Lagrangian Relaxation for the open pit block sequencing problem. Annals of Operations Research, 222(1), 419-438. Liu, S. Q., & Kozan, E. (2016). New graph-based algorithms to efficiently solve large scale open pit mining optimisation problems. Expert Systems with Applications, 43, 59-65. Montgomery, D. C., Peck, E. A., & Vining, G. G. (2012). Introduction to linear regression analysis (Vol. 821): John Wiley & Sons. Moreno, E., Espinoza, D., & Goycoolea, M. (2010). Large-scale multi-period precedence constrained knapsack problem: a mining application. Electronic notes in discrete mathematics, 36, 407414. Ramazan, S. (2007). The new fundamental tree algorithm for production scheduling of open pit mines. European Journal of Operational Research, 177(2), 1153-1166. Rodríguez-Martín, I., & Salazar-González, J. J. (2010). A local branching heuristic for the capacitated fixed-charge network design problem. Computers & Operations Research, 37(3), 575-581. Shishvan, M. S., & Sattarvand, J. (2015). Long term production planning of open pit mines by ant colony optimization. European Journal of Operational Research, 240(3), 825-836. Underwood, R., & Tolwinski, B. (1998). A mathematical programming viewpoint for solving the ultimate pit problem. European Journal of Operational Research, 107(1), 96-107. Zhu, G., Bard, J. F., & Yu, G. (2006). A branch-and-cut procedure for the multimode resourceconstrained project-scheduling problem. INFORMS Journal on Computing, 18(3), 377-390.

21