European Journal of Operational Research 257 (2017) 297–309

Contents lists available at ScienceDirect

European Journal of Operational Research

journal homepage: www.elsevier.com/locate/ejor

Innovative Applications of O.R.

Optimizing the open pit-to-underground mining transition

Barry King a, Marcos Goycoolea b, A. Newman c,∗
a
Colorado School of Mines, Operations Research with Engineering (ORwE), Golden, CO 80401, USA
b
School of Business, Universidad Adolfo Ibañez, Peñalolén, Santiago 7941169, Chile
c
Colorado School of Mines, Mechanical Engineering Department, Golden, CO 80401, USA

a r t i c l e i n f o a b s t r a c t

Article history: A large number of metal deposits are initially extracted via surface methods, but then transition under-
Received 24 August 2015 ground without necessarily ceasing to operate above ground. Currently, most mine operators schedule
Accepted 11 July 2016
the open pit and underground operations independently and then merge the two, creating a myopic so-
Available online 16 July 2016
lution. We present a methodology to maximize the NPV for an entire metal deposit by determining the
Keywords: spatial expanse and production quantities of both the open pit and underground mines while adhering
Mining/metals industries: optimal to operational production and processing constraints. By taking advantage of a new linear programming
extraction sequence solution algorithm and using an ad-hoc branch-and-bound scheme, we solve real-world scenarios of our
Production scheduling: transition problem transition model to near optimality in a few hours, where such scenarios were otherwise completely in-
Integer programming applications: exact tractable. The decision of where and when to transition changes the net present value of the mine by
and heuristic approaches hundreds of millions of dollars.
© 2016 Elsevier B.V. All rights reserved.

1. Introduction and literature review extraction method that results in the largest undiscounted profit
for each three-dimensional discretization of the ore body and sur-
The mining industry contributes trillions of dollars annually to rounding rock. Mine operators tend to delay the transition, leading
the global economy by providing minerals, metals, and aggregates. to NPV losses of up to hundreds of millions of dollars. We provide
This, and volatile metal prices, make it critical that mines pos- a systematic means by which a mine operator can determine the
sess an efficient production schedule, which can be categorized highest value of a combined open pit and underground design.
as: (i) short-term (days to months), (ii) long-term (years), and The most common method used to extract material is open pit,
(iii) strategic (life-of-mine) (Gershon, 1983). A short-term sched- or surface, mining. Open pit mines vary in both shape and size,
ule might determine what material to process on a given day; a and their design is based on the deposit’s block model, a model
long-term schedule may examine production rate changes (Alonso- which discretizes the orebody and surrounding rock, and assigns a
Ayuso et al., 2014; Epstein et al., 2014). Finally, a strategic schedule series of attributes, including mining cost, degree of mineralization
is used to evaluate large capital investments, and other decisions (referred to as grade), location, and the cost or profit associated
that have long-ranging impacts. Because the transition from open with processing the specific block. Blocks can be categorized us-
pit to underground extraction affects a mine for the remainder of ing a minimum cutoff grade; blocks at or above the cutoff grade
its operational life, it falls into the category of strategic scheduling. are sent to the processing plant, referred to as a mill, while those
At the time of this writing, a large number of metal deposits below the cutoff grade are sent to a waste dump. The slope angle
are being extracted via surface methods, but plan to transition to for the open pit mine, resulting from geotechnical constraints of
concurrently or exclusively extracting ore via underground mining the host rock, ensures the stability of the pit’s walls (Crawford &
methods. For safety reasons, the underground mine must be suf- Hustrulid, 1979).
ficiently geographically separated, with horizontally positioned in Given the block attributes and slope angles, mine planners de-
situ rock, from the open pit mine via what is typically referred to termine the largest economically viable pit for a given deposit,
as a crown pillar. Current industry practice places the crown pillar i.e., the ultimate pit limit (Lerchs, 1965; Underwood & Tolwinski,
based on: (i) largest economically viable open pit mine, or (ii) the 1998). However, while the solution to the ultimate pit limit prob-
lem yields the size of the open pit mine, it provides no indication
∗
of the extraction sequence required to maximize its discounted
Corresponding author.
value. Johnson (1968) originally formulated the open pit block
E-mail addresses: [email protected] (M. Goycoolea), [email protected]
(A. Newman). sequencing problem as an integer program that schedules the

http://dx.doi.org/10.1016/j.ejor.2016.07.021
0377-2217/© 2016 Elsevier B.V. All rights reserved.
298 B. King et al. / European Journal of Operational Research 257 (2017) 297–309

