Anderson Graduate School of Management ...

9 downloads 0 Views 663KB Size Report
Apr 1, 2010 - the conversion of Chiquicamata, Codelco's largest open pit mine, to underground ...... Underground planning at Stillwater Mining Company.
Anderson Graduate School of Management – Decisions, Operations, and Technology Management UC Los Angeles Peer Reviewed Title: Optimizing Long�Term Production Plans in Underground and Open Pit Copper Mines Author: Epstein, R., University of Chile Goic, M., University of Chile Weintraub, A., University of Chile Catalán, J., Fundación para la Transferencia Tecnológica Santibáñez, P., Fundación para la Transferencia Tecnológica Urrutia, R., Fundación para la Transferencia Tecnológica Cancino, R., CODELCO Chile Gaete, S., CODELCO Chile Aguayo, A., CODELCO Chile Caro, F., UCLA Anderson School Publication Date: 04-01-2010 Publication Info: Decisions, Operations, and Technology Management, Anderson Graduate School of Management, UC Los Angeles Permalink: http://escholarship.org/uc/item/475353w1 Abstract: We present a methodology for long�term mine planning based on a general capacitated multicommodity network flow problem formulated as a comprehensive math programming model. The planning methodology considers underground and open pit ore deposits sharing multiple downstream processing plants over a long horizon. The purpose of the model is to optimize several mines in an integrated fashion, though the real size instances are hard to solve due to the combinatorial nature of the problem. We tackle this by solving the relaxed problem of a tight linear formulation and we round the resulting near�integer solution with a customized procedure. The outcome of this project has been implemented at Codelco, the largest copper producer in the world. Since 2001, the system is used on a regular basis and has increased the net present value of the production plan for a single mine by 5%. Moreover, the integrated approach for multiple mines provided an additional increase of 3%. The system has allowed the planners to evaluate more scenarios. In particular, the model was used to study the option of delaying by four years the conversion of Chiquicamata, Codelco's largest open pit mine, to underground operations.

eScholarship provides open access, scholarly publishing services to the University of California and delivers a dynamic research platform to scholars worldwide.

Optimizing Long‐Term Production Plans in Underground and  Open Pit Copper Mines    Rafael Epstein, Marcel Goic, Andrés Weintraub  Industrial Engineering Department, University of Chile, Casilla 2777, Santiago, Chile  {repstein, mgoic, aweintra}@dii.uchile.cl 

  Jaime Catalán, Pablo Santibáñez, Rodolfo Urrutia  Fundación para la Transferencia Tecnológica, Avda. Beauchef 993, Santiago, Chile  {jcatalan, psantiba, rurrutia}@dii.uchile.cl 

  Raúl Cancino, Sergio Gaete, Augusto Aguayo  CODELCO Chile, División Norte – El Teniente – División Andina  {rcancino, sgaete, aaguayo}@codelco.cl 

  Felipe Caro1  UCLA Anderson School of Management, Los Angeles, CA 90095  [email protected] 

  April 22, 2010    ABSTRACT  We present a methodology for long‐term mine planning based on a general capacitated multi‐ commodity network flow problem formulated as a comprehensive math programming model.  The planning methodology considers underground and open pit ore deposits sharing multiple  downstream  processing  plants  over  a  long  horizon.  The  purpose  of  the  model  is  to  optimize  several mines in an integrated fashion, though the real size instances are hard to solve due to  the  combinatorial  nature  of  the  problem.  We  tackle  this  by  solving  the  relaxed  problem  of  a  tight  linear  formulation  and  we  round  the  resulting  near‐integer  solution  with  a  customized  procedure. The outcome of this project has been implemented at Codelco, the largest copper  producer in the world. Since 2001, the system is used on a regular basis and has increased the  net  present  value  of  the  production  plan  for  a  single  mine  by  5%.  Moreover,  the  integrated  approach for multiple mines provided an additional increase of 3%. The system has allowed the  planners to evaluate more scenarios. In particular, the model was used to study the option of  delaying  by  four  years  the  conversion  of  Chiquicamata,  Codelco's  largest  open  pit  mine,  to  underground operations.   

1

 Corresponding author. 



1. INTRODUCTION  Basically there are two major forms of mining: underground and open pit. For both, the use of  models to support extraction planning has a long practice in the mining industry. Traditionally,  systems based on average values have been used. A fundamental concept has been that of a  cutoff grade (Lane 1988) based on the idea of finding which is the minimum grade in a block  (i.e., the basic geological planning unit) which is profitable to mine.   Given the advances in software and hardware, lately the idea of modeling the mine in  finer detail and reality as is possible using mixed‐integer programming (MIP) models has gained  strength  (Alford  et  al.  2007,  Cacceta  2007).  This  approach  is  far  more  powerful,  as  it  allows  representing the mine in complete detail, where the block of rock is the minimum element that  can  be  defined.  In  this  form,  a  MIP  can  represent  the  extraction  problem  in  a  more  exact,  realistic  way.  These  models  can  incorporate  in  a  natural  form  constraints  related  to  geomechanical, economic or extraction sequencing conditions. The size of such models though  makes these formulations difficult to solve, and the literature so far presents few applications  of  large‐scale  mines.  This  paper  presents  and  solves  a  MIP  model  that  has  been  successfully  used  in  Chilean  copper  mines  by  state‐owned  Codelco,  the  largest  copper  enterprise  in  the  world, for both underground and open pit mines.    The paper has several novelties:  ƒ It integrates the mine extraction with downstream processes, showing the advantages of  integrating  planning  along  the  whole  production  chain,  until  the  final  product,  refined  copper  is  obtained  to  be  shipped.  Traditionally,  mining  and  plant  processes  have  been  solved separately, iterating between the two to match supply and demand.  ƒ It considers, as is the case in Codelco, the interaction of multiple mines, multiple products  and  multiple  downstream  processes.  In  particular,  it  incorporates  the  notion  of  prime  material  competing  for  limited  processing  capacity  as  the  mineral  goes  along  the  chain.  Traditionally, planning has been carried out not considering this multiplicity.  ƒ It develops a model that generalizes the concept of underground and open pit mining, as  part of a capacitated network representation.  This allows representing any form of mining  as part of a process in the production chain, and in particular the change from open pit to  underground becomes natural to introduce.  ƒ The  network  representation  allows  representing  directly:  a)  constraints  related  to  maximum  and  minimal  rates  of  extraction,  b)  smoothness  of  extraction  and  correct  sequencing  of  extraction  points  in  underground  mines,  c)  correct  breaking  of  rocks,  geomechanical  constraints  and  possible  interactions  between  sectors,  and  d)  gradually  sloping walls in open pit mining to avoid sliding of waste material. Also, in the network, the  use and investments in machinery, equipment and transportation are modeled in a natural  2 

way.  These  models  are  very  large  MIP's,  impossible  to  solve  using  standard  commercial  codes.  Heuristic  techniques  were  developed  which  allowed  solving  these  problems  with  relatively small error bounds.    The model to be shown in the next sections has been implemented at Codelco. First, a  model was developed and used for extraction and investment planning on a 25 year horizon at  El Teniente, a very large underground copper mine. Its use supported decisions on the timing,  sequencing of blocks to be mined, as well as deciding which blocks did not yield enough value  to be mined.  Improvements due to the use of the model led to an increase in value of the mine  of over $100 million.   Then,  the  model  was  expanded  to  incorporate  open  pit  mines.  It  is  being  used  for  planning  in  the  northern  groups  of  mines  in  Chile.  In  this  case,  several  mines  share  common  plant  installations,  which  include  processes  of  flotation,  leaching,  bioleaching  and  low  grade  sulfates. The need to have one model integrating into the production chain both underground  and open pit mines is due to increasingly common case of  mixed mines, where both types of  mines  feed  material  to  the  downstream  processes,  or  an  open  pit  mine  is  converted  to  an  underground one. We show specific cases of such mines.  An important result of the use of the model is in guiding the allocation of mineral to the  processing  plants,  taking  advantage  of  the  comparative  advantages  of  mineral  from  different  mines to be processed in the different plants. In this case, many instances have been run of the  model to test multiple scenarios. The solution of these runs is coordinated with mining analysis  to  determine  the  timing  and  sequencing  of  the  operation  in  order  to  be  integrated  with  downstream plant processes. We show an important impact of the introduction of our model  into  the  decision  process,  with  gains  between  5%  and  8%  in  net  present  value  confirmed  by  three independent case studies here described.  The reminder of the paper has the following structure. In the next section we present a  brief literature review. In section 3 we describe the problem and we provide some background  on mining operations and the planning process at Codelco. In section 4 we present our planning  methodology  based  on  a  general  capacitated  multi‐commodity  network  flow  problem.  In  section  5  we  introduce  the  customized  rounding  heuristics  used  to  generate  approximate  solutions.  In  section  6  we  provide  an  overview  of  the  model’s  implementation  and  impact  at  Codelco. In the final section we conclude. The theoretical proofs are given in the appendix.    2. LITERATURE REVIEW   Most optimization models available in the literature have been developed for open pit mines  and they solve only a partial version of the long‐term planning problem (Newman et al. 2009).  Historically, the problem has been divided in two, the ultimate pit problem, which determines  3 

