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 satisﬁed. Among all admissible extraction sequences, of particular economic interest are those which are the more proﬁtable.

Block model

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

Dynamics

We assume that blocks are extracted sequentially under the following hypothesis [2, 1]:

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 proﬁle

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

An admissible proﬁle is one that respects local angular constraints at each point, due to physical requirements. A state x = (x₁,…,x_C) is said to be admissible if the geotechnical slope constraints are respected in the sense that

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 ﬁctituous 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 proﬁle x(t) =

x₁(t),…,x_C(t)

becomes

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

Indeed, the top block of column c is no longer at depth x_c(t) but at x_c(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.

Decision constraints

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

When in state x ∈ 𝔸, the decision u ∈ℂ is admissible if the future proﬁle DYN(x,u) ∈ 𝔸, namely if it satisﬁes the geotechnical slope constraints. This may easily be transformed into a condition u ∈ 𝔹(x), where

Intertemporal proﬁt maximization

The open pit mine optimal scheduling problem consists in ﬁnding a path of admissible blocks which maximizes an intertemporal discounted extraction proﬁt. It is assumed that the value of blocks diﬀers in depth and column because richness of the mine is not uniform among the zones as well as costs of extraction. The proﬁt model states that each block has an economic value Rent(d,c) ∈ ℝ, supposed to be known.¹ 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 block² at depth x_u(t)(t) yields the value Rent

x_u(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 ﬁnite, bounded above by the number T − 1 of blocks. Thus, we shall rather consider³

Dynamic programming equation

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

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 ﬂavor of the complexity of the problem, notice that 4 columns (C = 4) each consisting of nine blocks (D = 9) give 10 000 states (𝕏^♯ = 10⁴), while this raises to 10^20 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₁,…,x_C)

φ(x) = (x₁,x₂ − x₁,…,x_C − x_C−1). Let x ∈ 𝔸 and y = φ(x). Since x satisﬁes the admissibility condition (1), y necessarily satisﬁes y₁ ∈{1, 2} and sup _c=2,…,C |yc|

≤ 1. Hence, card(𝔸) ≤ 2 × 3^C−1 and the dynamic programming algorithm will be released with the new state y = (y₁,…,y_C) ∈{1, 2}×{−1, 0, 1}^C−1 corresponding to the increments of the state x. Table 2 provides some ﬁgures.

A new incremental state

By construction, x = φ⁻¹(y) always satisﬁes the slope constraints (1). However, it does not always correspond to a mine proﬁle state. Indeed, some y ∈ 𝕐 give x = φ⁻¹(y) ⁄∈ 𝕏 = {1,…,D + 1}^C: for instance y = (1,−1,−1,…,−1) yields x = (1, 0,−1,−2,…,C − 2).

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

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

for i = 2,…,C. Thus (y₁,…,y_C) = λ⁻¹(z) is given by y₁ = z₁ + 1 ∈{1,…, 3^Q+2} and y_i = z_i − 1 ∈{−1, 0, 1} for i = 2,…,C. Denote by ℤ_a ⊂{0, 1,…, 3^Q+C+1} those z to which corresponds an admissible state x = φ⁻¹

λ⁻¹(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 = φ⁻¹

λ⁻¹(z)

∈ 𝔸, property that we shall denote by u ∈ 𝕌_a(z) ⊂{1,…,C}.

3 Simulating extraction strategies

3.1 Data

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

  // -----------------------------------------

  // 2D agregated Marvin mine
  RICHNESS=zeros(5,11);
  richness=fscanfMat("marvin_7440_valores.txt");
  for i=1:5 do
    for j=1:11 do
      RICHNESS(i,j)=sum(richness([((i-1)*4+1):mini((i*4),17)],[((j-1)*6+1):mini((j*6),61)]));
    end
  end
  richness=RICHNESS;

  // -----------------------------------------
  // RICHNESS MATRIX SIZE
  // -----------------------------------------

  [HH,NN]=size(richness);
  TT=HH*(HH+2)+1;// horizon
  TT=HH*NN;// horizon


  // -----------------------------------------
  // DISCOUNT FACTOR
  // -----------------------------------------

  interest=10;
  annual_discount=1/(1+(interest/100));

  discount=(annual_discount)^(1/1500);
  // 1500 extractions a year

  discount=(annual_discount)^(1/(6*4));
  // 6*4 extractions a year

3.2 Upper and lower bounds for the value of the mine

  //----------------------------------------
  // UPPER BOUND VALUE OF THE MINE
  //----------------------------------------

  blocks=prod(size(richness));

  rich=matrix(richness,1,blocks);
  sorted_richness=sort(rich);// richness sorted by decreasing order
  sorted_richness=maxi(0,sorted_richness);

  UB=sum((discount .^(0:(blocks-1))) .*sorted_richness);

3.3 Optimal trajectories simulation

State labelling

  // -----------------------------------------
  // LABELLING STATES AND DYNAMICS
  // -----------------------------------------

  MM=3^{NN};// number of states

  Integers=(1:MM)';// labelled states

  State=zeros(MM,NN);
  // 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.

  Increments=zeros(MM,NN);
  // 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}.

  Partial_Integer=zeros(MM,NN);
  // 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}.

  Future_Integer=zeros(MM,NN);
  // 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(:,1)=pmodulo(Integers,3^{1});
  // State(z,1)=s_1
  Partial_Integer(:,1)=State(:,1);
  // Partial_Integer(z,1)=s_1
  Future_Integer(:,1)=maxi(1,Integers+1-3^{1});
  // Dynamics (with a "maxi" because some integers in
  // Integers+1-3^{1} do not correspond to "mine profiles").
  Increments(:,1)=State(:,1)-1;

  for k=2:NN do
    remainder=(Integers-Partial_Integer(:,k-1))/3^{k-1};
    // s_{k} + s_{k+1}*3 + ...
    State(:,k)=pmodulo(remainder,3);
    // State(:,k)=s_{k}
    Increments(:,k)=State(:,k)-1;
    Partial_Integer(:,k)=Partial_Integer(:,k-1)+3^{k-1}*State(:,k);
    // Partial_Integer(z,k)=s_1+s_2*3^{1}+...+s_k*3^{k-1}
    Future_Integer(:,k)=maxi(1,Integers+3^{k-1}-3^{k});
    // Dynamics (with a "maxi" because some integers
    // in Integers+3^{k-1}-3^{k}
    // do not correspond to "mine profiles")
  end

  Future_Integer(:,NN)=mini(MM,Integers+3^{NN-1});
  // 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").



  // -----------------------------------------
  // FROM PROFILES TO INTEGERS
  // -----------------------------------------

  function z=profile2integer(p)
    // p : profile as a row vector
    yy=p-[0,p(1:($-1))];
    ss=yy+1;
    z=sum(ss .*[1,3 .^{[1:(NN-1)]}]);
  endfunction


  // -----------------------------------------
  // ADMISSIBLE INTEGERS
  // -----------------------------------------

  // 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.

  Admissible=zeros(MM,NN);

  Profiles=cumsum(Increments,"c");
  // 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.

  Admissible(:,1)=bool2s(Profiles(:,1)==0);
  Admissible(:,NN)=bool2s(Profiles(:,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.


  Stop_Integers=Integers(prod(1-Admissible,"c")==1);
  // Labels of states for which no decision is admissible,
  // hence the decision is the retirement option


  // -----------------------------------------
  // INSTANTANEOUS GAIN
  // -----------------------------------------

  Forced_Profiles=mini(HH,maxi(1,Profiles));
  // Each line contains a profile, forced to be realistic.
  // This trick is harmless and useful
  // to fill in the instantaneous gain matrix.

  instant_gain=zeros(MM,NN);

  for uu=1:NN do
    instant_gain(:,uu)=Admissible(:,uu) .* ...
                       richness(Forced_Profiles(Future_Integer(:,uu),uu),uu)+ ...
                       (1-Admissible(:,uu))*bottom;
  end
  // 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

  // -----------------------------------------
  // DYNAMIC PROGRAMMING ALGORITHM
  // -----------------------------------------

  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
    loc=[];
    // will contain the vector to be maximized
    loc_bottom=mini(VV(:,t+1));
    // 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)];
    end
    // 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.
    //
    loc=[loc,VV(:,t+1)+discount^t*0];
    // 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
    VV(:,t)=lhs;
    UUopt(:,t)=rhs;
    //
    UUopt(Stop_Integers,t)=(NN+1)*ones(Stop_Integers);
    // retire option
  end

Visualization

  // exec mine_visualize.sce

  lines(0);

  [HH,NN]=size(richness);

  // -----------------------------------------
  // Visualizing the richness matrix
  // -----------------------------------------

  Margin=0.01;
  // Margin=0.5 ;

  range=6;
  range=16;
  range=64;

  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
  prop=1/2;
  prop=0.2;
  rich_color=prop*range+(1-prop)*rich_color;


  vect_richness=sort(matrix(rich_color,HH*NN));
  // xset("window",10) ;xbasc(); plot(vect_richness);
  // xtitle("Distribution of colors");

  visual_range=range;
  visual_range=2*range;// to avoid white
  xset("colormap",hotcolormap(visual_range));
  // table of colors

  xset("window",0);
  xbasc();plot2d2(-Margin+0:(NN+2),-[0,zeros(1,NN),0,0],rect = [0,-HH-Margin,NN+1,Margin])
  Matplot1(rich_color,[0,-HH-Margin,NN+1,Margin]);
  xtitle("Mine richness");


  function display_profile(xx,time)
    // Displays graphics of mine profiles
    xset("window",time);
    xset("colormap",hotcolormap(visual_range));
    // xset("colormap",hotcolormap(64));
    xbasc();
    plot2d2(-Margin+0:(NN+2),-[0,xx(time,:),0,0],rect = [0,-HH-Margin,NN+1,Margin]);
    // Just to to set the frame
    Matplot1(rich_color,[0,-HH-Margin,NN+1,Margin]);
    // 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]);
    e=gce();
    p=e.children;
    p.thickness=3;
  endfunction

Optimal trajectories simulation

  // exec mine_DP_trajectories.sce

  // -----------------------------------------
  // OPTIMAL TRAJECTORIES
  // -----------------------------------------

  xx=zeros(TT,NN);
  zz=zeros(TT,1);
  uu=(NN+1)*ones(TT,1);
  vv=0;

  xx(1,:)=zeros(1,NN);
  // initial profile

  xset("window",999);
  // xset("colormap",hotcolormap(64));
  xset("colormap",hotcolormap(visual_range));
  // table of colors
  xbasc();plot2d2(-Margin+0:(NN+2),-[0,xx(1,:),0,0],rect = [0,-HH-Margin,NN+1,Margin])
  Matplot1(rich_color,[0,-HH-Margin,NN+1,Margin]);

  t=1;last_control=0;
  //
  while last_control < NN+1 do
    zz(t)=profile2integer(xx(t,:));
    uu(t)=UUopt(zz(t),t);
    xx(t+1,:)=xx(t,:);
    if uu(t) < NN+1 then
      xx(t+1,uu(t))=xx(t,uu(t))+1;
      vv=vv+discount^t*richness(xx(t+1,uu(t)),uu(t));
    end
    //  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);
    //
    last_control=uu(t);
    t=t+1;
  end
  //xtitle("Optimal extraction profile for annual interest rate "...
  //       +string(interest)+"%")

  final_time=t;

  time=1+final_time;
  display_profile(xx,time)
  xtitle("Optimal final extraction profile at time t="+string(time-1)+ ...
         " for annual interest rate "+string(interest)+"%")

  for time=1:final_time do
    display_profile(xx,time)
    xtitle("Optimal extraction profile at time t="+string(time)+" for annual interest rate "+ ...
           string(interest)+"%")
    xs2eps(time,'optimal_'+string(time)+'.eps');
  end

Optimal scheduling for open pit mine extraction

Contents

1 A two-dimensional open pit mine optimal extraction model