Fig. 1. (Open stoping) In this mining method, rib pillars provide stability, as does
the backfilling of open voids left by extracted stopes. Stope advance shows the di-
rection in which mining proceeds.
Fig. 2. (Transition zone) The transition zone is an area where it is economically
extraction of blocks such that the open pit’s value is maximized viable to extract material via open pit or underground mining methods. We see the
open pit, black, encroaching on the underground mine, gray, in the transition zone.
subject to resource and precedence constraints.
Solution techniques for open pit block sequencing problems are
still widely studied (Chicoisne, Espinoza, Goycoolea, Moreno, & Ru- The latter characteristic implies that our model contains no contin-
bio, 2012; Osanlo, Gholamnejed, & Karimi, 2008; Ramazan, 2007; uous variables. On the other hand, we determine sill pillar place-
Shishvan & Sattarvand, 2015; Souza, Coelho, Ribas, Santos, & Mer- ment, i.e., locations in which material is left in situ to allow for a
schmann, 2010; Topal & Ramazan, 2010). One such recent signif- change in mining direction, which adds a layer of complexity.
icant advance for the linear programming relaxation of a gen- An early transition model assigns large aggregated blocks to be
eral version of the so-called precedence constrained production extracted via open pit or underground mining methods in order to
scheduling problem (PCPSP), i.e., the open pit block sequencing maximize the value of the deposit (Bakhtavar, Shahriar, & Oraee,
problem, is with the use of an algorithm outlined in Bienstock and 2008). This idea was later improved to include the time element
Zuckerberg (2010), which exploits the problem structure (Muñoz and to capture underground capital costs (Newman, Yano, & Rubio,
et al., 2015). Lambert, Brickey, Newman, and Eurek (2014) present 2013). In both previous transition models, there is little differenti-
a guide to formulating and efficiently solving monolithic instances ation between the mining units used above and below ground. The
of the open pit block sequencing problem, i.e., without decompo- mining industry comments on the difficulty of modeling the transi-
sition. tion correctly (Finch, 2012); however, decisions regarding the tran-
Underground mining is used when an economically viable de- sition are becoming increasingly relevant (Araneda, 2015). Fig. 2
posit is situated sufficiently deep such that open pit mining is shows an open pit atop an underground mine. The transition zone
cost prohibitive. There exist many underground mining techniques: is depicted as the material that would be extracted were it done
(i) open stoping (Fig. 1), (ii) room-and-pillar, (iii) sublevel caving, via underground methods; the corresponding amount of material
(iv) drift-and-fill, (v) longwall, and (vi) block caving. Determining would be greater were open pit methods used in the transition
which method(s) to use is typically based on geotechnical con- zone.
straints, size, and shape of the deposit (Qinglin, Stillborg, & Li, We present a new model and corresponding solution tech-
1996). For the purpose of this paper, we confine our discussion to niques to determine the timing of a transition from open pit
open stoping mining and its associated sequencing options. to underground mining in both a spatial and a temporal sense.
A stope is a large, three-dimensional, mineable volume whose This transition incorporates a crown pillar placement that sepa-
maximum size is correlated with the geotechnical properties of the rates the open pit from the underground mine, and of the sill pil-
host rock, and is the basic unit for stoping methods. The void left lars, i.e., levels left in situ that can grant earlier access to stopes
by an extracted stope is sometimes filled with an aggregate to pro- by creating a false bottom. Our methodology is based on an ad-
vide structural stability, a process referred to as backfilling. Most hoc branch-and-bound approach that incorporates decomposition
underground stoping mines are separated into vertically spaced methods for solving PCPSP linear programming relaxations, and
levels based on the maximum stope height, creating a near-regular that includes rounding heuristics. We outline underlying models
grid of possible stope positions (Alford, 2006). for the transition in Section 2. Mathematical reformulations to
After determining possible locations from which the ore can enhance tractability are presented in Section 3, and the solution
economically be extracted, i.e., possible stope locations, mine plan- strategy in Section 4. Sections 5 and 6 provide the numerical re-
ners design the development (Alford, 2007; Brazil & Thomas, sults and conclusions, respectively.
2007), which is required to gain access to the ore, provide haulage
routes, and maintain proper ventilation within the underground 2. Underlying models
mine. All stoping activities require the completion of a specific set
of development activities before that stope’s extraction can com- In this section, we introduce three models that underlie our
mence. Underground sequencing constraints are created after the computationally tractable transition model. We first present a sur-
design, and provide rules for the order in which the development face extraction formulation, followed by an underground formula-
and stopes are extracted. Given a fixed design and sequencing tion, and conclude with a preliminary transition formulation which
method, we can schedule the underground mining activities to, is essentially a combination of the two.
e.g., maximize NPV, or minimize deviation from production tar-
gets (Brickey, 2015; Carlyle & Eaves, 2001; Martinez & Newman, 2.1. Surface model
2011; Newman & Kuchta, 2007; O’Sullivan & Newman, 2015). Trout
(1997) provides one of the first generalized formulations for un- We consider a surface model based on open pit mining with a
derground stope scheduling; our formulation is a bit more stream- multi-phase pit design (Fig. 3), in which a phase corresponds to
lined than his in that we do not differentiate between scheduled a sub-region of the pit. A block within a phase consists of all of
and actual decisions, and because we assume that once an activity the material in the phase that resides within a predefined verti-
commences, it must continue at a prescribed rate until finished. cal distance. (Note that some mine operators refer to our blocks
 B. King et al. / European Journal of Operational Research 257 (2017) 297–309 299

Fig. 3. (Phase-block-bin data aggregation) The phases are shown as sub-pits with a block occupying a small vertical space within the phase. Each block is separated into bins
based on grade. Material in the low- (stripes), medium- (checkered), and high- (waves) grade bins may go directly to the mill (dashed arrows), or to an individual stockpile
(solid arrows). Waste is sent to the dump. Although naturally occurring material of different grades is scattered within the block, for stylistic purposes, we group material of
each grade.

as benches.) Inside each block, there exists a series of bins that cial solution strategy, we omit rehandling costs from our formu-
are differentiated by grade, categorized as waste, low, medium, lations. For our transition model, these costs prove to be insignifi-
or high, and geological properties. This type of phase-block-bin cant; post-processing them into the objective function value results
scheduling is common in the mining industry and is the ba- in hundredths of a percentage change, an amount that is not likely
sis for the Minemax (2013) software package. Whittle Consulting to significantly increase were we to impose the rehandling costs a
(2013) have developed multiple products for this type of schedul- priori and certainly not sufficiently substantial to consider as part
ing. of strategic planning costs.
The objective of the surface model, (S), is to schedule the ex- We define the notation below. In general, use of lower case let-
traction, stockpiling, and mill feed in such a way that the NPV is ters is reserved for indices and parameters. Upper case letters in
maximized, while adhering to annual extraction and milling ca- roman font represent variables, and sets are given in calligraphic
pacity constraints. In addition, the desired shape of the open pit font. An S superscript on a parameter or variable denotes notation
is maintained by precedence constraints, which can be catego- specific to the surface model; we use hats to differentiate parame-
rized into two types: (i) intra-phase precedence expresses that the ters and variables that represent similar entities.
blocks inside the phase be extracted from the surface down and
(ii) inter-phase precedence expresses that blocks inside a phase be Indices and sets:
extracted only after a specific block in the predecessor phase has b∈B blocks b
been fully extracted. A maximum sinking rate restricts the number n ∈ Nb bins in block b
of blocks in each phase one can mine in a given time period based bˆ ∈ Bˆb blocks that must be mined directly before block b
on operating constraints. We also require that the material con- b ∈ Bp blocks in phase p
tained in each block and bin be extracted in equal proportion to d∈D bin destination (1 = mill, 2 = stockpile, 3 = waste )
prevent selective extraction within the block. Once extracted, in- p∈P phases p
dividual bins can be directed to one of three destinations: waste r∈R resources (1 = mine, 2 = mill, 3 = sinking rate )
dump, stockpile, or mill. All material that is below cutoff grade is t∈T time periods
sent to the waste dump. Data:
Stockpiling is important in our application, because processing S−
cnb mining cost for bin n in block b [$]
stockpiled material augments the underground production to en- S+
cnb revenue generated after having milled bin n of block
sure that the mill remains at maximum capacity after extraction
b [$]
ceases in the open pit mine. Bley, Boland, Froyland, and Zucker-
qSrnb quantity of resource r consumed by bin n of block
berg (2012) outline the formulation from which we construct our
b [1 and 2 = tonnes]
“warehouse-style” stockpiling strategy, i.e., each stockpile contains
r rt , r̄rt minimum, maximum amount of resource r available
only one block-bin combination; the objective function value cor-
in time t [1 and 2 = tonnes, 3 = blocks]
responding to an optimal solution subscribing to this strategy pro-
δt discount factor for time period t [fraction]
vides an upper bound on the NPV that can be obtained. Mate-
rial retrieved from the stockpile is identical to material placed in Decision variables:
the stockpile. Some authors have attempted to use “mixing con- XbtS 1 if block b has finished being extracted by the end of
straints” to more accurately model the characteristics of material time t; 0 otherwise
S
retrieved from a stockpile (Bley et al., 2012), but Moreno, Reza- Ynbdt fraction of bin n in block b extracted by the end of time
khah, Newman, and Ferreira (2016) show that there are both more t and sent to destination d
S
accurate and more tractable methods for modeling inventory; they Inbt fraction of bin n from block b in the stockpile at the
also show that the warehouse model indeed provides a reason- end of time t
S−
able approximation of reality, within a few percentage points of Inbt fraction of bin n from block b sent to the mill from
the “real” net present value for representative data sets (not unlike stockpile at the beginning of time t
our own). For ease of exposition and to enable us to use a spe-
300 B. King et al. / European Journal of Operational Research 257 (2017) 297–309


(S ) max δt cnb
S+
((Ynb1
S
t − Ynb1,t−1 ) + Inbt )
S S− must be completed before extraction on the given level can be-
b∈B n∈Nb t∈T gin. The method advances such that mining proceeds away from
   an initial stope determined a priori. The ore contained within a sill
− S−
cnb (Ynbdt
S S
− Ynbd,t−1 ) (1a) pillar (Fig. 4) is partially sterilized and can only be recovered, with
b∈B n∈Nb d∈D t∈T
significant dilution, at the end of the mine life. Sill pillar place-
  ment must balance the sterilization of ore with the increase in net
s.t. S
Ynbd,t−1 ≤ S
Ynbdt ∀b ∈ B, n ∈ Nb , t ∈ T (1b) present value gained by earlier access to stopes.
d∈D d∈D The precedence relationships for an underground mine that
uses this sequencing method can be categorized as follows: (i)
S
Xb.t−1 ≤ XbtS ∀b ∈ B, t ∈ T (1c) fixed predecessors include the development required to access the
stope, and stopes on the same level. These predecessors ensure
  that with respect to the mining direction, the adjacent stope on the
YnS−1,bdt = S
Ynbdt ∀b ∈ B, n ∈ Nb , t ∈ T (1d) same level be fully extracted and backfilled before extraction of the
d∈D d∈D next stope can commence. If no adjacent stope on the level exists,
then the stope has only development activities as predecessors. (ii)

XbtS ≤ S
Ynbdt ∀b ∈ B, n ∈ Nb , t ∈ T (1e) Conditional stope predecessors require that the stope directly be-
d∈D
low and the stopes on either side of the given stope on the level
below must be fully extracted and backfilled before the given stope
 can be extracted. If the stopes on the level below act as a sill pillar,
S
Ynbdt ≤ XˆS ∀b ∈ B, n ∈ Nb , bˆ ∈ Bˆb , t ∈ T (1f)
bt then the conditional predecessors are omitted (Fig. 5). In addition,
d∈D
every underground activity has a set of predecessor activities that
are dictated by the mine design.
S
Inb,t+1 S
= Inbt S−
− Inbt + (Ynb2
S
t − Ynb2,t−1 )
S
∀b ∈ B, n ∈ Nb , t ∈ T (1g) The underground mine scheduling model, (U), determines sill
pillar placement and a life-of-mine schedule consisting of devel-
 opment, stoping, and backfilling activities to maximize the under-
r rt ≤ qSrnb (Ynbdt
S S
− Ynbd,t−1 ) ≤ r̄rt ground mine’s NPV. This model precludes specific pairs of activities
b∈B n∈Nb d∈D from being completed in the same time period, and resource con-
∀r ∈ R  r = 1 , t ∈ T (1h) straints limit development, extraction, and backfilling. We assume
fixed activity rates.
 We maintain the same notation style as in the surface model,