the  economic  mineral  reserve  in  an  open  pit  mine,  and  the  production‐scheduling  problem,  which  considers not only which blocks to remove, but also when to remove these blocks.   There  are  two  main  approaches  to  solve  the  ultimate  pit  problem.  One  is  based  on  cutoff grades and the other one is based on Operations Research (OR) techniques. The optimal  cutoff  grade  methodology  was  made  popular  through  the  work  of  Lane  (1988).  Fytas  et  al.  (1987),  Kim  and  Zhao  (1994),  and  Poniewierski  et  al  (2003)  have  extensively  discussed  this  process. The basic premise of this approach is that one can determine a cutoff grade policy to  maximize  net  present  value  subject  to  capacity  constraints,  with  higher  cutoff  grades  in  the  initial  years  of  the  project  leading  to  higher  overall  net  present  value.  It  has  important  operational  advantages  and  it  is  embedded  in  the  background  of  all  mining  practitioners.  However,  the  assumption  of  a  fixed  cutoff  grade,  which  depends  on  an  arbitrary  delineation  between ore and waste, generates suboptimal solutions because it ignores that the value of a  block  is  not  inherent  to  the  block  but  rather  depends  on  the  interaction  with  the  rest  of  the  blocks and the capacity of the downstream processes.  The  use  of  OR  techniques  to  solve  the  ultimate  pit  problem  started  with  the  classic  “moving cone” heuristic, see for example Kim (1978) or Laurich (1990). This approach takes a  block  as  a  reference  point  and  expandes  the  pit  upwards  according  to  pit  slope  rules.  This  solution  can  be  suboptimal  but  it  is  intuitively  appealing.  Among  the  algorithms  with  a  convergence  proof  to  optimality,  historically  the  most  important  are  Lerchs  and  Grossmann  (1965)  and  Picard  (1976).  The  first  one  is  based  on  graph  theory  but  it  can  be  shown  that  its  structure  is  very  similar  to  the  dual  simplex  algorithm,  see  Underwood  and  Tolwinski  (1998).  The  second  paper  shows  the  equivalence  of  the  ultimate  pit  problem  and  the  well‐studied  problem of finding a maximum closure in a graph. Therefore, the ultimate pit problem reduces  to a maximum flow problem that can be solved, for example, with a push‐relabel procedure. A  comparison of both procedures together with a complete literature review of the ultimate pit  problem  can  be  found  in  Hochbaum  and  Chen  (2000).  Other  relevant  papers  in  this  vein  are  Wright  (1989),  who  suggests  a  dynamic  programming  approach,  Frimpong  et  al.  (2002),  who  use  neural  networks,  and  Jalali  et  al.  (2006),  who  propose  the  use  of  Markov  chains.  Even  though  these  algorithms  are  efficient,  they  do  not  consider  what  happens  with  the  mineral  beyond  the  extraction  stage  and,  for  instance,  do  not  provide  a  production  schedule  that  considers all the downstream capacities.    OR  techniques  have  also  been  used  to  solve  the  production‐scheduling  problem  for  a  single  mine.  Among  the  papers  based  on  optimization  methods,  Gershon  (1983)  uses  mixed  integer programming, Dagdelen and Johnson (1986) suggest a Lagrangian relaxation, Caccetta  and Hill (2003) use branch and cut, while Tolwinski and Underwood (1996) propose a dynamic  programming approach.  Given the complexity of the problem, several papers have resorted to  heuristic  methods  such  as  simulation  (Fytas  et  al.  1993,  Erarslan  and  Celebi  2001),  neural  4 

networks  (Denby  et  al.  1991),  simulated  annealing  (Kumral  and  Dowd  2005)  and  genetic  algorithms  (Samanta  et  al.  2005).    For  other  approaches  that  combine  optimization  and  heuristic methods see also Sevim and Lei (1988), Sarin and West‐Hansen (2005), Busnach et al.  (1985),  Dowd  and  Onur  (1993),  Elevli  (1995),  Wang  and  Sun  (2001),  Smith  et  al.  (2003),  Ramazan (2007), and Chicoisne et al. (2009).   The literature on underground mining is more recent, partially due to the complicated  nature  of  underground  operations.  No  parallel  to  the  Lerchs‐Grossman  algorithm  exists  for  underground mines, so the early optimization work does not rely on network‐based methods.  Strategic decisions in underground mining may consider how to geographically position various  facilities,  e.g.,  mills,  on  the  mine  site.  Qinglin  et  al.  (1996)  use  neural  networks  based  on  a  deposit’s  geotechnical  and  economic  factors  to  optimally  select  an  underground  mining  method.  Yun  et  al.  (1990)  use  a  genetic  algorithm  to  determine  the  number  and  spacing  of  openings  given  restrictions  on  their  relative  placement.  Alford  (1995)  describes  the  floating  stope  method  as  a  tool  for  analysis  of  mineral  reserves  and  stope  geometry.  Barbaro  and  Ramani (1986) formulate a mixed integer programming model that can be used to determine  whether  to  produce  in  a  given  time  period.  Kuchta  et  al.  (2004)  and  Kuchta  et  al.  (2003)  describe  the  implementation  of  a  production  scheduling  mixed  integer  program  at  LKAB's  Kiruna  underground  mine.  The  objective  is  to  minimize  the  deviation  from  a  pre‐specified  demand  target  per  period,  which  contrasts  with  our  model  that  maximizes  net  present  value  (NPV). The sizes of the instances are not comparable either but there are some similarities in  terms the methods used to reduce the size of the model.    As noted in Hochbaum and Chen (2000), a comprehensive model of mining operations  has seldom been addressed in the literature. Even fewer implementations of such models have  been reported in practice. The work by Carlyle and Eaves (2001) is a noteworthy attempt. Their  mixed  integer  model  was  developed  for  an  underground  platinum  and  palladium  mine  and  it  considers several planning decisions, but it was applied to only one sector of the mine and near  optimal integer solutions were obtained with a commercial software package. In our case, we  consider  all the  sectors simultaneously  and  we even  allow  for  multiple mines.  In  fact,  we are  unaware  of other  implementations that  solve  an  integrated  long‐term planning  problem  with  the level of detail considered here.     3. PROBLEM DESCRIPTION  We now provide background information and give more details of the problem tackled in the  paper.  In  section  3.1  we  introduce the  long‐term  planning  problem  in  the  mining  industry.  In  section  3.2  we  provide  some  background  on  underground  and  open  pit  mining.  Finally,  in  section 3.3 we describe the planning process at Codelco and how our methodology fits in it.    5 

3.1 Long‐Term Planning   A  long‐term  mining  plan  consists  of  a  detailed  specification  of  each  aspect  related  to  investments  and  production  in  the  mines  and  plants.  The  importance of  the  plan  stems  from  the fact that most decisions are final and cannot be reversed, meaning that they can drastically  affect the future of the mine. For example, poor capacity choices at the concentration plant or  selecting an inappropriate extraction sequence can trigger a significant yield decrease.    Since most mines have a long life cycle, the planning horizon is set somewhere between 10  and 30 years, with shorter periods at the beginning when more detail is required and deviations  in the production have a major impact in the results. The plan not only must be economically  profitable, but it also must be technically impeccable. This means, among other things, that is  has  to  comply  with  the  geotechnical  rules  that  govern  the  mine,  it  must  consider  constraints  due  to  the  chemical  processes  at  the  plants,  it  must  meet  environmental  regulations,  and  it  should guarantee that downstream capacities will not be exceeded.    The  long‐term  planner  must  solve  three  main  problems  taking  into  account  the  considerations  mentioned  above:  (i)  Investment:  determine  the  selection  and  timing  of  investments;  (ii)  Extraction:  determine  the  production  in  the  mine;  and  (iii)  Processing:  determine the operation of the plants. Solving each one of these problems, even separately, is  already  a  complex  task.  However,  the  need  of  an  integrated  planning  approach  cannot  be  dismissed.  For  instance,  the  real  value  of  an  investment  project  is  only  appreciated  once  its  coherence is verified with respect to the extraction and processing decisions. Our methodology  explicitly solves problems (ii) and (iii), i.e., extraction and processing. The investment decisions  are usually restricted to a small set of options and can be evaluated with our model by running  different scenarios.    3.2 Mining Operations  3.2.1 Underground mining  We here present a brief description of the operations at an underground mine. The first step is  the selection of reserves. This consists in defining the boundaries of the orebody, delimiting the  material  to  be  removed.  This  is  done  according  to  economic  criteria  that  determine  what  is  profitable based on the grade distribution. Large underground mines are typically subdivided in  sectors by design, and the selection of reserves is done for each sector.      For  planning  and  operational  purposes,  the  geological  model  of  the  whole  deposit  is  expressed  through  a  block  model.  Each  block  is  uniquely  identified  together  with  its  geologic  characteristic,  in  particular  the  ore  grades.  These  values  are  estimated  using  geostatistic  procedures (kriging). A column (or drawpoint) is the vertical aggregation of blocks, and a mine  or sector corresponds to a set of neighboring columns. 



  Once the reserves are identified, the following step is to program the ore extraction on a  time scale taking into account geological restrictions and the downstream capacities. There are  several  mining  methods  for  underground  mines  (Newman  et  al.  2009).  In  our  model,  we  consider  the  block  caving  method  which  is  the  prevalent  underground  mining  method  at  Codelco. In a nutshell, the block caving method consists in creating a void at each drawpoint so  that the rock breaks and falls due to gravity and the weight of the overburden material. For this  to happen, the rock extraction pattern must follow specific rules: (i) the columns have to enter  production in a particular sequence to generate a "wave" that breaks the rock; (ii) the wave has  to advance smoothly, which requires regularity in heights among neighboring columns; and (iii)  at  each  drawpoint  there  is  maximum  extraction  rate  to  prevent  the  roof  from  collapsing  and  there is a minimum number of blocks that must be removed in order to avoid solid pillars.    Figure 1 shows the typical flow in an underground mine. The broken rock is removed from  the drawpoints which are arranged along parallel crosscuts. Each crosscut has enough space for  an  LHD  machine  that  hauls  the  material  to  a  dumping  point.2  The  material  is  then  funneled  through ore passes and reaches the internal crusher, where the rock is crushed to a size that  can  be  transported  by train  to  the  processing stages  outside  the mine.  For  our  purposes,  the  flow  finishes  at  the  concentration  plants  where  the  profitable  minerals  are  obtained  by  eliminating  waste.  Each  one  of  these  intermediate  processes  can  be  characterized  through  technical coefficients such as its capacity and variable operational cost.  

  Figure 1: Description of underground mining operations.  2

 LHD stands for Loading, Hauling & Dumping. 



