Fermer X


1 A two-dimensional open pit mine optimal extraction model

Open pit mines are generally designed with the help of so-called block models, or resource models. These latter represent the material inside the pit using millions of blocks with given shape. Open pit mines extraction is submitted to physical restrictions: only blocks at the surface may be extracted; a block cannot be extracted if the slope made with one of its neighbors is too high, due to geotechnical constraints on mine wall slopes.

In essence, exploiting an open pit mine can be seen as selecting an admissible extraction sequence of blocks where, by admissible, we mean that the above physical restrictions are satisfied. Among all admissible extraction sequences, of particular economic interest are those which are the more profitable.

Block model

For simplicity, we shall consider a two-dimensional mine [1] as in Fig. 1. Each block is a two-dimensional rectangle identified by its vertical position d ∈{1,,D} (d for depth) and by its horizontal position c ∈{1,,C} (c for column).


Figure 1: An extraction profile in an open pit mine


We assume that blocks are extracted sequentially under the following hypothesis [21]:

  • it takes one time unit to extract one block;
  • only blocks at the surface may be extracted;
  • a block cannot be extracted if the slope made with one of its two neighbors is too high, due to geotechnical constraints on mine wall slopes (the rule is mathematically decribed in the sequel);
  • a retirement option is available where no block is extracted.

States and admissible states

Denote discrete time by t = 0, 1,,T, where an upper bound for the number of extraction steps is the number C × D of blocks (it is in fact strictly lower due to slope constraints). At time t, the state of the mine is a profile

x (t) = (x1(t),...,xC(t)) ∈ 𝕏 = {1, ...,D + 1}C

where xc(t) ∈{1,,D + 1} is the vertical position of the top block with horizontal position c ∈{1,,C} (see Fig. 1). The initial profile is x(0) = (1, 1,, 1) while the mine is totally exhausted in state x = (D + 1,D + 1,,D + 1).

An admissible profile is one that respects local angular constraints at each point, due to physical requirements. A state x = (x1,,xC) is said to be admissible if the geotechnical slope constraints are respected in the sense that

|  x1 = 1    or    2    (border  slope )
   |xc+1 − xc| ≤ 1,   for   c = 1,...,C −  1
   xC =  1    or   2   (border slope).

Denote by 𝔸 ⊂{1,,D + 1}C the set of admissible states satisfying the above slope constraints (1).

Controlled dynamics

A decision is the selection of a column in = {1,,C}, the top block of which will be extracted. A decision may also be the retirement option, that we shall identify with an additional fictituous column denoted . Thus, a decision u is an element of the set = ∪{∞} = {1,,C,∞}.

At time t, if a column u(t) ∈{1,,C} is chosen at the surface of the open pit mine, the corresponding block is extracted and the profile x(t) = (x1(t),,xC(t)) becomes

x (t + 1) =   xc(t) + 1    if  c = u(t)
 c            xc(t)        else.

In case of retirement option u(t) = , then x(t + 1) = x(t) and the profile does not change. In other words, the dynamics is given by x(t + 1) = DYN(x,u) where

DYN   (x,u) =    xc + 1     if     c = u ∈ {1,...,C }
     c           xc         if     c ⁄= u    or   c = ∞.

Indeed, the top block of column c is no longer at depth xc(t) but at xc(t) + 1, while all other top blocks remain.

Starting from state x = (2, 3, 1, 4, 3) in Figure 2 and applying control u = 3, one obtains the following state (2, 3, 2, 4, 3) as in Figure 3.


Figure 2: State (2,3,1,4,3)


Figure 3: State (2,3,2,4,3)

Decision constraints

Not all decisions u(t) = c are possible either because there are no blocks left in column c (xc = D + 1) or because of slope constraints.

When in state x 𝔸, the decision u is admissible if the future profile DYN(x,u) 𝔸, namely if it satisfies the geotechnical slope constraints. This may easily be transformed into a condition u 𝔹(x), where

𝔹 (x ) := {u ∈ ℂ | DYN (x, u) ∈ 𝔸}.

Intertemporal profit maximization