r rt ≤ qSrnb ((Ynb1
S
t − Ynb1,t−1 ) + Inbt ) ≤ r̄rt
S S−
but use the superscript U to represent underground-specific pa-
b∈B n∈Nb
rameters and variables; we use checks and bars as accents.
∀r ∈ R  r = 2 , t ∈ T (1i)
Indices and sets:
 a∈A set of all activities
( XbtS − S
Xb,t−1 ) ≤ r̄rt ∀ p ∈ P, r ∈ R  r = 3, t ∈ T (1j)
s∈S⊂A set of stoping activities
b∈B p
ǎ ∈ Ǎa set of fixed predecessors for activity a
S S S− ā ∈ Āa set of fixed predecessor activities ā that must be
0 ≤ Ynbdt , Inbt , Inbt ≤ 1; XbtS binary
completed one time period
∀b ∈ B, n ∈ Nb , d ∈ D, t ∈ T (1k) in advance of activity a
The objective (1a) maximizes discounted revenue associated š ∈ Šs ⊂ Ǎs set of conditional predecessors š below stope s
with mill profits, and mining costs. Constraints (1b) and (1c) en- s̄ ∈ S̄s ⊂ Ās set of conditional predecessor stopes s̄ that must
sure that once a bin-block combination is completed, it remains be completed one or more
completed. Constraint (1d) precludes selective mining of any bin time periods in advance of stope s
in a block, i.e., the constraint forces all bins to be mined in equal l∈L levels in the mine
proportion. Constraint (1e) relates the fractional and binary ex- l ∈ Ľs level on which stope sexists (set has cardinality
traction variables. Constraint (1f) enforces precedence by prevent- of one)
ing the extraction of a block until its predecessors’ blocks have r∈R resources (4 = mine/mill capacity, 5 = backfill
been fully extracted. Constraint (1g) balances the inventory in the capacity, 6 = development capacity )
stockpile at the end of every time period. Constraint (1h) limits Data:
the capacity for extraction tonnes in each time period. Constraint caU monetary value associated with completing activity
(1i) bounds processing at the mill in each time period. Constraint a [$]
(1j) prevents mining too rapidly in one phase. Constraint (1k) en- qU
ra quantity of resource r associated with activity a [2, 4
forces non-negativity and integrality of the decision variables, as and 5 = tonnes, 6 = meters]
appropriate. r rt , r̄rt minimum, maximum level of resource r in time t
[4 − 5 = Million tonnes/year, 6 = months/year]
2.2. Underground model δt discount factor for time period t [fraction]
Decision variables:
Our underground formulation incorporates a mine design based
on pre-constructed stope shapes, organized into vertical levels.
XatU 1 if activity a is finished by the end of time t;
Drifts, i.e., tunnels that are only open at one end, are used to ac-
0 otherwise
cess the mine. A vertical decline is a drift that descends from the
WlU 1 if level l serves as a sill pillar; 0 otherwise
surface to the lowest underground level. On each level, horizon-

tal drifts are constructed from the decline to the stope locations. (U ) max δt caU (XatU − Xa,t−1
U
) (2a)
Our method sequences stopes from the bottom up such that ex- a∈A t∈T
traction and backfilling on the level underneath the given level
 B. King et al. / European Journal of Operational Research 257 (2017) 297–309 301

Fig. 4. (Regular grid of stopes) Sill pillar levels are black, and each block constitutes a designed stope. Any horizontal level of stopes shown in the figure could act as a sill
pillar.

completed. Constraint (2c) enforces fixed precedence, and con-

straint (2d) enforces conditional precedence based on sill pillar
placement. Constraint (2e) ensures that at least one time period
elapses between the completion of the specified pair of activities.
Constraint (2f) forces at least one time period to elapse between
the completion of two given stoping activities unless a sill pillar
exists on a level in between them; these stoping activities need
not be on consecutive levels because the precedence might actu-
ally allow stopes on consecutive levels to be mined in the same
time period. Constraint (2g) prevents mining the stopes on level
l, if level l acts as a sill pillar. Constraint (2h) bounds extraction,
backfill, and development resource use. (Mill processing capacity
is essentially unconstrained underground because of the low pro-
Fig. 5. (Sequencing method) The initial stope location is labeled, with a black ar- duction rate.) All variables are required to be binary by (2i).
row, and stope X’s predecessors are denoted by gray arcs. On the left is an example
of standard precedence, i.e., both fixed and conditional predecessors. On the right
Delay constraints (2e) and (2f) capture sub-annual detail in our
is a stope that has only fixed predecessors, because the sill pillar eliminated all of model with annual time fidelity, and exclude two specified activi-
the conditional predecessors. ties (a , a) from being completed in the same time period, where
a is a predecessor of a, and the minimum time required to elapse
between the start of a and the completion of a is greater than the
U
s.t. Xa,t−1 ≤ XatU ∀ a ∈ A, t ∈ T (2b) time fidelity of the model. It is important to construct a minimal
number of delay constraints so as not to unnecessarily increase the
number of precedence constraints.
XatU ≤ XǎUt ∀ a ∈ A, ǎ ∈ Ǎa , t ∈ T (2c)

2.3. Basic transition model

XstU ≤ XšUt + WlU ∀ s ∈ S, š ∈ Šs , l ∈ Ľš , t ∈ T (2d)
The basic transition model, (Tb ), may be formulated by com-
bining of the surface, (S), and underground, (U), models. The ob-
XatU ≤ XāU,t−1 ∀ a ∈ A, ā ∈ Āa , t ∈ T (2e) jective function maximizes the NPV of the combined open pit and
underground operations, i.e., the sum of (1a) and (2a). A vast ma-
 jority of the constraints, (1b)–(1h), (1j), (1k), and (2b)–(2i), remain
XstU ≤ Xs̄U,t−1 + W jU ∀ s ∈ S, s̄ ∈ S̄s , i ∈ Ľs̄ , k ∈ Ľs , t ∈ T (2f) the same as in their respective models. The precedence constraints
i≤ j≤k in both the open pit and underground mine do not need to be
changed based on the crown pillar location, because of the direc-
XstU + WlU ≤ 1 ∀ s ∈ S, l ∈ Ľs , t ∈ T (2g) tion of extraction for each method. Any non-zero lower bounds
on the underground mine’s resource constraints are removed to
allow for a delayed start of the underground mine. The resource

∀ r ∈ R  r ≥ 4, t ∈ T constraints must be altered to accurately reflect that the open
ra (Xat − Xa,t−1 ) ≤ r̄rt
qU U U
r rt ≤ (2h)
a∈A
pit, stockpile, and/or underground mine may be sending material
to the mill in the same time period. Constraints associated with
the additional variables, i.e., that represent the crown pillar loca-
XatU , WlU binary ∀a ∈ A , l ∈ L, t ∈ T (2i)
tion, preclude any open pit or underground extraction of mate-
The objective function (2a) maximizes net present value. Con- rial located in the crown pillar. Additional notation is shown be-
straint (2b) ensures that once an activity is completed, it remains low with the superscript Tb representing transition model-specific
302 B. King et al. / European Journal of Operational Research 257 (2017) 297–309

variables.

Indices and sets:

v∈V set of crown pillar elevations
b˜ ∈ B˜v set of blocks that exist below the crown pillar if the
crown pillar is located at elevation v Fig. 6. (Special knapsack reformulation) An additional precedence arc, dashed, is
a˜ ∈ A˜v set of activities that exist above the crown pillar if the added to prevent block p˜ and the successor block b from being completed in the
crown pillar is located at elevation v same time period, because their precedence separates them by more blocks than
can be completed in a time period. Immediate precedence is shown with solid arcs.
Decision variables:
WvT
b
1 if the crown pillar is located a elevation v; 0 otherwise
constraints, orders of magnitude faster than simplex-based meth-
ods. With reformulation and an ad-hoc branch-and-bound strategy,

(Tb ) max δt cnb
S+
((Ynb1
S
t − Ynb1,t−1 ) + Inbt )
S S− we are able to identify open pit-to-underground transition options
b∈B n∈Nb t∈T with near-optimal NPVs.
   We reformulate the basic transition model (Tb ) by transforming
− δt cnb
S−
(Ynbdt
S S
− Ynbd,t−1 ) some of the side constraints into precedence constraints, specifi-
b∈B n∈Nb d∈D t∈T
 cally, a special knapsack, (1j), and the “warehouse-style” inventory,