3.2.2 Open pit mining  Following the discussion of underground mining we review here the operations at an open pit  mine. Open pit mining is conceptually similar to underground mining, except that the orebody  is reached from above, which requires removing plenty of overburden material from the soil. A  preliminary definition of the regions to be mined is done using geological models that take into  account  the  geometry  of  the  pit  and  technical  requirements  (e.g.,  the  maximum  slope  to  prevent the walls of the pit from collapsing). The regions are called expansions which are also  known as pushbacks or phases. In figurative terms, an expansion can be viewed as a "slice" of  the  pit.  Each  expansion  is  subdivided  in  benches  of  a  pre‐determined  height  that  must  be  extracted  in  order,  from  top  to  bottom.  Figure  2  illustrates  the  geometry  of  expansions  and  benches and Figure 3 shows a view from above.   

  Figure 2: Description of open pit mining operations. Within an expansion, benches are  extracted from top to bottom.    E x p a n s ió n 4 8 N E

E x p a n s ió n 4 1 E

E x p a n s ió n 5 5 S

Figure 3: Expansions of an open pit mine, viewed from above.      The bench is the minimal extraction decision unit. When a bench is mined all its material  needs  to  be  removed  from  the  pit  including  waste  without  economic  value.  There  are  three  processes  involved  in  removing  a  bench:    Perforation  and  blasting,  where  the  bench  is  separated from the soil using dynamite; loading, where all the material is loaded into large size  trucks; and transportation, where the material is taken to either a metallurgic process or to a  8 

waste dump. It is worth to mention that perforation and blasting imposes severe constraints on  the  sequence  in  which  expansions  are  extracted.  For  example,  there  is  a  minimum  distance  requirement  between  two  expansions  in  operation  on  the  same  wall  of  the  pit  to  avoid  rock  spillage.    Like  in  underground  mining,  the  sequence  of  downstream  processes  involves  rock  size  reduction and ore recovery. However, there are some differences that are worth mentioning.  First, the existence of stockpiles plays an important role in balancing the plant inflow and they  are  held  for  several  years  before  entering  the  metallurgic  processes.  Second,  a  significant  fraction of the material extracted from the mine is taken to the waste dumps which might be  located  several  miles  away  from  the  pit.  This  introduces  an  important  trade‐off  between  sending low grade ore to the downstream processes or to the costly dumps.    The existence of more than a single ore type, sulfides and oxides, is a key element of an  integrated  plan.  Although  in  the  geological  reserves  of  a  copper  mine  we  typically  find  both  types, in the traditional planning scheme, mining plans are defined considering only one main  metallurgic process for a particular type. Thus, an integrated plan considers different ore types  in alternative concentration lines that would be otherwise discarded or sent to processes with  low yields.    3.3 Codelco's Production Planning Approach  Recently Codelco has migrated from a traditional planning scheme based exclusively on cut‐off  grades to a methodology based on the model presented in this paper. In order to compare the  results we outline the planning approach. Note that overall the steps have not changed. What  has been modified is the way steps d) and e) are carried out, as explained below.  a) The available resources are determined considering the current geologic characteristics of  the deposit and the mine’s history.  b) A  selection  of  resources  is  performed.  This  is  done  considering  Macro  Options  (i.e.  alternative  technologies,  economic  indicators,  management  directives,  etc)  in  order  to  delimit  the  floor,  perimeter  and  ceiling  of  the  resources  to  extract  during  the  planning  horizon. The outcome is known as the Model of Mining Reserve that contains the resources  with positive benefit taking into account the processing costs and the projected recovery.   c) The  phases  and  benches  of  extraction  are  designed  on  the  basis  of  physical  y  geomechanical principles.   d) The Consumption Strategy defines a Mining Plan, which contains for every period: tons of  mineral to extract, productive sectors involved, ore grades and level of uncertainty based  on the drillings and samples.  



e) The  Mining  Plan  is  economically  evaluated,  providing  Profit  Indicators,  and  measures  of  Technical  Risk  and  Vulnerability.  If  the  result  is  not  satisfactory,  the  Planning  Parameters  are modified and the process is repeated from point c).   f) If an external change occurs affecting the Macro Options, these are incorporated into the  model in point b), and the process is repeated from that stage.   g) The final result is a plan that provides the relevant information for operating each sector,  including the starting period, operation rates and termination period.      In the legacy planning approach steps d) and e) where carried out by looking separately at  the different parts of the mine. This scheme allowed evaluating different options by means of  alternative plans, but not in a unified (global) approach, as each part was treated separately. In  addition,  due  to  the  high  complexity  of  the  procedure  and  the  involved  calculations,  several  weeks (if not months) were spent in obtaining one particular plan, preventing the planner from  evaluating more options, going back and forth between the different steps. Using the model in  steps d) and e) allows planning in an automated fast way. It also leads to better solutions as the  whole  system  is  planned  in  an  integrated  fashion,  and  the  planner  can  evaluate  multiple  scenarios and options.    4. MODEL FORMULATION AND DISCUSSION   In  this  section  we  present  a  math  programming  approach  to  optimize  long‐term  production  plans  in  open  pit  and  underground  mines  that  share  downstream  processes.  A  graphical  representation  of  the  different  stages  in  the  model  is  given  in  Figure  4  where  material  flows  from left to right. 

Figure 4: Network flow representation for mining operations where several open pit and  underground mines share downstream processes.  10 

  The network of mining operations begins with the production or rock extraction phase that  takes places at sectors and expansions for underground and open pit mines respectively. This  phase corresponds to the first stage on the left side of Figure 4. Then the extracted rock is fed  into the network of downstream processes that involves size reduction and chemical reactions.  The  last  stage  is  the  concentration  phase  which  yields  the  three  final  products:  copper,  molybdenum  and  arsenic.  The  first  two  have  commercial  value  while  the  third  one  is  a  contaminant on which environmental restrictions are imposed.     In what follows, we present the model with its essential components. When necessary we  use  the  sub‐indices O and U to  denote  the  two  type  of  mines  open  pit  and  underground  respectively. In section 4.1 we first introduce a unified model for the network of downstream  processes. This part of the model is shared by all the mines considered regardless of their type.  In  contrast,  we  present  a  different  model  of  the  production  phase  for  each  type  of  mine  in  order  to  capture  the  unique  characteristic  of  the  respective  mining  methods.  This  is  done  in  sections 4.2 and 4.3 for underground and open pit respectively.     4.1 Unified Model Formulation for Downstream Processes  As  shown  in  Figure  4,  the  downstream  processes  are  modeled  as  a  capacitated  multi‐ commodity network flow problem. The starting nodes  s ∈ S = SU ∪ SO  represent the different  sectors  and  expansions  for  underground  ( SU )  and  open  pit  mines  ( SO )  respectively.  The  intermediate nodes  u , v ∈ V  represent stocking points and processing stages that precede the  concentration plants, which are the final nodes ( F ) in the network. Several different products  flow  through  the  network.  A  given  product  k ∈ K   represents  material  (i.e.,  rock)  with  an  average  rock  size  and  a  grade  of  copper,  molybdenum,  and  arsenic  within  a  certain  range.  Depending on the characteristics of the ore deposits, the product definition might also include  other  relevant  attributes  such  as  the  chemical  nature  (sulfide  or  oxide)  or  certain  geological  properties (e.g., rock hardness). The set of commercial products is denoted COM , whereas the  set of contaminants is denoted CONT .   

The flow of product  k ∈ K  measured in tons and going from node  u to node  v  in period 

t ∈ Τ = {1,..., T }  is denoted fuvkt . As the material flows through the downstream processes, the  average  rock  size  is  reduced  and  the  valuable  minerals  are  separated  from  waste.  Let  TCvkk '   denote the transformation coefficient at node v ∈ V . That is, the amount of output  k '  obtained  for  each  ton  of  input k ,  with ∑ k '∈K TCvkk ' = 1 ∀v ∈ V , k ∈ K .  Let cvkt denote  the  cost  of  processing a ton of product  k at node  v in period t ∈ Τ . The maximum processing capacity at  node  v  in period  t is given by CPvt . For nodes  v ∈ V that can hold inventory, let  yvkt represent  the amount of product  k available at the beginning of period  t , and let  CSvt be the maximum   11 

stock  level  per  period  in  tons  (if  the  node  cannot  hold  inventory,  then  CSvt = 0 ).  If k   is  a  contaminant,  then  the  maximum  amount  that  can  be  released  per  period  is  Ekt   (otherwise, 

Ekt = ∞ ). We assume that the firm is a price‐taker and that all the production of product k in  period  t  can be sold at the market price  pkt (clearly,  pkt = 0 for non‐commercial products).   

