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 extractionsequence 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).

Figure 1: An extraction proﬁle in an open pit mine

Dynamics

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

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

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 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 = (x1,…,xC) is said to be admissible if the geotechnical slope
constraints are respected in the sense that

(1)

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) = x1(t),…,xC(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

(2)

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

(3)

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.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 Rentxu(t)(t),u(t). With discount rate rf≥ 0 and discount factor
0 ≤≤ 1, the optimization problem is sup u(⋅)∑t=0+∞()tRentxu(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
consider3

(4)

Dynamic programming equation

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

(5)

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 (𝕏♯ = 104), while this raises to 1020000 if we assume that
the surface consists of 100 × 100 columns (C = 10000) 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 = (x1,…,xC)φ(x) = (x1,x2− x1,…,xC− xC−1).
Let x ∈ 𝔸 and y = φ(x). Since x satisﬁes the admissibility condition (1), y necessarily satisﬁes
y1∈{1, 2} and sup c=2,…,C≤ 1. Hence, card(𝔸) ≤ 2 × 3C−1 and 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. Table 2 provides some ﬁgures.

C

D

𝕏♯

𝔸♯

𝔸♯∕𝕏♯

4

9

104

≤ 270

≤ 3%

10 000

99

1020000

≤ 2 × 1010000

2×≤ 10−10000

Table 1: Cardinals of states and acceptable states sets

A new incremental state

The dynamic programming algorithm will be released with the new state

(6)

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

3.3 Optimal trajectories simulation

State labelling

Copy the following code in a ﬁle "mine_DP_label.sce"

// -----------------------------------------// 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:NNdo 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:NNdo 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

Copy the following code in a ﬁle "mine_DP_algorithm.sce"

// -----------------------------------------// 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 matrixfor t=(TT-1):(-1):1do 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:NNdo 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 optionend