+ δt caU (XatU − Xa,t−1
U
) (3a) (1g), constraints in the surface model, (S). Mathematical proofs
a∈A t∈T showing that these reformulations are no weaker than the origi-
nal formulations can be found in Appendix A.
s.t. Xb˜St ≤ 1 − WvT ∀b˜ ∈ B˜v , v ∈ V, t ∈ T
b
(3b)
3.1. Special knapsack reformulations

Xa˜Ut ≤ 1 − WvT ∀a˜ ∈ A˜v , v ∈ V, t ∈ T

b
(3c) We show how to transform sinking rate constraints, (1j), into
precedence constraints by exploiting the facts that: (i) the blocks
within a phase are required to be completed in a fixed order,

WvT = 1
b
(3d) i.e., blocks must be extracted in sequential order from the surface
v∈V downwards, and (ii) the left-hand side is 0 for all time periods
in all scenarios. Constraint (1j) from the initial surface model (S)
 appears as follows:
r rt ≤ qSrnb ((Ynb1
S
t − Ynb1,t−1 ) + Inbt )
S S−

b∈B n∈Nb (XbtS − Xb,t−1
S
) ≤ r̄rt ∀ p ∈ P, r ∈ R  r = 3, t ∈ T (4)

+ ra (Xat − Xa,t−1 ) ≤ r̄rt
qU U U
∀r ∈ R  r = 2 , t ∈ T (3e) b∈B p

a∈A The reformulation of constraint (4) prevents a block b that is r̄3t

successor blocks away from the selected block p˜ in the phase from
WvT binary ∀v ∈ V
b
(3f) being completed in the same time period as block b. This requires
Retained constraints from (S ): (1b )–(1h ), (1j ), (1k ) the following set definition:

Retained constraints from(U ): (2b)–(2i ). p˜ ∈ P˜ b predecessors for block b that must be completed at
The objective function (3a) maximizes net present value of the least one time period prior to block b
entire deposit and replaces (1a) and (2a). Constraint (3b) allows for The reformulation is shown in (5):
open pit mining to only occur above the crown pillar. Constraint
(3c) restricts underground mining to only occur below the crown XbtS ≤ X pS˜,t−1 ∀b ∈ B, p˜ ∈ P˜ b , t ∈ T (5)
pillar. Constraint (3d) forces the placement of a crown pillar. Con-
This constraint set, (5), has far greater cardinality than (4), but
straint (3e) replaces constraint (1i) with respect to mill capacity.
possesses precedence structure. Fig. 6 shows an example of the
constraint construction.
3. Reformulations
3.2. Inventory balance reformulations
The basic transition model, (Tb ), is theoretically NP-hard, and,
in practice, real-world size problems are intractable with current
In this section, we present a reformulation of the inventory bal-
computer hardware and software. Our scenarios contain nearly
ance constraint (1g):
50,0 0 0 variables and more than 1.5 million constraints, even after
efficient variable elimination techniques are used (Lambert et al., S
Inb,t+1 S
= Inbt S−
− Inbt + (Ynb2
S
t − Ynb2,t−1 )
S
∀b ∈ B, n ∈ Nb , t ∈ T (6)
2014; O’Sullivan, 2013).
Bienstock and Zuckerberg (2010) provide an algorithm, the “BZ Our reformulation implies that material must be placed in in-
algorithm,” for efficiently solving the LP relaxation of problems ventory before it is processed in the same or in a later time period,
with the math structure seen in PSPCP, i.e., a model in which a and is mathematically equivalent to the original under the assump-
majority of the constraints are precedence, rather than “side,” e.g., tion that there is no value lost for placing material in inventory,
knapsack, constraints. In practice, the BZ algorithm’s solution time i.e., there is no mixing, degradation, or rehandling cost associated
is more sensitive to the latter type of constraints than to the num- with placing or retrieving material. We require the following vari-
ber of precedence constraints. Muñoz et al. (2015) provide an im- able definitions:
plementation framework for solving the LP relaxation of open pit YˆbtS fraction of block b extracted and able to be processed
mining problems using the BZ algorithm, and show that it is pos- by the end of time t
sible to obtain LP relaxation solutions to PCPSPs with millions of S
Znbt fraction of bin n in block b sent to the mill by the end
variables and precedence constraints, but fewer than 200 “side” of time t
 B. King et al. / European Journal of Operational Research 257 (2017) 297–309 303

S
We replace all instances of variables Ynbdt with the appropri- XbtS ≤ YˆbtS ∀b ∈ B, t ∈ T (10f)
ˆ S S
ate Ybt or Znbt variables, where the former newly introduced vari-
able represents the fraction of a block extracted in time t, without
recognizing the destination. The latter newly introduced variable YˆbtS ≤ XˆS
bt
∀b ∈ B, bˆ ∈ Bˆb , t ∈ T (10g)
S tracks both the processing time period and destination of each
Znbt
bin-block combination. If both variables for a given bin-block com-
bination assume a value of 1 in the same time period, that bin-
S
Znbt ≤ YˆbtS ∀b ∈ B, n ∈ Nb , t ∈ T (10h)
block combination is immediately sent to the mill for processing
after extraction. For all periods in which a specific bin-block com-
bination is in the stockpile, the variable representing that block’s
XatU ≤ X p̄Ut ∀a ∈ A, p̄ ∈ P̄a , t ∈ T (10i)
extraction, YˆbtS , assumes a value of 1 and the corresponding vari-
able representing processing, ZnbtS , assumes a value of 0. Any bin-
YˆbtS ≤ Yˆp˜S,t−1 ∀b ∈ B, p˜ ∈ P˜b , t ∈ T (10j)
block combination that is extracted and not processed is sent to
the waste dump, resulting in all corresponding ZnbtS variables pos-
sessing a value of 0. The reformulation of (6) is shown in (7):
XatU ≤ X p̄U,t−1 ∀a ∈ A, p̄ ∈ P̄a , t ∈ T (10k)
S
Znbt ≤ YˆbtS ∀b ∈ B, n ∈ Nb , t ∈ T (7)
Constraint (7) allows only material that has been extracted to 
r rt ≤ qSrnb (YˆbtS − Yˆb,t−1
S
) ≤ r̄rt ∀r ∈ R  r = 1 , t ∈ T (10l)
be sent to the mill. This reformulation also requires substituting
S S− b∈B n∈Nb
the variables Ynbdt and Inbt in the objective function and mill ca-
pacity constraints with Yˆ S and Z S :  
bt nbt r rt ≤ qSrnb (Znbt
S S
− Znb,t−1 )+ ra (Xat − Xa,t−1 ) ≤ r̄rt
qU U U

max δ S+
t cnb ( S
Znbt − S
Znb,t−1 ) b∈B n∈Nb a∈A

b∈B n∈Nb t∈T ∀r ∈ R  r = 2 , t ∈ T (10m)


− δ S−
t cnb (
YˆbtS S
− Yˆb,t−1 ) (8)

b∈B n∈Nb t∈T
r rt ≤ ra (Xat − Xa,t−1 ) ≤ r̄rt
qU U U
∀r ∈ R  r ≥ 4 , t ∈ T (10n)
and a∈A

r rt ≤ qSrnb (Znbt
S S
− Znb,t−1 ) ≤ r̄rt ∀r ∈ R  r = 2 , t ∈ T (9)
b∈B n∈Nb
XatU , XbtS binary ∀a ∈ A, b ∈ B, t ∈ T (10o)

For all other constraints, the variable YˆbtS replaces Ynbdt

S ; and, the
S and I S− are eliminated.
variables Inbt nbt
0 ≤ YˆbtS , Znbt
S
≤1 ∀b ∈ B, n ∈ Nb , t ∈ T (10p)
The objective function (10a) maximizes net present value, and
3.3. Enhanced transition model
replaces the objective function (3a). Constraints (10b)–(10e) ensure
that each completed activity or block remains completed, and are
The enhanced transition model is the combination of the refor-
a substitute for (1b), (1c), and (2b). Constraints (10f) and (10g) en-
mulated surface (S) and underground (U) mine scheduling models,
force the precedence structure for the open pit mine, and replace
in which the crown and sill pillar placements are fixed a priori. All
(1e) and (1f). Note that the replacement constraints do not sum
other constraints are similar to those in the basic transition model,
on the destination index because all material is sent to a stock-
(Tb ). The hat and tilde accents are reserved for open pit sets, and
pile, even if just instantaneously. Constraint (10h) ensures that the
bar accents for underground sets. Additional notation and the en-
fraction of a bin that is sent to the mill is no greater than the frac-
hanced transition model (Te ) follow:
tion extracted from the corresponding block, and is a reformulation
Indices and sets: of (1g). Constraint (10i) enforces underground mine precedence,
p̄ ∈ P̄a predecessors for activity a that must be completed and is used instead of constraints (2c), (2d), and (2g). Constraints
at least one time period in advance of activity a (10j) and (10k) ensure that one time period elapses between the
completion of two specific activities or blocks, and are a replace-