The link of the downstream processes with the upstream production phase is established 

through  the  auxiliary  variables ProdCostst and xskt .  The  former  denotes  the  total  production  cost,  while  the  latter  represents  the  tons  of  product k produced  in  sector/expansion s in  period t . The mathematical formulation of the downstream processes is the following:   



T

max

∑∑ ∑

t =1 k∈K u∈V ∪ S

βt ⎜ ∑



pk 'tTCvkk ' fuvkt −

⎝ v∈F k '∈COM

⎞ T cvkt fuvkt ⎟ − ∑∑ βt ProdCostst v∈V ∪ F ⎠ t =1 s∈S



(1)

subject to

∑ f = ∑ f

xskt =

∑ ∑ TC

u∈V ∪ S k∈K

f

vkk ' uvkt

+ yvk 't −1

∑ ∑f ∑y

u∈V ∪ S k∈K

∑ ∑ ∑ TC

u∈V ∪ F

svkt

vuk ' t

+ yvk 't

∀s ∈ S , k ∈ K , t ∈ Τ

(2)

∀v ∈ V , k ' ∈ K , t ∈ Τ

(3)

uvkt

≤ CPvt

∀v ∈ V ∪ F , t ∈ Τ

(4)

vkt

≤ CSvt

∀v ∈ V , t ∈ Τ

(5)

≤ Ek 't

∀k ' ∈ CONT , t ∈ Τ

(6)

k∈K

u∈V ∪ S v∈F k∈K

v∈V ∪ F

f

vkk ' uvkt

fuvkt , yvkt ≥ 0

∀u, v ∈ V , k ∈ K , t ∈ Τ. (7)

    The objective of the model is to maximize the net present value of the mines over the next  T periods.  The  first  term in  the  objective  function  (1)  corresponds  to  the  discounted  benefits  obtained from selling the final products. The second and third terms represent the discounted  processing and production costs respectively (here the parameter  β t < 1  is the discount factor).  Constraints  (2)  feed  the  production  from  each  sector/expansion  into  the  network  of  downstream processes. Constraints (3) ensure the flow conservation and inventory balance at  each  node.  Constraints  (4)  and  (5)  impose  maximum  processing  capacities  and  stock  levels  respectively.  Finally,  constraints  (6)  limit  the  release  of  pollutants  at  the  final  stage,  and  constraints (7) are the usual non‐negativity conditions on flows and inventory levels.    4.2 Underground Extraction  To model the production phase for underground mines, we assume that the mining method is  block caving, as explained in section 3.2.1. Let the index  i ∈ I ( s)  denote a drawpoint or column  12 

in  sector s ∈ SU   and  let  the  index  n ∈ N (i )   denote  the  blocks  in  that  column  which  are  numbered 0 to N (i ) − 1 from bottom to top. The main decision in the production phase is when  to  remove  each  block.    For  that,  we  introduce  the  binary  variable zint ,  that  equals  one  if  the  extraction  of  block i   in  column n   is  initiated  in  period t ,  and  the  continuous  variable eint ,  which represents the fraction of the block that is removed in each period.    

As part of the model's input, we need the following sets and parameters. Let TON ink be the 

amount (tons) of product k contained in block (i, n) . Let EMIN (i ) denote the minimum set of  blocks  of  column i   that  must  be  removed  once  production  at  that  column  is  initiated.  Let 

IU = ∪ s∈SU I ( s )   be  the  set  of  drawpoints/columns  in  all  the  underground  mines  considered.  Let Δ denote  the  maximum  number  of  periods  a  drawpoint  can  remain  open.  Let  α in   be  the  height of block (i, n) and let the variable hit denote the height of column i  in period t , which is  given  by  the  blocks  that  have  been  removed  so  far.  Let δ s   denote  the  maximum  height  differential  allowed  between  two  neighboring  columns  in  sector s ,  and  let IJ be  the  set  of  column pairs that are close enough to be considered neighbors.    

The  extraction  precedence  relationship  among  columns  is  given  by  the  set SECU , 

where (i, j ) ∈ SECU   if  the  extraction  of  column i   must  begin  before  it  takes  places  at  column j . Initiating the production at column i  incurs a fixed cost ai  and removing block (i, n)   takes  a  minimum  of γ in   days.  There  is  also  a  variable  cost bst per  ton  of  rock  extracted  from  sector s in period t . Let Dt denote the duration of period t  measured in days, and for sector s ,  the maximum (minimum) incorporated area and extracted tons allowed in that period are given  by  MAX_AREAst   ( MIN_AREAst )  and  MAX_EXTst   ( MIN_EXTst )  respectively.  Finally,  let 

WIN ( s ) denote  the  time  window  in  which  sector s can  be  in  production.  The  following  constraints complete the model for underground mines.   

xskt =

ProdCostst =

∑z t∈Τ t

int

∑e g =1

ing

∑ ∑ TON

i∈I ( s ) n∈N ( i )

∑ az

i∈I ( s )

i i 0t

≤1

e

ink int

+ bst ∑ xskt k∈K

∀s ∈ SU , k ∈ K , t ∈ Τ

(8)

∀s ∈ SU , t ∈ Τ

(9)

∀i ∈ IU , n ∈ N (i ) t

≤ ∑ zing g =1

13 

(10)

∀i ∈ IU , n ∈ N (i ), t ∈ Τ (11)

T

T

∑ ∑

t =1 n∈EMIN ( i )

eint ≥ ∑ EMIN (i ) zi 0t t =1 t

t

∑z g =1

in +1 g

t

∑z g =1

zi 0t +

j0g

∑e



n∈N ( i )

ing

≤ ∑ ei 0 g

∀(i, j ) ∈ SECU , t ∈ Τ

≤1

∀i ∈ IU , n ∈ N (i ), t ∈ Τ (15)

g =1

γ in eint ≤ DPt hit = hit −1 +

∑α

n∈N ( i )

in

hit − h jt ≤ δ s MIN_AREAst ≤



i∈I ( s )

eint

(14)

∀i ∈ IU , t ∈ Τ

(16)

∀i ∈ IU , t ∈ Τ

(17)

∀i, j ∈ IJ , s ∈ SU , t ∈ Τ (18)

AREAi zi 0t ≤ MAX_AREAst

MIN_EXTst ≤ ∑ xskt ≤ MAX_EXTst k∈K

∑z

(12)

∀i ∈ IU , n ∈ N (i ), t ∈ Τ (13)

g =1 t

T

g =t+Δ

≤ ∑ eing

∀i ∈ IU

∀s ∈ SU , t ∈ Τ

(19)

∀s ∈ SU , t ∈ Τ

(20)

=0

∀s ∈ SU , t ∈ Τ \ WIN ( s ) (21)

zint ∈ {0,1}, eint , hit ≥ 0

∀i ∈ IU , n ∈ N (i ), t ∈ Τ (22)

i∈I ( s )

i 0t

    Constraints (8) define the production per sector, product and period. Constraints (9) define  the production cost for underground mines which has a fixed cost component for each column  that starts production and a variable cost per ton of material removed. Constraints (10) ensure  that  each  block  is  removed  at  most  once.  Constraints  (11)  establish  the  logical  link  between  variables eint  and zint  so that the extraction of a block cannot occur in a period prior to when it is  initiated. Constraints (12) ensure that the minimum set of blocks is removed if a column enters  production. Constraints (13) require the extraction of block n  to finish before the extraction of  block n + 1   can  start  at  each  drawpoint.  Constraints  (14)  restrict  the  order  in  which  the  drawpoints  can  enter  production.3  Constraints  (15)  impose  the  maximum  duration  of  a  drawpoint. Constraints (16) limit the extraction rate of each column. Constraints (17) define the  height of the column in each period and constraints (18) guarantee that neighboring columns  have  similar  heights  so  that  the  rock  breaks  in  a  smooth  fashion.  Constraints  (19)  and  (20)  impose  bounds  on  the  area  incorporated  and  the  rock  extracted  per  period  respectively.  Constraints (21) ensure that a sector is not extracted outside its time window, and constraints  (22) restrict the domain of the decision variables.  3

 Note that the left hand side of constraints (11), (13) and (14) is formulated as a cumulative sum over time, which  helps to tighten the LP relaxation. 

14 

 

In  the  model,  we  need  constraints  (11)  and  (13)  together  with  the  binary  requirement 

on zint to ensure that the blocks of a column are extracted one by one from the bottom to the  top. However, if the blocks at the bottom have the highest grades, then the binary requirement  can be relaxed since the sole fact that this is a maximization problem will force the extraction to  occur sequentially in the desired order. Hence, in that case the LP relaxation provides a near‐ integer  solution.  This  is  an  important  observation  that  is  relevant  to  our  solution  method.  Therefore, we state it formally in the following lemma.    Lemma  1.    If  the  blocks n ∈ N (i )   in  column i   are  identical  except  for  the  grade  of  the  commercial  minerals,  which  are  decreasing,  then  the  optimal  solution  to  the  LP  relaxation  satisfies  zint ∈ {0,1}  for all blocks except the first and last one extracted in period t .    4.3 Open Pit Extraction  In contrast with underground mines, the mining method for open pit excavations allows for a  more aggregated model since there is no need of getting into the details at the block level. Here  the index i ∈ I ( s )  represents a drawpoint or bench in expansion s ∈ SO , which are numbered 1  to I ( s )   from  top  to  bottom.  The  decision  variables  are: zit ∈ {0,1} ,  which  equals  one  if  the  extraction  of  bench i   is  initiated  in  period t ;  and eit ≥ 0 ,  which  is  a  continuous  variable  that  represents  the  amount  (tons)  extracted  per  period.  Let  I O = ∪ s∈SO I ( s )   be  the  set  of  drawpoints/benches  in  all  the  open  pit  mines  under  consideration.  Let S (m) be  the  set  of  all  expansions in the open pit mine m ∈ M O . Let SECO  denote the set of all pairs (i, j )  such that  the  extraction  at  bench i   must  be initiated  before  it  takes  place  at  bench j .  Note  that  these  could be benches from different expansions. All the other parameters are the same as in the  underground case. The following constraints complete the production plan optimization model  for open pit mines.   