The open pit mine optimal scheduling problem consists in finding a path of admissible blocks which maximizes an intertemporal discounted extraction profit. It is assumed that the value of blocks differs in depth and column because richness of the mine is not uniform among the zones as well as costs of extraction. The profit model states that each block has an economic value Rent(d,c) , supposed to be known.1 By convention Rent(d,) = 0 when the retirement option is selected. Selecting a column u(t) at the surface of the open pit mine, and extracting the corresponding block2 at depth xu(t)(t) yields the value Rent(xu(t)(t),u(t)). With discount rate rf 0 and discount factor 0 -1--
1+rf 1, the optimization problem is sup u() t=0+(--1-
1+rf)tRent(x u(t)(t),u(t)). Notice that the sum is in fact finite, bounded above by the number T 1 of blocks. Thus, we shall rather consider3

    sup        (---1-- )tRent(xu (t)(t),u(t)).
u(0),...,u(T−1)t=0 1 + rf

Dynamic programming equation

The maximization problem (4) may theoretically be solved by dynamic programming. The value function 𝒱(t,x) solves 𝒱(T,x) = 0 and

𝒱 (t,x) = max  ((---1--)tRent(x  ,u) + 𝒱 (t + 1,DYN (x,u ))).
          u∈𝔹(x)  1 + rf         u

However, due to numerical considerations (curse of dimensionality), we shall have to introduce a new, and more parcimonious, state before applying the dynamic programming algorithm.

2 Dynamic programming algorithm

Numerical considerations

To give a flavor of the complexity of the problem, notice that 4 columns (C = 4) each consisting of nine blocks (D = 9) give 10 000 states (𝕏 = 104), while this raises to 1020 000 if we assume that the surface consists of 100 × 100 columns (C = 10 000) with ninety-nine blocks (D = 99) .

However, the set 𝔸 of acceptable states has a cardinal 𝔸 which is generally much smaller than 𝕏. To see this, let us introduce the mapping x = (x 1,,xC)↦→φ(x) = (x1,x2 x1,,xC xC1). Let x 𝔸 and y = φ(x). Since x satisfies the admissibility condition (1), y necessarily satisfies y1 ∈{1, 2} and sup c=2,,C|yc|1. Hence, card(𝔸) 2 × 3C1 and the dynamic programming algorithm will be released with the new state y = (y1,,yC) ∈{1, 2}×{−1, 0, 1}C1 corresponding to the increments of the state x. Table 2 provides some figures.


4 9 104 270 3 %

10 000991020 0002 × 1010 0002×≤ 1010 000

Table 1: Cardinals of states and acceptable states sets

A new incremental state

The dynamic programming algorithm will be released with the new state

y = (y1,...,yC ) ∈ 𝕐 :=  {1,2} × {− 1,0,1}C −1

corresponding to the increments of the state x given by the inverse mapping

                         −1                                            C
y = (y1,...,yC ) ∈ 𝕐 ↦→ φ   (y) = (y1,y1 + y2,...,y1 + y2 + ⋅⋅⋅ + yC ) ∈ ℕ .

By construction, x = φ1(y) always satisfies the slope constraints (1). However, it does not always correspond to a mine profile state. Indeed, some y 𝕐 give x = φ1(y) ⁄∈ 𝕏 = {1,,D + 1}C: for instance y = (1,1,1,,1) yields x = (1, 0,1,2,,C 2).