(Te ) max δt cnb
S+
(Znbt
S S
− Znb,t−1 ) ment for constraints (1j), (2e), and (2f). Constraint (10l) bounds
b∈B n∈Nb t∈T
open pit-specific resource use, and is a substitute for constraint
 (1h). Constraint (10m) bounds the mill capacity, and is a substitute
− δt cbS− (YˆbtS − Yˆb,t−1
S
) (10a) for constraint (3e). Constraint (10n) bounds underground-specific
b∈B t∈T resource consumption, and is equivalent to constraint (2h). Con-
 straints (10o) and (10p) enforce binary and variable bounds, where
+ δt caU (XatU − Xa,t−1
U
) appropriate.
a∈A t∈T
S
s.t. Xb,t−1 ≤ XbtS ∀b ∈ B, t ∈ T (10b) 4. Solution strategy

We obtain near-optimal solutions for the enhanced transition

S
Yˆb,t−1 ≤ YˆbtS ∀b ∈ B, t ∈ T (10c) model, (Te ), presented in Section 3.3, by: (i) exhaustively search-
ing possible crown and sill pillar placement options using an ad-
hoc branch-and-bound strategy and solving the resulting LP re-
S
Znb,t−1 S
≤ Znbt ∀b ∈ B, n ∈ Nb , t ∈ T (10d) laxations, (ii) using a rounding heuristic to convert the LP relax-
ation solutions with favorable objective function values into inte-
ger solutions, and (iii) using integer solutions to eliminate a num-
U
Xa,t−1 ≤ XatU ∀a ∈ A , t ∈ T (10e)
ber of possible crown and sill pillar placement options to reduce
the amount of computation required in (ii).
304 B. King et al. / European Journal of Operational Research 257 (2017) 297–309

The reformulation in Sections 3.1 and 3.2 reduces the num- on processing properties, and a waste bin. This results in a total of
ber of side constraints in the basic transition model, (Tb ), but the 1312 bin-block combinations ranging from 250 to 1,10 0,0 0 0 tonnes.
model is still not in the desired form for obtaining an efficient LP Extraction activities may possess as many as three immediate pre-
relaxation solution using the BZ algorithm. By fixing, i.e., branch- decessor activities. The cost of extraction increases as the depth of
ing on, all of the variables associated with the placement of the the open pit increases. Section 2.1 provides a detailed description
b
crown and sill pillars, WlU and WvT , respectively, we convert all of the open pit precedence and physical representation of the data.
of the conditional precedence constraints, (2d) and (2f), in the un- Our basic underground model dataset consists of 1123 devel-
derground model, (U), to standard precedence constraints, and the opment activities, 351 stoping activities and an equal number of
basic transition model at each node to a model with a PCPSP math- backfilling activities. Stopes range from approximately 50 0 0 to
ematical structure. We branch as follows: for each possible crown 40,0 0 0 tonnes, resulting in 17 levels in the underground mine that
pillar placement (which, for our data set, is twelve), we consider are a maximum of 40 meters in height. The required development
all sill pillar placements consisting of between zero and three such and backfilling is estimated based on the stope properties. Each
pillars, where three would be a maximum operationally feasible activity has up to 12 immediate predecessors and up to 100 de-
number. The total number of viable crown and sill pillar placement lay constraints. Section 2.2 provides a description of the under-
options numbers in the thousands. ground mine’s precedence structure. For our analysis, we construct
We first solve the LP relaxation of the transition model for each ten distinct scenarios, each for a 24-year time horizon with de-
set of reasonable crown and sill pillar placements, i.e., for each cisions made at yearly fidelity, and each defined as a set of up-
branch, using the BZ algorithm. We then sort these LP relaxation per bounds on the resource constraints and a given discount rate
solutions, decreasing by objective function value. Solutions with (Table 1). We set the underground backfilling capacity equal to the
the best LP relaxation objective function values are transformed underground extraction capacity.
into IP solutions using TopoSort (Chicoisne et al., 2012), which has
been shown to provide near-optimal solutions quickly for open pit 5.2. Numerical results
mine scheduling problems that only have non-zero upper bounds
on resource constraints, and which is based on the premise that We compare the performance of the OMP Solver (Muñoz et al.,
the earlier the expected completion time of a block or activity in 2015) to that of AMPL/CPLEX, (AMPL Optimization LLC, 2014; IBM
the LP solution, the earlier the block or activity is scheduled in CPLEX Optimizer, 2013), using a Dell PowerEdge R410 with 16 pro-
the IP solution. The algorithm maintains precedence constraints cessors (2.72 gigahertz each) and 28 gigabytes of RAM. OMP, Ver-
by the fact that the expected completion time of a block or ac- sion 1509 is an academic, customized solver that uses standard
tivity in the LP relaxation is always greater than or equal to that preprocessing and exploits the mathematical structure of PCPSP
of its predecessors. Also employed in our variant of TopoSort is to solve the LP relaxation quickly using the BZ algorithm; then,
an “alpha points” procedure in which activities are ordered not by we execute the TopoSort heuristic eleven times, each with a dif-
their expected completion time, but by the time period in which ferent alpha point value between 0 and 1, inclusive, incremented
a specified fraction, i.e., alpha point, of the activity has been com- by 0.1; this procedure transforms the LP relaxation to an integer-
pleted. Therefore, an alpha point of 0.7 would set the order based feasible solution, of which we choose the best one. All other pa-
on the first time period in which the “by” variable obtains a value rameter settings are default. CPLEX 12.6.0.0 uses default parameter
larger than 0.7. The TopoSort heuristic allows for us to match the settings other than memory emphasis, and a 40,0 0 0-second time
(Te ) formulation exactly, i.e., create a mixed integer solution to the limit. Variable elimination techniques are employed before passing
open pit portion, and a fully integral solution to the underground the model to CPLEX (Lambert et al., 2014; O’Sullivan, 2013). Both
portion. Once an IP solution is obtained from the LP relaxations solvers provide solutions to the enhanced transition model, (Te ),
with the largest objective function values, we use bound domi- with the same crown and sill pillar placement options available
nance to eliminate a significant number of the crown and sill pillar in each scenario. Depending on the crown and sill pillar place-
placement options. Specifically, every LP relaxation whose objec- ment, for our dataset, the enhanced transition model, (Te ), av-
tive function value is less than that of an existing feasible IP so- erages 50,0 0 0 variables and 1.5 million constraints. (The numer-
lution’s cannot correspond to an optimal integer placement of the ical performance of (Tb ) is dominated by that of (Te ) using our
crown and sill pillar. methodology; see Appendix A.)
We first compare the performance using CPLEX to solve the en-
5. Data and numerical results hanced transition model (Te ) for a fixed crown and sill pillar loca-
tion (giving CPLEX the benefit of the faster LP solver) against that
We introduce the data required for the enhanced transition of the OMP Solver. For a representative crown and sill pillar place-
model. Computational results highlight the speed, effectiveness, ment given as the ordered pair [(820), (460)], where these eleva-
and robustness of the methodology which yields consistent near- tions are relative to sea level, the enhanced transition model, (Te ),
optimal solutions to our multiple scenarios of the transition model. contains approximately 60,0 0 0 variables and 1,20 0,0 0 0 constraints,
of which 120 are “side” constraints. CPLEX averages 163.77 seconds
5.1. Data with the faster LP solver for each LP relaxation over the ten scenar-
ios, and produces slightly better integer solutions in only two of
Our industry partner provided all of the data required for the the ten scenarios (while CPLEX is unable to find an integer-feasible
transition model from an active mine in Africa; grade and cost data solution in the other scenarios due either to memory or time lim-
are confidential. The deposit is known to extend over a large ver- itations). By contrast, OMP is able to solve the LP relaxations in
tical expanse, and the overlap between the upper-most designed fewer than ten seconds, regardless of the scenario, and produces
stope and lowest-planned open pit extraction elevation is over 400 an integer solution within 6 percent of optimality or better in just
meters; nearly 80 percent of the remaining recoverable material is a few additional seconds (Table 2).
located in this overlap, or transition zone. Fig. 7 depicts the LP relaxation objective function value and
The open pit dataset consists of a four-phase design for a par- the best known IP objective function value for each reasonable
tially extracted open pit mine with a total of 336 blocks ranging set of crown and sill pillar placements for Scenario 1 using the
in weight from 20,0 0 0 to 5,50 0,0 0 0 tonnes. Blocks may contain a OMP Solver. Both the LP relaxation objective function value and
high-, medium-, and low-grade bin for two material types based the best-known IP objective function value follow the same trend
 B. King et al. / European Journal of Operational Research 257 (2017) 297–309 305