xskt = ProdCostst =

∑z t∈Τ t

it

∑ TON

i∈I ( s )

∑ az

i∈I ( s )

i it

≤1 t

∑ eig ≤ ∑ zig g =1

g =1

15 

e

ik it

+ bst ∑ xskt k∈K

∀s ∈ SO , k ∈ K , t ∈ Τ

(23)

∀s ∈ SO , t ∈ Τ

(24)

∀i ∈ I O

(25)

∀i ∈ I O , t ∈ Τ

(26)

T

T

∑e ≥ ∑ z it

t =1

t

∑z g =1

i +1 g

∀i ∈ I S

(27)

∀i ∈ I O , t ∈ Τ

(28)

∀(i, j ) ∈ SECO , t ∈ Τ

(29)

∀s ∈ SO , t ∈ Τ

(30)

≤ MAX_EXTmt

∀m ∈ M O , t ∈ Τ

(31)

=0

∀s ∈ SO , t ∈ Τ \ WIN ( s ) (32)

t =1 t

it

≤ ∑ eig g =1

t

t

g =1

g =1

∑ z jg ≤ ∑ eig ∑γ

i∈I ( s )

MIN_EXTmt ≤

i

eit ≤ DPt

∑ ∑x

s∈S ( m ) k∈K

skt

∑z

i∈I ( s )

it

zit ∈ {0,1}, eit ≥ 0

∀i ∈ I O , t ∈ Τ.

(33)

    Constraints (23) and (24) define the production and costs variables respectively, similar  to constraints (8) and (9) in the underground case. Constraints (25) ensure that each bench is  removed at most once. Constraints (26) establish the logical link between variables eit  and zit  so  that the extraction of a bench cannot occur in a period prior to when it is initiated.  Constraints  (27) ensure that bench is not partially extracted. Constraints (28) require the extraction to take  place from top to bottom within an expansion. Constraints (29) impose the extraction sequence  among  benches  to  comply  with  geomechanical  requirements.  Constraints  (30)  limit  the  extraction rate per expansion and  constraints (31) impose bounds  on  the total extraction per  period. Finally, constraints (32) impose an extraction time window and constraints (33) restrict  the domain of the decision variables.    4.4 Model Discussion and Strengthened Formulation  We  begin  this  section  stating  the  model's  complexity.  The  following  proposition  justifies  the  heuristic methods described in the next section.     Proposition 2: The long‐term production planning problem for underground and open pit mines  given by Equations (1)‐(33) is strongly NP‐hard.    Given Proposition 2, we turn to the LP relaxation of the model to generate near optimal  solutions. In the case of underground mines, from Lemma 1, if ore grades are decreasing, then  most  of  the zint variables  in  the  LP  solution  have  integer  values.  There  is  a  natural  rounding  heuristic that follows from this observation which is described in the next section. Underground  mines are usually designed in such a way that the lower blocks of a column have higher grades.  16 

This  provides  support  for  using  our  heuristic  which  is  confirmed  by  its  good  performance  in  practice (see section 6.1).  For  open  pit  mines,  there  is  no  equivalent  to  Lemma  1.  In  fact,  several  benches  of  overburden material must be removed before the ore deposit is reached, and even then, the  best ore grades are usually deeper in the ground. In terms of our model, this means that the  solution of the LP relaxation is highly fractional since it is optimal to reach the lower benches as  soon  as  possible.  Indeed,  most  of  the zit variables  have  small  fractional  values  very  close  to  zero.  This  makes  it  hard  to  construct  feasible  integer  solutions  based  on  the  LP  relaxation.  In  particular, branch and bound or rounding heuristics become very inefficient. To overcome this  obstacle, we strengthen the formulation for open pit mines by adding a network structure to  the production phase. This is explained next.  For each expansion, we generate a directed graph where the nodes represent a bench‐ period grid, as shown in the left hand side of Figure 5.  At a given node, there is an outgoing arc  to all the benches that can be reached in the next period, according to the conditions imposed  by constraints (28) and  (30)‐(32). An artificial staring bench and period is added to the upper  level of the grid, and similarly an artificial ending bench and period is added to the lower part  (see the right hand side of Figure 5). Hence, timing the extraction of the expansion is equivalent  to moving one unit of flow from the upper left corner to the lower right of the grid.   

Figure 5: Extraction graph of the production phase for open pit mines.    To  formulate  the  extraction  graph  for  an  expansion s ∈ SO ,  we  first  need  to  introduce  some  additional  notation.  Let t = 0 and t = T + 1   denote  the  artificial  starting  and  ending  periods  respectively.  Similarly,  let Or ( s )   and De( s )   denote  the  artificial  starting  and  ending  benches  respectively.  Let AN (i, t )   and SU (i, t )   denote  the  antecessors  and  successors  of  node (i, t ) in the extraction graph respectively. Let wijt be a binary decision variable that equals  one if the unit of flow (i.e., the extraction) goes from bench i  to bench j  in period t . In other  17 

words,  if wijt = 1 ,  then  it  implies  that  the  extraction  of  benches i   and j   must  begin  in  periods t − 1   and t   respectively,  and  given  the  precedence  constraints  (29)  all  the  benches  in  between must be fully extracted in period t . Note that if there is an arc from node (Or ( s ), T )   to ( De( s), T + 1) , then the model could still decide not extract the expansion at all (see Figure  5). The mathematical formulation of the extraction graph is the following:   



wOr ( s ) j1 = 1

∀ s ∈ SO

(34)

w jDe ( s )T +1 = 1

∀ s ∈ SO

(35)

∀i ∈ I O , t ∈ Τ

(36)

∀i ∈ I O , t ∈ Τ

(37)

∀i, j ∈ I O , t ∈ Τ.

(38)

j∈SU ( Or ( s ),0)



j∈AN ( De ( s ),T +1)



j∈AN ( i ,t )

w jit =



j∈SU ( i ,t )

i −1

zit = ∑

wij (t +1)



h = 0 j∈SU ( h ,t −1) j ≥i

wijt ∈{0,1}

whjt

  Constraints (34) and (35) impose that a unit of flow starting from the upper left node of  the extraction graph in period t = 0 must finish at the lower right in period t = T + 1 . Constraints  (36)  establish  the  flow  conservation  at  each  node  of  the  graph.  Constraints  (37)  link  the  flow  variables wijt with  the  decision  variable zit from  the  original  formulation,  and  constraints  (38)  impose the binary condition.  By construction, constraints (34)‐(38) are valid for the original formulation of the open  pit production phase. It is well known that in the absence of other constraints, the LP relaxation  of the network formulation given by constraints (34)‐(38) yields an integer solution. When the  precedence  constraints  and  the  downstream  processes  are  considered,  some  of  the wijt variables  are  fractional.  However,  the  LP  relaxation  is  tighter  than  the  original  formulation, and more importantly, the fractional values are closer to one. The latter provides a  good starting point for the rounding heuristic that we describe in the next section.    5. SOLUTION APPROACH  Our approach to solve the model described in section 4 consists in solving the LP relaxation and  then using a heuristic rounding procedure in order to find a feasible integer solution. In sections  5.1 and 5.2 we provide the details for the underground and open pit cases respectively. Note  that in the relaxed problem we can assume with no loss of optimality that the constraints (11)  hold  as  an equality  so  that  the zint and eint variables  take  the  same  values.  Therefore,  we  omit 

18 

the  latter  (idem  for zit   and eit ).  The  procedures  described  here  below  could  apply  to  other  minerals besides copper.    5.1 Underground Mining Rounding Heuristic  In the case of underground mines, the LP relaxation does not present many fractional values.  This can be explained primarily by two reasons: (i) the model formulation itself is very tight; and  (ii) the copper grade in a column is typically decreasing, then larger benefits are obtained when  the inferior blocks are extracted first (see Lemma 1). To obtain an integer solution, the rounding  heuristic  proceeds  iteratively  from  the  inferior  blocks  to  the  higher  ones  fixing  those  that  already have an integer value. During this process the tons of ore extracted from each column is  fixed to the value obtained in the relaxed solution. Figure 6 shows the passages of the heuristic  for a particular column.   

  Figure 6: The rounding heuristic for underground mining.    5.2 Open Pit Mining Buffer‐based Heuristic  The heuristic to solve the open pit problem is also performed in an iterative fashion. As before,  the starting point is the solution of the LP relaxation. Then, a subset of extraction variables is  fixed,  and  the  LP  relaxation  is  solved  again.  This  procedure  continues  until  there  are  no  fractional variables remaining. The overall structure of the heuristic is the following:    Step 0: Initialization. Solve the LP relaxation of the open pit problem and let  zitLP ( q )  denote the  solution for the extraction variables where  q  represents a counter. Initialize  q = 1  and  initialize  Λ (m)  to contain all the expansions in the open pit mine m ∈ M O .  Step  1:  Extraction  simulation.  For  each  expansion s ∈ ∪ m∈M O Λ (m) ,  apply  Algorithm  1  to  simulate  its  extraction.  Start  with  the  inner  expansions  of  the  open  pit  and  continue  outwards.  Algorithm  1  is  described  below  in  pseudo‐code.    Note  that  IND( A)   is  an  indicator function that equals one if the argument  A  is true and is zero otherwise. The  sets  Buf1 ( s )  and  Buf 2 ( s )  are used to store the last bench that is reached in each period  of the simulated extraction of expansion s .   19 