The dynamics (x,u)↦→DYN(x,u) now becomes (y,u)↦→G(y,u) = φ(DYN(φ1(y),u)):

           {  yc + 1     if        c = u
G (y,u )c =    yc − 1     if        c = u + 1
           (  y          else.

Labelling the states

The new state set 𝕐 is a product space: we shall label its elements with integers. Let Q ∈{−1, 0, 1, 2,} be the unique integer such that

 Q+1             Q+2
3    < D  + 1 ≤ 3    .

Consider the following labelling mapping4

y = (y1,...,yC) ∈ 𝕐 ↦→  λ(y ) := (y1 − 1) +   (yi + 1)3Q+i ∈ ℕ.

The inverse mapping of λ is given as follows. Any integer z ∈{0, 1,, 3Q+C+1} may be decomposed in the basis {1, 3Q+2,, 3Q+C} as

z = z1 +     zi3Q+i    with    z1 ∈ {0,1,...,3Q+2  − 1}    and    zi ∈ {0,1,2}

for i = 2,,C. Thus (y1,,yC) = λ1(z) is given by y 1 = z1 + 1 ∈{1,, 3Q+2} and yi = zi 1 ∈{−1, 0, 1} for i = 2,,C. Denote by a ⊂{0, 1,, 3Q+C+1} those z to which corresponds an admissible state x = φ1(λ1(z)) 𝔸.

The dynamics (8) is translated in

                                     {  z + 1 − 3Q+2           if    u =  1
(z,u) ∈ ℕ ×  {1,...,C } ↦→ H (z,u ) :=   z + 3Q+u − 3Q+u+1      if    u ∈ {2, ...,C − 1}
                                     (  z + 3Q+C               if    u =  C

because of the expression (11) for z.

Let z a. Of course, z= H(z,u) corresponds to an admissible state (z′∈ a) if and only if u is admissible for the state x = φ1(λ1(z)) 𝔸, property that we shall denote by u 𝕌a(z) ⊂{1,,C}.

3 Simulating extraction strategies

3.1 Data

Copy the file "marvin_7440_valores.txt" in your local directory. Then, create a file "data.sce" in your local directory containing the following code.

  // -----------------------------------------
  // 2D agregated Marvin mine
  for i=1:5 do
    for j=1:11 do
  // -----------------------------------------
  // -----------------------------------------
  TT=HH*(HH+2)+1;// horizon
  TT=HH*NN;// horizon
  // -----------------------------------------
  // -----------------------------------------
  // 1500 extractions a year
  // 6*4 extractions a year

3.2 Upper and lower bounds for the value of the mine

The so called value of the mine is 𝒱(0,x(0)) given by (4).

Question 1 Rearrange the elements of the matrix Rent(d,c) and deduce an upper bound for the value of the mine.

  sorted_richness=sort(rich);// richness sorted by decreasing order
  UB=sum((discount .^(0:(blocks-1))) .*sorted_richness);

3.3 Optimal trajectories simulation

State labelling

Copy the following code in a file "mine_DP_label.sce"

  // -----------------------------------------
  // -----------------------------------------
  MM=3^{NN};// number of states
  Integers=(1:MM)';// labelled states
  // Will contain, for each integer z in Integers,
  // a sequence s=(s_1,...,s_NN)  with
  // s_1 \in \{1,2\}
  // and s_k \in \{0,1,2\} for k \geq 2, such that
  // 1) z = s_1 + s_2*3^{1} + ... + s_NN*3^{NN-1}
  // 2) a mine profile p_1,...,p_NN is given by
  // p_k = y_1 + ... + y_k where y_1= s_1-1 \in \{0,1\}
  // and y_j = s_j - 1 \in \{-1,0,1\} for j>1.
  // Will contain, for each integer z in Integers,
  // the sequence y=(y_1,...,y_NN).
  // The initial profile is supposed to be p(0)=(0,0,...,0)
  // to which corresponds y(0)=(0,0,...,0) and
  // s(0)=(1,1,...,1) and z(0)=1+ 3^{1} + ... + 3^{NN-1}.
  // Will contain, for each integer z in Integers,
  // a lower aproximation of z in the basis
  // 1, 3^{1},..., 3^{NN-1}
  // Partial_Integer(z,k)=s_1 + s_2*3^{1} +...+ s_k*3^{k-1}.
  // Will contain, for each integer z in Integers,
  // the image by the dynamics under the control consisting in
  // extracting block in the corresponding column.
  // State(z,1)=s_1
  // Partial_Integer(z,1)=s_1
  // Dynamics (with a "maxi" because some integers in
  // Integers+1-3^{1} do not correspond to "mine profiles").
  for k=2:NN do
    // s_{k} + s_{k+1}*3 + ...
    // State(:,k)=s_{k}
    // Partial_Integer(z,k)=s_1+s_2*3^{1}+...+s_k*3^{k-1}
    // Dynamics (with a "maxi" because some integers
    // in Integers+3^{k-1}-3^{k}
    // do not correspond to "mine profiles")
  // Correction for the dynamics
  // when the last column NN is selected.
  // Dynamics (with a "mini" because some integers in
  // Integers+3^{NN-1} do not correspond to "mine profiles").
  // -----------------------------------------
  // -----------------------------------------
  function z=profile2integer(p)
    // p : profile as a row vector
    z=sum(ss .*[1,3 .^{[1:(NN-1)]}]);
  // -----------------------------------------
  // -----------------------------------------
  // Mine profiles are those for which
  // HH \geq p_1 \geq 0,..., HH \geq p_NN \geq 0
  // that is, HH \geq y_1 + ... + y_k \geq 0 for all k
  // Since, starting from the profile p(0)=(0,0,...,0)),  the
  // following extraction rule will always give "mine profiles",
  // we shall not exclude other unrealistic profiles.
  // Each line contains a profile, realistic or not.
  adm_bottom=bool2s(Profiles < HH);
  // A block at the bottom cannot be extracted:
  // an element in adm_bottom is one if and only if
  // the top block of the column is not at the bottom.
  // Given a mine profile, extracting one block at the surface
  // is admissible if the slope is not too high.
  // Extracting block in column 1 is admissible
  // if and only if p_1=0.
  // Extracting block in column j 1<j<NN) is not admissible
  // if and only if y_j=1 or y_{j+1}=-1 that is,
  // (s_j-1)=1 or (s_{j+1}-1)=-1.
  // Extracting block in column j is admissible
  // if and only if s_j < 2 and s_{j+1} > 0.
  // Extracting block in column NN is admissible
  // if and only if p_{NN}=0.
  // Corresponds to side columns 1 and NN, for which only the
  // original top block may be extracted:
  // an element in columns 1 and NN of Admissible is one
  // if and only if the pair (state,control) is admissible.
  Admissible(:,2:($-1))=bool2s(State(:,2:($-1)) < 2 & State(:,3:$) > 0);
  // An element in column 1<j<NN of AA is one if and only
  // s_j < 2 and s_{j+1} > 0.
  Admissible=Admissible .*adm_bottom;
  // An element in column j of admissible is zero
  // if and only if
  // extracting block in column j is not admissible,
  // else it is one.
  // Labels of states for which no decision is admissible,
  // hence the decision is the retirement option
  // -----------------------------------------
  // -----------------------------------------
  // Each line contains a profile, forced to be realistic.
  // This trick is harmless and useful
  // to fill in the instantaneous gain matrix.
  for uu=1:NN do
    instant_gain(:,uu)=Admissible(:,uu) .* ...
                       richness(Forced_Profiles(Future_Integer(:,uu),uu),uu)+ ...
  // When the control uu is admissible,
  // instant_gain is the richness of the top block of column uu.
  // When the control uu is not admissible, instant_gain
  // has value "bottom", approximation of -infinity.

Dynamic programming algorithm

Copy the following code in a file "mine_DP_algorithm.sce"

  // -----------------------------------------
  // -----------------------------------------
  VV=zeros(MM,TT);// value functions in a matrix
  // The final value function is zero.
  UUopt=(NN+1)*ones(MM,TT);// optimal controls in a matrix
  for t=(TT-1):(-1):1 do
    // will contain the vector to be maximized
    // The value attributed to the value function VV(:,t)
    // when a control is not admissible.
    for uu=1:NN do
      loc=[loc, ...
           Admissible(:,uu) .* ...
           (discount^t*instant_gain(:,uu)+VV(Future_Integer(:,uu),t+1)), ...
           (1-Admissible(:,uu)) .*(discount^t*bottom+loc_bottom)];
    // When the control uu is admissible,
    // loc is the usual DP expression.
    // When the control uu is not admissible,
    // loc is the  DP expression
    // with both terms at the lowest values.
    // Adding an extra control/column which provides zero
    // instantaneous gain and does not modify the state:
    // retire option.
    [lhs,rhs]=maxi(loc,"c");// DP equation
    // retire option


Copy the following code in a file "mine_visualize.sce"

  // exec mine_visualize.sce
  // -----------------------------------------
  // Visualizing the richness matrix
  // -----------------------------------------
  // Margin=0.5 ;
  rich_color=floor((range-1)*(maxi(richness)-richness) ./(maxi(richness)-mini(richness)))+1;
  // varies between 1 and range:
  // maxi(richness) corresponds to 1
  // mini(richness) corresponds to range
  // To change the color convention, do
  // rich_color=range-rich_color +1 ;
  // gives too much color
  // xset("window",10) ;xbasc(); plot(vect_richness);
  // xtitle("Distribution of colors");
  visual_range=2*range;// to avoid white
  // table of colors
  xbasc();plot2d2(-Margin+0:(NN+2),-[0,zeros(1,NN),0,0],rect = [0,-HH-Margin,NN+1,Margin])
  xtitle("Mine richness");
  function display_profile(xx,time)
    // Displays graphics of mine profiles
    // xset("colormap",hotcolormap(64));
    plot2d2(-Margin+0:(NN+2),-[0,xx(time,:),0,0],rect = [0,-HH-Margin,NN+1,Margin]);
    // Just to to set the frame
    // The value of the mine blocks in color.
    // Will be in background.
    plot2d2(-Margin+0:(NN+2),-[0,xx(time,:),0,0],rect = [0,-HH-Margin,NN+1,Margin]);

Optimal trajectories simulation

Copy the following code in a file "mine_DP_trajectories.sce"

  // exec mine_DP_trajectories.sce
  // -----------------------------------------
  // -----------------------------------------
  // initial profile
  // xset("colormap",hotcolormap(64));
  // table of colors
  xbasc();plot2d2(-Margin+0:(NN+2),-[0,xx(1,:),0,0],rect = [0,-HH-Margin,NN+1,Margin])
  while last_control < NN+1 do
    if uu(t) < NN+1 then
    //  halt()
    //  xbasc();
    //  Matplot1(rich_color,[0,-HH-Margin,NN+1,Margin]);
    //  plot2d2(-Margin + 0:(NN+2),-[0 xx(t+1,:) 0 0],...
    //   rect=[0,-HH-Margin,NN+1,Margin]);
    //  e=gce();
    //  p=e.children;
    //  p.thickness=3;
    // xtitle("Optimal extraction profile for annual interest rate "...
    //       +string(interest)+"%")
    //  xpause(400000);
  //xtitle("Optimal extraction profile for annual interest rate "...
  //       +string(interest)+"%")
  xtitle("Optimal final extraction profile at time t="+string(time-1)+ ...
         " for annual interest rate "+string(interest)+"%")
  for time=1:final_time do
    xtitle("Optimal extraction profile at time t="+string(time)+" for annual interest rate "+ ...

Question 2 Copy the following code in a file "mine_DP.sce", then execute it. Change the economic model and the discount rate.

  // exec mine_DP.sce
  // try higher values
  exec data.sce
  // -----------------------------------------
  // -----------------------------------------
  exec mine_DP_label.sce
  exec mine_DP_algorithm.sce
  // do it once, then comment the above line
  // -----------------------------------------
  // -----------------------------------------
  exec mine_visualize.sce
  exec mine_DP_trajectories.sce


[1]   M. De Lara and L. Doyen. Sustainable Management of Natural Resources. Mathematical Models and Methods. Springer-Verlag, Berlin, 2008.

[2]   G.C. Goodwin, M.M. Seron, R.H. Middleton, M. Zhang, B.F. Hennessy, P.M. Stone, and M. Menabde. Receding horizon control applied to optimal mine planning. Automatica, 42(8):1337–1342, 2006.

L'École des Ponts ParisTech est membre fondateur de

L'École des Ponts ParisTech est certifiée 

ParisTech ParisEst ParisTech


Copyright 2014
École des Ponts ParisTech
All rights reserved