Table 1
(Scenario summary) Capacities and discount rates used in each scenario, with equation numbers from the enhanced transition
model also given in the column headers.

Scenario Annual capacities Annual

Extraction Development discount

Open pit (tonnes) (10l) Underground (tonnes) (10n) (meters) (10n) Mill (tonnes) (10m) rate (percent)

1 50,0 0 0,0 0 0 2,0 0 0,0 0 0 50 0 0 8,0 0 0,0 0 0 9

2 50,0 0 0,0 0 0 2,0 0 0,0 0 0 50 0 0 7,0 0 0,0 0 0 9
3 50,0 0 0,0 0 0 2,0 0 0,0 0 0 2500 8,0 0 0,0 0 0 9
4 50,0 0 0,0 0 0 1,50 0,0 0 0 50 0 0 8,0 0 0,0 0 0 9
5 40,0 0 0,0 0 0 2,0 0 0,0 0 0 50 0 0 8,0 0 0,0 0 0 9
6 40,0 0 0,0 0 0 2,0 0 0,0 0 0 50 0 0 7,0 0 0,0 0 0 9
7 50,0 0 0,0 0 0 2,0 0 0,0 0 0 50 0 0 6,0 0 0,0 0 0 9
8 50,0 0 0,0 0 0 1,50 0,0 0 0 2500 8,0 0 0,0 0 0 9
9 50,0 0 0,0 0 0 2,0 0 0,0 0 0 50 0 0 8,0 0 0,0 0 0 1
10 50,0 0 0,0 0 0 2,0 0 0,0 0 0 50 0 0 8,0 0 0,0 0 0 15

Table 2
(OMP and CPLEX comparison) Comparison of solution times and optimality gaps between the OMP Solver and CPLEX for (Te ). All sce-
narios are run with a crown pillar located at elevation 820 and a sill pillar located at level 460.

Scenario CPLEX OMP solver

LP solution time (seconds) IP solution Optimality LP solution TopoSort solution Optimality

Barrier Simplexa time (seconds) gap (percent) time (seconds) time (seconds) gap (percent)

b
1 163.75 488.41 – 9.72 3.59 2.78
b
2 163.68 490.24 – 5.23 3.13 3.56
b
3 177.85 1631.15 – 8.43 3.29 4.23
c
4 183.08 577.07 – 6.42 3.67 3.90
5 146.24 613.93 26,100 2.66 6.28 3.62 4.49
6 172.88 783.70 34,728 2.67 6.60 3.44 5.86
c
7 152.41 416.32 – 5.44 3.12 4.76
c
8 152.00 1300.20 – 8.57 2.14 4.71
b
9 163.26 707.38 – 4.80 2.41 0.64
b
10 162.59 653.63 – 5.10 2.87 4.24
  IP Obj. Func. Value 
Note: Optimality gaps are calculated as 100% · 1 − .
LP Relaxation Value
a
CPLEX is allowed to choose the variant of simplex to use, which results in employing dual simplex on the dual problem.
b
CPLEX was unable to produce an integer solution within a 5 percent gap before running out of memory.
c
CPLEX was unable to produce an integer solution within a 5 percent gap before the 40,0 0 0 second limit.

Fig. 7. (LP and IP comparison) Left: LP relaxation values for all feasible crown and sill pillar placement options. Right: best-known IP objective function value for the
corresponding crown and sill pillar placement options. The vertical band of crosses at each crown pillar elevation is associated with the scaled NPV corresponding to all
viable sill pillar location combinations. A horizontal line indicates overall best-known IP objective function value for Scenario 1.

as we exhaustively enumerate all of the 3500 crown and sill pillar largest objective function value produces the largest IP objective
placement options. The average gap between the LP relaxation function value, suggesting empirically that our solution methodol-
objective function value and the best-known IP objective function ogy provides consistently high-quality IP solutions relative to the
value is, on average, 3.91 percent, and the time to obtain the LP solutions for the enhanced transition model, (Te ).
integer solution is, on average, 9.78 seconds. This gap also appears Our ad-hoc branch-and-bound strategy supplies a wealth of in-
to be relatively consistent across all of the crown and sill pillar formation for the mine operator: crown pillar placement affects
placement options for this scenario. The LP relaxation with the the NPV significantly more than sill pillar placement. Additional in-
306 B. King et al. / European Journal of Operational Research 257 (2017) 297–309

Fig. 8. (Zoom comparison) Left: LP relaxation values associated with crown and sill pillar placements that produce a high objective function relaxation value, and a horizontal
line representing the best-known IP objective function value. Right: corresponding IP objective function values for models whose LP relaxation objective function value is
greater than the best known IP objective function value. Note: circled is the best LP relaxation objective function value and its corresponding IP objective function value.

sights might involve geology: if it is undesirable to have a crown Table 3

(Scenario summary) Optimal solution, integrality gap, and total solution
pillar located at elevation 820, moving the crown pillar to elevation
time for each scenario if enumerated crown and sill pillar placement
780 would have the least impact on the mine’s NPV (Fig. 7). solves are performed in serial.
It is possible to heavily prune our ad-hoc branch-and-bound
Scenario Optimal crown and Integrality Total solution
tree using bound dominance. For example, for Scenario 1, after we
number sill pillar placement gap (percent) time (seconds)
obtain an IP objective function value associated with the crown
and sill pillar placement option that has the highest LP relaxation 1 [(820), (500)] 2.51 29,652
objective function value, we can eliminate solving the integer pro- 2 [(820), (420)] 3.28 28,220
3 [(820), (500)] 3.25 30,642
gram corresponding to all crown and sill pillar placements whose
4 [(820), (660)] 3.42 31,426
LP relaxation objective function value is lower. Only 40 of the over 5 [(820), (500)] 4.40 25,993
3500 crown and sill pillar placement options have an LP relax- 6 [(820), (420)] 4.97 27,404
ation objective function value greater than the best known IP ob- 7 [(820), (460)] 4.76 36,312
jective function value (Fig. 8). The mine planner interested in ro- 8 [(820), (500)] 3.45 41,652
9 [(700),(420)] 0.70 26,016
bust solutions might note that of those options, only one of them
10 [(820), (500,660)] 3.65 23,825
is not associated with a crown pillar located at elevation 820, and
the corresponding LP relaxation’s objective function value is only Notes: Crown and sill pillar placement option format is [(Crown pillar
elevation), (sill pillar elevation(s))]. Total solution time is the time re-
0.13 percent greater than the best-known IP objective function quired for the LPs associated with all possible crown and sill pillar lo-
value. cations, and the additional time to obtain an IP solution for the non-
Table 3 summarizes the near-optimal crown and sill pillar dominated LP relaxations.
placement options associated with each scenario. The average gap
between the LP and IP objective function values is 5.55 percent. For
any scenario with a discount rate of 9 percent, the crown pillar as- narios would require fewer than ten seconds to solve with the ap-
sociated with the highest IP objective function value is located at propriate hardware; our methodology efficiently provides a way to
the same elevation, 820. (Changing the discount rate affects the identify near-optimal crown and sill pillar placements where no
best-known crown and sill pillar locations.) Techniques may be such methodology had existed.
employed to reduce the gaps, but, for the scenarios we tested, it Qualitatively, we can establish some generalities about the
is unlikely that such refinements would lead to solutions with a schedules for each scenario. The crown pillar location that is cho-
change to the crown pillar placement, because of the 40 solutions sen in a majority of the scenarios, 820, contains the fifth-highest
not eliminated by bound dominance, only one had a placement at amount of metal, and, as such, does not correspond to an intuitive
a level other than 820 (see Fig. 8). We report solution time as the solution of minimizing lost metal. Additionally, for the best sched-
CPU time required to solve the LP relaxations associated with all ule we report for each scenario, underground construction and
reasonable crown and sill pillar placement options, plus that re- production begins as soon as possible owing to that fact that all
quired to solve the necessary integer programs, i.e., those not ex- underground mine production is sufficiently high grade that it dis-
cluded by bound dominance. However, our procedure is massively places material from the open pit at the mill. We do observe some
parallelizable in that all LPs can be solved simultaneously, as can fluctuations in both the open pit and underground production,
all relevant IPs. Hence, on average, even the longest-running sce- which is undesirable from an operational standpoint, and would
 B. King et al. / European Journal of Operational Research 257 (2017) 297–309 307