let t = 1 let i = Or(s) let j = Or(s)+1 let Buf1(s) = Buf2(s) = φ while (t≤T and jt {z LP } jg jg

if (f1 = f2 = f3 =1) (q) let z SIM =  1 jt update downstream flows let j = j+1 else (q) let z SIM =  0 jt let Buf1 ( s ) = Buf1 ( s ) ∪ { j} and Buf 2 ( s ) = Buf 2 ( s ) ∪ {( j , t )} let i = j - 1 let t = t+1 end end Algorithm 1: Extraction simulation for an expansion s ∈ ∪ m∈M O Λ (m) .    Step  2:  Expansion  selection.  For  each  mine m ∈ M O ,  select  the  expansion  that  minimizes  the  squared difference between the extraction simulation (from the previous step) and the  solution of the  q ‐th LP relaxation. In formal terms, let  sq∗ (m) = arg min s∈Λ ( m )



i∈ I ( s ),  t∈Τ

( zitSIM ( q ) − zitLP ( q ) )2 . 

Step 3: Tighten LP formulation. Add the following two constraints to the open pit problem: 

zit = zitSIM ( q ) zit + zit +1 = 1  

∀i ∈ I ( sq* (m)) \ Buf1 ( sq* (m)), m ∈ M O , t ∈ Τ (39)

∀(i, t ) ∈ Buf 2 ( sq* (m)), m ∈ M O LP ( q +1) it

  Solve  the  q + 1 ‐th  LP  relaxation  and  let  z

(40)

  denote  the  solution  for  the  extraction 

variables.  For  each m ∈ M O ,  redefine Λ(m) = Λ(m) \{sq* (m)} .  If  ∪ m∈M O Λ (m)   is  empty,  then STOP. Otherwise, increment the counter  q = q + 1  and go to Step 1.      The premise of the heuristic, as in most rounding procedures, is that the fractional solution  of  the  LP  relaxation  is  a  good  starting  point  to  build  near‐optimal  integer  solutions.  After  initializing  the  variables  (Step  0),  the  outcome  of  the  q ‐th  LP  relaxation  is  used  in  Step  1  to  build  an  integer  solution  for  each  expansion,  which  is  denoted zitSIM ( q ) .   The  superscript  20 

SIM stands  for  "simulation"  of  the  extraction  phase.  This  is  performed  in  Algorithm  1,  which  traverses the extraction network from the upper left corner down to the lower right exit node  (see Figure 5). At a given node ( j , t ) , the algorithm must decide whether to move down (i.e., 

extract bench j + 1 ) or to the right (i.e., transition to period t + 1 ). The latter occurs by default  unless three conditions are met, which are given by the flags f1 , f 2 , and  f3  respectively.    

Flag  f1  checks whether node  j  is accessible from the node reached in the previous period, 

i.e.,  node (i, t − 1) .  Flag f 2 checks  the  precedence  constraints  (29)  and  the  availability  of  downstream capacity to process the material from bench j in period t . The last flag  f 3  checks  whether in the  q ‐th LP solution the amount extracted from bench  j  up to period  t  is greater  than the maximum of the fractional values in the remaining periods. The latter measures the  desirability  of  extracting  bench  j   in  period t .  This  rule  contrasts  with  the  usual  rounding  criterion that simply picks the largest fractional value and rounds it up to one. Figure 7 shows  four  examples  that  illustrate  the  two  criterions  for  a  given  bench j .  Each  example  shows  a  different solution of the LP relaxation. In examples (a), (b) and (c), the rounding criterion based  on  f 3   yields  the  same  outcome  as  the  criterion  that  rounds  up  the  largest  fractional  value.  However, in example (d),  f3  would allow the extraction to start in the second period, whereas  the  criterion  that  considers  the  largest  fractional  value  would  delay  the  extraction  until  the  fourth period just as in example (b). In other words, the criterion based on  f3  attempts to start  extracting  the  bench  earlier  and  will  wait  only  if  the  fractional  values  of  the  LP  solution  are  considerably skewed.  f3 holds (a)

0.7

0.1

f3 holds 0.1

0.1

bench j

0.1

periods (c)

0.4

0.2

0.2

0.1

0.1

0.7

(b)

0.4

(d)

periods 0.2

bench j

f3 holds

0.2

0.2

0.2

f3 holds

  Figure 7: Examples of the rounding criterion for open pit mining.     

In Step 2 of the heuristic, the expansions to be fixed are selected based on how much the 

extraction  simulation  zitSIM ( q )   differs  from  the  LP  solution zitLP ( q ) .  Finally,  in  Step  3,  additional  constraints are added to the open pit problem to fix the extraction variables for the expansions  selected  in  Step  2.  Constraint  (39)  sets  the  variable  zit equal to  0‐1  according  to zitSIM ( q ) .  For  a  bench that is extracted in two consecutive periods, constraint (40) is imposed, which allows the  LP relaxation to decide what fraction of the bench to extract in each period. These benches act  21 

as  buffers  between  periods  and  give  the  LP  relaxation  more  flexibility  to  satisfy  all  the  constraints. Note that when there are few expansions left to be fixed, usually it is possible to  solve the MIP formulation to optimality in reasonable time. Therefore, Algorithm 1 can be seen  as a way of reducing the MIP until it can be solved by a commercial solver.    6. IMPLEMENTATION AND IMPACT AT CODELCO  All  the  data  need  by  the  model  was  extracted  directly  from  Codelco's  databases.  The  equipment,  machinery  and  facilities  with  a  lifespan  shorter  than  the  planning  horizon  were  prorated  and  included  as  variable  costs.  Some  standard  pre‐processing  routines  similar  to  Kuchta  et  al.  (2003)  were  performed  in  order  to  reduce  the  number  of  binary  variables.  Representative instance sizes after pre‐processing are shown in Table 1.    

Model  N° Constraints  N° Variables  0‐1 Variables  Underground mine  446,521  535,639  196,386  Open pit mine  245,391  898,742  160,386  Table 1: Example of instance sizes after pre‐processing.      The model was implemented using the modeling language GAMS and was solved using the  parallel  barrier  interior  point  algorithm  of  CPLEX  on  a  computer  running  Microsoft  Windows.  Given the large timeframe spanned by the project, it has benefited from several upgrades in the  commercial software and the hardware performance. In sections 6.1 and 6.2 we describe large‐ scale case studies performed at Codelco El Teniente (underground) and Codelco North Division  (open  pit)  respectively.  In  section  6.3  we  summarize  the  current  uses  of  the  model  and  its  impact on Codelco's operations.    6.1 Case Study at El Teniente4  El Teniente is the largest underground mine in the world. It is located in the Andes Mountains  at  2,200  meters  above  the  sea  level  and  80  kilometers  south  of  Santiago,  Chile's  capital.  Its  origins date back to pre‐colonial times but officially its exploitation started in 1904. It currently  has nearly 2,400 kilometers of underground tunnels. Our case study considered all the sectors  that were in operation during the 2002 planning cycle. The horizon was divided in six periods,  the first two of 1 and 4 years respectively, and the last four of 5 years each, spanning 25 years  overall. The annual discount rate was 10% and all sunk costs were excluded. A fixed price of 85¢  per pound of fine copper was used throughout the horizon.5    Table 2 shows the results of the legacy planning approach explained in section 3.3 and the  solution  based  on  the  optimization  model.  The  comparison  considered  the  same  investment  4 5

 A preliminary version of this case study without the financial impact was presented in Epstein et al. (2003).   As of February 2010 the price was $3.4/lb. 

22 

projects in both cases. Hence, the differences come from the selection of reserves. The mining  plan obtained with the model‐based approach increased El Teniente's NPV by 5.1%.      

Model‐based vs. Legacy  Percentage Difference  Revenues  2.4%  Production Costs  1.2%  Investments  0.0%  NPV  5.1%  Table 2: Legacy vs. model‐based planning approach at El Teniente.       A  detailed  comparison  of  the  solutions  showed  that  both  provided  production  schedules  that  were  very  similar  in  terms  of  total  tons  extracted  per  period  and  sector,  but  the  model  obtained higher grades of copper and molybdenum, which translated into a higher income, as  seen in Table 2. This improvement came mostly from considering the mine as a whole and from  timing the production of each sector. The legacy approach that looked at each sector separately  did not maximize the throughput of the downstream processes.     The model‐based solution consistently reached the capacity of the intermediate processes  and  it  provided  an  appropriate  ore  blend  from  all  sectors  so  that  the  concentration  plant  received a constant supply. Moreover, the model was solved with an explicit constraint to limit  the total amount of arsenic produced (as a byproduct) each period. The legacy approach had  difficulties incorporating such kind of environmental restrictions.     6.2 Case Study at the North Division  Codelco  North  Division  is  located  1,650  kilometers  north  of  Santiago  and  at  approximately  3,000 meters above the sea level. It comprehends several open pit mines at different stages of  their  life‐cycle.  Currently,  the  two  most  important  mines  are  Chuquicamata  (better  known  as  "Chuqui") and Radomiro Tomic (RT). Chuquicamata has been in operation for nearly a century  and  its  pit  today  covers  more  than  1,300  hectares  and  is  almost  one  kilometer  deep.  In  contrast, the exploitation of RT only started in 1995, but it ramped up production very quickly.  While Chuquicamata is mainly composed of sulfides, RT's reserves are mostly oxides.     The case study analyzed at Codelco North Division considered both mines, Chuquicamata  and RT, and we used the geological, technical and economic information available in the 2004  planning  period.  We  considered  a  10‐year  planning  horizon  divided  in  yearly  periods  with  a  fixed price for copper of 85¢/lb.  As before, sunk costs were excluded from the evaluation, so  investments were mostly the acquisition of major extraction machinery including trucks, drills  and loaders.  

23 

  Like in the case study at El Teniente, we conducted an analysis to quantify the impact of  the  implementation  of  our  proposed  methodology.  However,  a  simple  comparison  of  the  net  present values is not enough to assess the economic impact. Since one of the advantages of our  approach  is  that  it  can  be  used  to  simultaneously  plan  several  mines  and  plants,  the  implementation  of  our  methodology  facilitated  the  first  joint  plan  for  Chuquicamata  and  RT.  Therefore,  to  isolate  the  impact  of  our  planning  methodology  from  the  benefits  due  to  integrating both mines we compared our proposed solution against three benchmarks:    i. Legacy‐based independent plans ( l − ind ): This is the optimized baseline where we impose  the flow levels determined by the legacy planning approach and then we solve the model  to  find  the  best  solution  under  those  conditions.    Here,  different  mines  are  planned  independently so the downstream processes are not shared.  ii. Model‐based independent plans ( m − ind ): In this instance the two mines are not allowed  to  interact  (as  in l − ind )  but  we  no  longer  impose  the  flow  levels  from  the  legacy  iii.

approach, so the production plan for each mine is dictated by the model.  Legacy‐based  integrated  plan  ( l − int ):  In  this  instance  the  production  plan  is  determined  by the legacy planning process (as in l − ind ) but the material is allowed to flow between  the two mines, so it provides an estimation of the value of integrated planning under the  legacy approach.  

   

We denoted our proposed solution m − int , which stands for model‐based integrated plan. 

In  Table  3  we  compare m − int   against  the  three  benchmarks.  The  respective  percentage  increases are shown in Figure 8. First, we note that the total economic benefit from using the  model‐based  methodology  to  obtain  an  integrated  plan  corresponds  to  a  gain  of  8.2%  with  respect  to  the  optimized  baseline l − ind .  From  Figure  8,  the  total  impact  has  the  following  breakdown:  3.2  percentage  points  (pp)  are  due  to  integrated  planning  (see l − int vs. l − ind ),  4.7pp are due to the model‐based approach (see m − ind vs. l − ind ), 0.2pp can be attributed to  synergies  between  integrated  planning  and  the  model‐based  approach  (to  see  this,  compare  the  parallel  arrows),  and  finally  0.1pp  is  the  multiplicative  effect.  In  other  words,  the  specific  impact of the model‐based approach is in the order of 5% which is consistent with the result  obtained  in  the  case  study  done  at  El  Teniente  (see  section  6.1).  The  fact  that  there  are  synergies is another interesting observation meaning that integrated planning and the model‐ based approach reinforce each other.   

24 

  m − int   l − ind   m − ind   l − int   Revenues  1.4%  0.0%  1.7%  1.1%  Production Costs  ‐1.6%  0.0%  0.8%  0.0%  Investments  ‐27.3%  0.0%  ‐31.8%  0.0%  NPV   8.2%  0.0%  4.7%  3.2%  Table 3: Economic evaluation of the model‐based plan ( m − int ) vs. benchmarks at Codelco  North Division. Figures are shown as percentage changes with respect to the baseline ( l − ind ).    Integrated Planning Yes

No

Production Plan

 

 

Legacy Approach

l − ind 4.7%

Model-Based Approach

m − ind

3.2% 8.2%

3.4%

l − int 4.9%

m − int

  Figure 8: Impact of the model‐based integrated production plan ( m − int ) on NPV with  respect to the three benchmarks: l − ind , m − ind , and l − int .    In Table 4 we compare our solution m − int  with the three benchmark plans. We observe 

that  the  total  mined  material  is  almost  identical  in  all  instances,  so  the  total  costs  are  comparable (see Table 3). This is intuitive because the ore with higher grades is located at the  bottom  of  the  pit  and  to  get  there  all  the  material  above  it  needs  to  be  extracted.  The  total  extraction might be different if the model decides not to extract a whole expansion, but that  did  not  happen  in  this  exercise.    The  economic  improvement  of  the  model‐based  over  the  legacy approach comes mainly from two sources. First, the model‐based approach requires less  investment  in  machinery  because  it  can  achieve  better  utilizations  by  smoothing  production  whereas  the  legacy  approach  needs  to  make  additional  investments  to  accommodate  higher  production peaks.6  Second and more importantly, the model‐based approach produces more  copper than the legacy approach. In particular, with respect to the baseline l − ind , the solution 

m − int   sends  5.5%  more  ore  to  the  plants  reducing  the  average  grade  by  2.6%.  Unlike  traditional  methodologies  that  use  a  cutoff  grade  per  mine  and  period,  the  model‐based  approach  simultaneously  considers  multiple  elements  to  decide  when  and  where  to  process  each  bench  of  the  pit  resulting  in  a  better  usage  of  the  available  capacity.  For  instance,  we 

6

 The reduction in investments can be seen in Table 3. Percentagewise, these are important reductions. However,  in absolute terms they are small compared to total revenue and costs. 

25 

verify that in our solutions, when the benches are close to the plant it decides to send ore of  relatively lower grade than those generated from benches that are away from the plants.    

  m − int l − ind m − ind   l − int Total Extraction (Kton)  2,882,457  2,889,825  2,874,050  2,889,825  Ore to Process (Kton)  1,468,104  1,391,384  1,456,711  1,418,244  Average Ore Grade (%)  0.660  0.678  0.663  0.674  Average Recovery Rate (%)  80.4  81.0  82.1  80.5  Copper produced (Kton)  7,792  7,643  7,926  7,693  Table 4: Mining plan summary of the model‐based vs. benchmarks at North Division. Material  flows are shown in kilo‐ton, grades shown in percentage.       An interesting observation comes from comparing the independent and integrated model‐ based  solutions  ( m − ind and m − int respectively).  The  former  has  a  better  recovery  rate  and  produces more copper but has a lower NPV than the latter. The reason is that the integrated  solution  redirects  mineral  from  Chuquicamata  to  a  plant  next  to  RT  which  had  some  slack  capacity ( 0 and ein* < 1. (41) Due to the decreasing grades, the tonnage from block n  that reaches the concentration plant is  more profitable than an equivalent amount from block n + 1 . Except for this difference at the  commercial level, the blocks are identical. Hence, shifting production from block n + 1  to block  n  is feasible since the downstream processes and the other constraints remain unaffected. In 

other words, if (41) holds, it is possible to shift a certain amount ε > 0  from  ein* +1 to ein*  and the  objective value would improve. This means that {ein* }n∈N ( i )  cannot be optimal, but that would be  a contradiction. So condition (41) cannot hold, and ein* +1  can only be positive if ein*  equals one.  Therefore,  ein* ∈ {0,1}  for all blocks except the last one extracted (i.e., except for the highest n   such  that ein* > 0 ).    From  constraints  (11)  and  (13)  we  have  that ein ≤ zin ≤ ein −1 ∀n > 0 .  This  implies  that  the  optimal  solution  of  the  LP  relaxation  must  satisfy zin ∈ {0,1}   for  all  blocks  except the first and last one extracted, which concludes the proof.  Proof of Proposition 2  Consider a single period, a single product, and a single downstream process with finite capacity.  Let  the  columns  in  the  underground  mines  have  a  single  block.  Set TCvkk ' = 1 , EMIN (i ) = 1 , 

Δ = T = 1 , WIN ( s) = Τ , δ s = ∞ , γ in = 0 , Ekt = ∞ , pkt = 1   and  set  all  the  cost  parameters  equal  to zero. Set the minimum incorporated area and extraction requirements equal to zero, and set  the  maximum  values  equal  to  infinity.  Then,  the  production  planning  model  reduces  to  the  precedence  constraint  knapsack  problem  in  which  the  profits  are  equal  to  the  weights.  The  latter  is  strongly  NP‐hard  (see  section  13.2  of  Kellerer  et  al.  2004).  Since  we  have  shown  a  polynomial reduction, the long‐term production planning problem is strongly NP‐hard as well.   30 

REFERENCES  Alford,  C.,  1995.  Optimization  in  underground  mine  design.  In  25th  APCOM,  Australian  Instit.  Mining and Metallurgy (AusIMM), Brisbane, Australia, 213–218.  Alford,  C.,  M.  Brazil,  D.H.  Lee.  2007.  Optimisation  in  Underground  Minining,  in  Handbook  of  Operations Research in Natural Resources, A. Weintraub, C. Romero, T. Bjorndal, R. Epstein,  Eds.,  561‐578  Barbaro,  R.,  R.  Ramani.  1986.  Generalized  multiperiod  MIP  model  for  production  scheduling  and processing facilities selection and location. Mining Engineering 38(2) 107–114.  Brazil,  M.,  D.  Lee,  M.  Van  Leuven,  J.  Rubinstein,  D.  Thomas,  N.  Wormald.  2003.  Optimising  declines in underground mines. Trans. IMM, Sect. A: Mining Tech. 112(3) A164–A170.  Busnach,  E.,  A.  Mehrez,  Z.  Sinuany‐Stern.  1985.  A  production  problem  in  phosphate  mining.  J.Oper. Res. Soc. 36(4) 285–288.  Caccetta,  L.,  S.P.  Hill.  2003.  An  Application  of  Branch  and  Cut  to  Open  Pit  Mine  Scheduling.  Journal of Global Optimization 27 349–365.  Caccetta, L. 2007. Applicationn of Optimisation Techniques in Open Pit Mining, in Handbook of  Operations Research in Natural Resources, A. Weintraub, C. Romero, T. Bjorndal, R. Epstein,  Eds., 546‐560.  Carlyle,  W.,  B.  Eaves.  2001.  Underground  planning  at  Stillwater  Mining  Company.  Interfaces  31(4) 50–60.  Chicoisne, R., D. Espinoza, M. Goycoolea, E. Moreno, E. Rubio. 2009. A new algorithm for the  open‐pit  mine  scheduling  problem.  Working  Paper,  Universidad  Adolfo  Ibañez  and  Universidad de Chile.  Dagdelen,  K.,  T.  Johnson.  1986.  Optimum  open  pit  mine  production  scheduling  by  Lagrangian  parameterization. In 19th APCOM, University Park, PA, 127–141.  Denby, B., D. Schofield, S. Bradford. 1991. Neural Network Applications in Mining Engineering.  Department of Mineral Resources Engineering Magazine, University of Nottingham, 13–23.  Dowd,  P.,  A.  Onur.  1993.  Open‐pit  optimization  part  1:  optimal  open‐pit  design.  Trans.  IMM,  Sect. A: Mining Indust. 102 A95–A104.  Economía  y  Negocios,  online  version  (in  Spanish).  January  15,  2010.  Accessed  on  January  18,  2010. http://www.economiaynegocios.cl/noticias/noticias.asp?id=71136.  Elevli, B. 1995. Open pit mine design and extraction sequencing by use of OR and AI concepts.  Internat. J. Surface Mining 9(4) 149–153.  Epstein,  R.,  F.  Caro,  J.  Catalán,  P.  Santibáñez,  A.  Weintraub,  S.  Gaete.  2003.  Optimizing  long‐ term  planning  for  underground  copper  mines.  In  Copper  2003  ‐  Cobre  2003,  5th  Internat.  Copper Conf., Vol I. Plenary Lectures, Economics and Appl. of Copper Proc., Canadian Instit.  of Mining, Metallurgy and Petroleum, Santiago, Chile, 265–279. 

31 

Erarslan, K., N. Celebi. 2001. A simulative model for optimum pit design. Canadian Institute of  Mining (CIM) Bull. 94(1055) 59–68.  Frimpong, S., E. Asa, J. Szymanski. 2002. Intelligent modeling: advances in open pit mine design  and optimization research. Internat. J. Surface Mining Reclamation 16(2) 134–143.  Fytas,  K.,  J.  Hadjigeorgiou,  J.  Collins.  1993.  “Production  scheduling  optimization  in  open  pit  mines.” Internat. J. Surface Mining 7(1) 1–9.  Gershon,  M.  1983.  Mine  Scheduling  Optimization  with  Mixed  Integer  Programming.  Mining  Engineering 35 351–354.  Hochbaum, D., A. Chen. 2000. Performance analysis and best implementations of old and new  algorithms for the open‐pit mining problem. Operations Research 48(6) 894–914.  Jalali,  S.,  M.  Ataee‐pour,  K.  Shahriar.  2006.  Pit  limits  optimization  using  a  stochastic  process.  CIM Magazine 1(6) 90.  Kellerer, H., U. Pferschy, D. Pisinger. 2004. Knapsack problems. Springer, 546 pages.  Kim, Y.C. 1978. Ultimate Pit Limit Design Methodologies Using Computer Models – The State of  the Art. Mining Engineering 30 1454 –1459.  Kim, Y., Y. Zhao. 1994. Optimum open pit production sequencing ‐ the current state of the art.  In  Preprint  ‐  Soc.  Mining,  Metallurgy  and  Exploration  Annual  Meeting  Proc.,  SME  of  the  AIME, Littleton, CO.  Kuchta,  M.,  A.  Newman,  E.  Topal.  2003.    Production  Scheduling  at  LKAB's  Kiruna  Mine  Using  Mixed Integer Programming.   Mining Engineering April 35‐40.  Kuchta, M., A. Newman, E. Topal. 2004. Implementing a Production Schedule at LKAB’s Kiruna  Mine. Interfaces 34 124–134.  Kumral, M., P. A. Dowd. 2005. A Simulated Annealing Approach to Mine Production Scheduling.  The Journal of the Operational Research Society 56(8) 922 – 930.  Lane, K., 1988. The Economic Definition of Ore: Cutoff Grade in Theory and Practice, Mining J.  Books Limited, London.  Laurich,  R.  1990.  Planning  and  design  of  surface  mines.  In  Surface  Mining,  B.  Kennedy,  ed.,  Society for Mining, Metallurgy and Exploration, Port City Press, Baltimore, MD, chap. 5.2,  465–469.  Lerchs,  H.,  I.  Grossmann.  1965.  Optimum  design  of  open‐pit  mines.  Canadian  Mining  and  Metallurgical Bull. LXVIII 17–24.  Lizotte, Y., J. Elbrond. 1985. Optimal layout of underground mining levels. CIM Bull. 78(873) 41– 48.  Newman, A., E. Rubio, R. Caro, A. Weintraub,  K. Eurek. 2009. A Review of Operations Research  in Mine Planning. Interfaces (to appear).  Picard,  J.C.  1976.  Maximal  Closure  of  a  Graph  and  Applications  to  Combinatorial  Problems.  Management Science 22 1268‐1272.  32 

Poniewierski, J., G. MacSporran, I. Sheppard. 2003. Optimisation of Cut‐Off Grade at Mount Isa  Mines  Limited’s  Enterprise  Mine.  Proceedings  of  12th  International  Symposium  on  Mine  Planning and Equipment Selection. Kalgoorlie, WA, 531–538.  Qinglin,  C.,  B.  Stillborg,  C.  Li.  1996.  “Optimization  of  underground  mining  methods  using  grey  theory and neural networks.” In 5th MPES, W. Hennies, L. Ayres da Silva and A. Chaves, eds.,  AusIMM, Sao Paulo, Brazil, 393–398  Rahal,  D.,  M.  Smith,  G.  Van  Hout,  A.  Von  Johannides.  2003.  The  use  of  mixed  integer  linear  programming for long‐term scheduling in block caving mines. In 31st APCOM, F. Camisani‐ Calzolari, ed., SAIMM, Cape Town, South Africa, 123–131.  Ramazan, S. 2007. The new fundamental tree algorithm for production scheduling of open pit  mines. Eur.J. Oper. Res. 177(2) 1153–1166.  Samanta,  B.,  A.  Bhattacherjee,  R.  Ganguli.  2005.  A  genetic  algorithms  approach  for  grade  control  planning  in  a  bauxite  deposit.  In  32nd  APCOM,  R.  Ganguli,  S.  Dessureault,  V.  Kecojevic and J. Dwyer, eds., Taylor & Francis Gp., Tucson, AZ, 337–342.  Sarin, S., J. West‐Hansen. 2005. The long‐term mine production scheduling problem. IIE Trans.  37(2) 109–121.  Sevim, H., D. Lei. 1998. The Problem of Production Planning in Open Pit Mines. INFOR 36 1–12.   Smith, M., I. Sheppard, G. Karunatillake. 2003. Using MIP for strategic life‐of‐mine planning of  the lead/zinc stream at Mount Isa Mines. In 31st APCOM, F. Camisani‐Calzolari, ed., SAIMM,  Cape Town, South Africa, 465–474.  Tolwinski, B., R. Underwood. 1996. A Scheduling Algorithm for Open Pit Mines. IMA Journal of  Mathematics Applied in Business and Industry 7 247‐270.  Trout, L. 1995. Underground mine production scheduling using mixed integer programming. In  25th APCOM, AusIMM, Brisbane, Australia, 395–400.  Underwood,  R.,  B.  Tolwinski.  1998.  A  mathematical  programming  viewpoint  for  solving  the  ultímate pit problem. Eur. J. Oper. Res. 107(1) 96–107.  Wang, Q., H. Sun. 2001. A theorem on open‐pit planning optimization and its application. In 29th   APCOM, H. Xie, Y. Wang and Y. Jiang, eds., Swets & Zeitlinger, Beijing, China, 295–298.  Wright, E. 1989. Dynamic programming in open pit mining sequence planning: a case study. In  21st APCOM, A. Weiss, ed., SME of the AIME, Littleton, CO, 415–422.  Yun,  Q.,  J.  Liu,  Y.  Chen,  G.  Huang.  1990.  Optimization  of  planning  and  design  in  underground  mines. In 22nd APCOM, Technische Universität Berlin, Berlin, Germany, 255. 

33