require smoothing to create an operationally feasible schedule. is weak relative to that of the enhanced transition model (Te ), and
However, these fluctuations are not uncommon in a strategic plan. we compare LP solution times across standard algorithms. We also
show that CPLEX is unable to solve the basic transition model (Tb )
6. Conclusions and future work for any scenario, i.e., that the basic transition model is intractable
when solved with CPLEX. Furthermore, two proofs – one each
The methodology developed in this paper provides a robust for the special knapsack and inventory balance constraints (see
framework for solving a linear-integer program representing an Section 3) – show that our reformulations are no weaker than the
open-pit-to-underground transition model involving scenarios that original ones; moreover, computational results show that the re-
contain 50,0 0 0 variables and over 1.5 million constraints. An ad- formulations are strictly stronger.
hoc branch-and-bound scheme fixes the variables that destroy the
PCPSP structure without compromising optimality. This methodol- A1. Solutions to the basic transition model (Tb )
ogy permits us to test a wide variety of scenarios quickly and pro-
vides a better understanding of how crown and sill pillar place- Table 4 provides specific information regarding algorithmic per-
ment affects NPV. formance and solution quality for our ten test scenarios. We first
With our specialized technique, we are able to solve a relevant compare solution times for standard LP algorithms when used to
and economically significant problem for the mining industry. As solve (Tb ). All computations were run on the same server as the
current open pit mines are required to extract an increasing num- enhanced transition model (Te ). The barrier outperforms the best
ber of tons of waste material for every ton of ore, it becomes cru- version of simplex in all cases; the academic nature of the OMP
cial to identify the proper transition location. Although, for confi- solver precludes it from handling the complexity in the basic tran-
dentiality reasons, the exact NPVs are not given, our results show sition model, (Tb ), in particular, the decisions regarding crown and
that the NPV can change by hundreds of millions of dollars de- sill pillar placement (hence, our ad-hoc branch-and-bound strat-
pending on the crown pillar placement, and by tens of millions egy). Despite the ability of CPLEX (version 12.6.0.0 using default
based on the sill pillar placement. Our model provides not only a parameter settings other than turning on memory emphasis) to
near-optimal solution, but identifies the economic outcome of all solve the LP relaxation of (Tb ), the complexity of the problem
possible crown and sill pillar placements. proves intractable when seeking a good, integer-feasible solution.
Many mine operators defer underground mining until the open For all our scenarios, CPLEX exhausts a 40,0 0 0-second time limit
pit has finished production, resulting in insufficient cash flow to or runs out of memory before a solution within 5 percent of op-
justify an underground mine and an unmined portion of the de- timality is found. We attribute this poor performance, in part, to
posit that could have been extracted economically. By develop- the weak LP bound (as seen with a comparison of the scaled net
ing an efficient and tractable solution methodology, we can pro- present values in the penultimate and last columns of Table 4). By
vide mine operators with a tool to better understand the benefits contrast, the results we obtain from using our reformulations and
of each transition elevation, and the ability to confidently make a solution procedure on (Te ) result not only in tighter bounds, but
timely decision. also in near-optimal integer solutions, the latter stemming in large
Additional work could address: (i) accuracy, (ii) applicability, part from the mathematical structure of (Te ) our TopoSort heuris-
and (iii) optimality gap. The accuracy of the model would be im- tic is able to exploit. The values for the tight LP relaxation solution
proved with a better representation of the stockpiles. Since stock- we obtain from (Te ) in Table 4 result from the fixed crown and sill
piles contribute significantly to the NPV, it would be beneficial to pillar combination that gives the best objective function value for
include mixing of ore in the stockpile and the degradation of ore that scenario.
grade over time. The applicability of the model could be improved
by adding blending requirements at the mill and non-zero lower A2. On the strength of the special knapsack reformulation
bounds on the knapsack constraints, which can be vital to maintain
proper mill feed, but that would destroy the mathematical struc- In this section, we prove that the reformulation of the special
ture that the TopoSort heuristic relies on. Finally, we wish to in- knapsack constraint (1j), presented in Section 3.1, is not weaker
corporate a branch-and-bound algorithm within the OMP solver to than the original formulation.
reduce the optimality gap for a fixed crown and sill pillar place- For simplicity, but without loss of generality, we consider the
ment option. case of a single phase. For notation purposes, assume that the
blocks in the phase are numbered 1, . . . , |B | and that there is a
Acknowledgments limit of k blocks that can be extracted in a given time period. (In
Section 3.1, this term is represented as r̄.) We assume variables xbt
Alexandra Newman and Barry King received funding from the are defined as before:
Center for Innovation in Earth Science and Engineering at the Col- 
1 if block b is extracted by time t
orado School of Mines, and from Alford Mining Systems. Marcos xbt =
0 otherwise.
Goycoolea received funding from FONDECYT grant #1151098 and
CONICYT PIA Anillo grant 1407. The authors acknowledge Mon- We can model the bound of k extractable blocks in a time pe-
ica Dodd, Chris Alford, Xiaolin Wu, Hongliang Wang, Ralf Kintzel, riod either by adding the knapsack constraint:
Conor Meagher, and many others for their continued support of
|B |

and advocacy for this project. This research benefited significantly
from suggestions made by Daniel Espinoza (Universidad de Chile),
(xbt − xb,t−1 ) ≤ k, ∀t ∈ T  t > 1 (11)
b=1
Eduardo Moreno (Universidad Adolfo Ibañez), Orlando Rivera (Uni-
versidad Adolfo Ibañez), and Andrea Brickey (South Dakota School or by adding precedence constraints:
of Mines).
xbt ≤ xb−k,t−1 ∀b ∈ B  b ≥ k + 1 , t ∈ T  t > 1 . (12)

Appendix A We now show that (12) is at least as strong as (11).

Lemma. Let X be the set of xbt variables such that:

We demonstrate here computationally that, for the scenarios
we examine, the LP relaxation of the basic transition model (Tb ) xbt ≤ xb−1,t ∀t, b ∈ B  b ≥ 1 (13a)
308 B. King et al. / European Journal of Operational Research 257 (2017) 297–309

Table 4
(Basic and enhanced comparison) Comparison of solution times and solutions gaps between the
basic transition model (Tb ) and the enhanced transition model (Te ) using CPLEX.

Scenario LP solution time for (Tb ) (seconds) IP solution (Tb ) LP Largest (Te )

Barrier Simplexa time (seconds) value LP value

b
1 373.98 10343.09 3.20 2.72
b
2 359.53 8222.94 3.11 2.66
c
3 992.07 13821.83 2.71 2.55
c
4 333.09 13163.29 3.11 2.70
b
5 300.37 7451.82 3.20 2.70
b
6 338.96 8018.63 3.09 2.64
b
7 407.66 7792.96 3.01 2.58
c
8 756.81 25526.81 2.71 2.55
c
9 383.78 8240.55 6.24 4.86
b
10 369.52 7738.03 2.19 2.06
a
CPLEX is allowed to choose the variant of simplex to use, which results in employing dual
simplex on the dual problem.
b
CPLEX was unable to produce an integer solution within a 5 percent gap before running out
of memory.
c
CPLEX was unable to produce an integer solution within a 5 percent gap before the 40,0 0 0
second limit.

A3. Inventory balance reformulation as variable substitution

xbt ≤ xb,t+1 ∀b ∈ B, t ≤ |T | − 1 (13b)

In this section, we show that the reformulation of the inven-
tory balance constraints (6) and the corresponding variable substi-
tutions into expressions (7)–(9) presented in Section 3.2, are not
0 ≤ xbt ≤ 1 ∀b ∈ B, t ∈ T . (13c)
weaker than the original formulation. To this end, it suffices to
For sets: show that for every integer-feasible solution of the reformulation,
|B | there exists a corresponding feasible solution to the original for-
 mulation having the same objective function value. Given a solu-
P1 = {x ∈ X : (xbt − xb,t−1 ) ≤ k ∀t ∈ T }
tion YˆbtS , Znbt
S of the reformulation, we construct a feasible solution
b=1
in the original space via the following linear mapping:
P2 = {x ∈ X : xbt ≤ xb−k,t−1 ∀b ∈ B  b ≥ k + 1, ∀t ∈ T } S
Inbt S
= Yˆb,t−1 S
− Znb,t−1 ∀b ∈ B, n ∈ Nb , t ∈ T
we wish to show that P2 ⊆ P1 . S−
Inbt = S
Znbt − S
Znb,t−1 ∀b ∈ B, n ∈ Nb , t ∈ T
Proof. We show that if xbt satisfies (12), then xbt satisfies (11). For S
Ynb2t ∀b ∈ B, n ∈ Nb , t ∈ T
= YˆbtS
each t ∈ T : S
Ynb1t =0 ∀b ∈ B, n ∈ Nb , t ∈ T
|B |
 
k |
B |−k |B |
 |B |

(xbt − xb,t−1 ) = xbt − xb,t−1 + xbt − xb,t−1
S
Ynb3t = 0 ∀b ∈ B, n ∈ Nb , t ∈ T
b=k+1 b=|B |−k+1
b=1 b=1 b=1
That is, substituting into (S ) the expressions on the right-hand-

k |
B |−k |B |
 |B |
 side of the mapping for the variables listed on the left-hand side
≤ xbt − xb,t−1 + xb−k,t−1 − xb,t−1 results in true statements for each relevant set of constraints. (Con-
b=1 b=1 b=k+1 b=|B |−k+1 straints containing only variables not involved in the mapping re-

k |
B |−k |
B |−k |B |
 main unchanged.) This implies that the solution involving YˆbtS and
S is feasible, and therefore valid, for (S ). The same type of sub-
= xbt − xb,t−1 + xb,t−1 − xb,t−1 Znbt
b=1 b=1 b=1 b=|B |−k+1 stitution, and the correct interpretation of the new variables Yˆ S bt
S yields the same objective function value as in the original
and Znbt

k |B |

= xbt − xb,t−1 formulation.
b=1 b=|B |−k+1
References

k
≤ xbt
Alford, C. (2006). Optimization in underground mine design Ph.D. thesis. University of
b=1 Melbourne.
≤k Alford, C. (2007). Optimisation in underground mining. Handbook of Operations Re-
search in Natural Resources.
 Alonso-Ayuso, A., Carvallo, F., Escudero, L. F., Guignard, M., Pi, J., Puranmalka, R.,
& Weintraub, A. (2014). Medium range optimization of copper extraction plan-
Combining the first and last expressions implies that ning under uncertainty in future copper prices. European Journal of Operational
|B| Research, 233(3), 711–726.
(x − xb,t−1 ) ≤ k. The reformulation of the special knap-
b=1 bt AMPL Optimization LLC (2014). Version 20140908. www.ampl.com.
sack constraints is done not only to enable the model to be more Araneda, O. (2015). Opportunities and challenges of the transition from an open pit
to an underground operation in the Chuquicamata mine. In Proceedings of the
easily solved within a framework the OMP Solver can handle, but
mine planning 2015. Antofagasta, Chile.
also to improve the upper bound. For (S ), the reformulation does Bakhtavar, E., Shahriar, K., & Oraee, K. (2008). A model for determining the optimal
not improve the LP solution time over that obtained with the transition depth over from open-pit to underground mining. In Proceedings of
original model when using the OMP Solver for both formulations; the fifth international conference and exhibition on mass mining.
Bienstock, D., & Zuckerberg, M. (2010). Solving LP relaxations of large-scale prece-
however, for the scenarios we test, the LP bound for (S ) improves dence constrained problems. Integer Programming and Combinatorial Optimiza-
with the reformulation by approximately 10 percent. tion, 6080(1), 1–14.
 B. King et al. / European Journal of Operational Research 257 (2017) 297–309 309

Bley, A., Boland, N., Froyland, G., & Zuckerberg, M. (2012). Solving mixed integer Newman, A., & Kuchta, M. (2007). Using aggregation to optimize long-term produc-
nonlinear programming problems for mine production planning with stock- tion planning at an underground mine. European Journal of Operational Research,
piling. Optimization Online, Preprint, 1–30. https://scholar.google.com/scholar? 176(2), 1205–1217.
cluster=4004785229956384832&hl=en&as_sdt=0,6. Newman, A., Yano, C. A., & Rubio, E. (2013). Mining above and below gorund: Tim-
Brazil, M., & Thomas, D. A. (2007). Network optimization for the design of under- ing the transition. IIE Transactions, 45(8), 865–882.
ground mines. Networks, 49(1), 40–50. Osanlo, M., Gholamnejed, J., & Karimi, B. (2008). Long-term open pit mine produc-
Brickey, A. (2015). Underground production scheduling optimization with ventilation tion planning: a review of models and algorithms. International Journal of Min-
constraints Ph.D thesis. Colorado School of Mines. ing, Reclamation, and Environment, 22(1), 3–35.
Carlyle, M., & Eaves, C. (2001). Underground planning at Stillwater Mining Company. O’Sullivan, D. (2013). An optimization-based decomposition heuristic for solving com-
Interfaces, 31(4), 50–60. plex underground mine scheduling problems Ph.D thesis. Colorado School of
Chicoisne, R., Espinoza, D., Goycoolea, M., Moreno, E., & Rubio, E. (2012). A new Mines.
algorithm for the open-pit mine production scheduling problem. Operations Re- O’Sullivan, D., & Newman, A. (2015). Optimization-based heuristics for underground
search, 60(3), 517–528. mine scheduling. European Journal of Operational Research, 241(1), 248–259.
Crawford, J., & Hustrulid, W. (1979). Open pit mine planning and design (1st). Society Qinglin, C., Stillborg, B., & Li, C. (1996). Optimisation of underground mining meth-
of Mining Engineers. ods using grey theory and neural networks. Mine Planning and Equipment Selec-
Epstein, R., Goic, M., Weintraub, A., Catalàn, J., Santibáñez, P., Urrutia, R., tion Conference. Rotterdam: Balkema.
et al. (2014). Optimizing long-term production plans in underground and Ramazan, S. (2007). The new fundamental tree algorithm for production scheduling
open-pit copper mines. Operations Research, 60(1), 4–17. of open pit mines. European Journal of Operational Research, 177(2), 1153–1166.
Finch, A. (2012). Open pit to underground (pp. 88–89). International Mining. Shishvan, M. S., & Sattarvand, J. (2015). Long term production planning of open
Gershon, M. E. (1983). Mine scheduling optimization with mixed integer program- pit mines be any colony optimization. European Journal of Operational Research,
ming. (Littleton, Colo.). Minerals Engineering, 35(4), 351–354. 240(3), 825–836.
IBM CPLEX Optimizer (2013). Version 12.4.0.0. http://www-01.ibm.com/software/ Souza, J. F., Coelho, I. M., Ribas, S., Santos, H. G., & Merschmann, L. H. C. (2010). A
commerce/optimization/cplex-optimizer/. hybrid heuristic algorithm for the open-pit-mining operational planning prob-
Johnson, T. (1968). Optimum open pit mine production scheduling. University of Cali- lem. European Journal of Operational Research, 207(2), 1041–1051.
fornia, Berkeley Technical report. Topal, E., & Ramazan, S. (2010). A new MIP model for mine equipment schedul-
Lambert, B., Brickey, A., Newman, A., & Eurek, K. (2014). Open pit sequencing for- ing by minimizing maintenance cost. European Journal of Operational Research,
mulations: A tutorial. Interfaces, 62(3), 127–142. 207(2), 1065–1071.
Lerchs, H., & Grossman, I. (1965). Optimum design of open-pit mines. Canadian In- Trout, L. (1997). Formulation and application of new underground mine scheduling
stitute of Mining Bulletin, 58, 47–54. models Ph.D thesis. The University of Queensland.
Martinez, M., & Newman, A. (2011). A solution approach for optimizing long- and Underwood, R., & Tolwinski, B. (1998). A mathematical programming viewpoint
short-term production scheduling at LKABs Kiruna mine. European Journal of for solving the ultimate pit problem. European Journal of Operational Research,
Operational Research, 211(1), 184–197. 107(1), 96–107.
Minemax (2013). Mine planning and scheduling solutions. www.minemax.com. Whittle Consulting (2013). Whittle consulting global optimization software. http:
Moreno, E., Rezakhah, M., Newman, A., & Ferreira, F. (2016). Linear models for stock- //www.whittleconsulting.com.au/.
piling in open-pit mining production scheduling problems. Working paper
Muñoz, G., Espinoza, D., Goycoolea, M., Moreno, E., Queyranne, M., &
Rivera, O. (2015). Production scheduling for strategic open pit mine planning,
part I: A mixed integer programming approach. Working paper