Constructing Branching Trees of Geostatistical Simulations

Abstract

This paper proposes the use of multi-stage stochastic programming with recourse for optimised strategic open-pit mine planning. The key innovations are, firstly, that a branching tree of geostatistical simulations is developed to take account of uncertainty in ore grades, and secondly, scenario reduction techniques are applied to keep the trees to a manageable size. Our example shows that different mine plans would be optimal for the downside case when the deposit turns out to be of lower grade than expected compared to when it is of higher grade than expected. Our approach further provides the probabilities of these outcomes; that is, the idea is to move toward adaptive mine planning rather than just producing a single mine plan.

Appendices

Appendix 1: Expected NPVs for Each Scenario, for the Three Prices

See Table 8.

Table 8 The NPVs of the 45 scenarios for the three commodity prices used

Appendix 2: Integral Range

The variability in average grades from one realisation to another and in grade–tonnage curves in Fig. 4 is due to the small size of the orebody (1,110 m $\times$ 90 m) relative to the range of the variogram (150 m). The orebody is not big enough to reproduce the expected grade in each realisation. That is, there is a lack of ergodicity: the average grade calculated over a realisation (spatial average) does not coincide with the mathematical expectation of the grade (average over realisations).

To understand what is happening, we went back to first principles and reread the two works by Lantuéjoul (1991, 2002) that relate the concepts of ergodicity to those of integral range. Let Z(V) be the average value of a spatial variable (ore grade) over region V. The integral range, IR, is defined as

$$\begin{aligned} IR=\lim _{[V]\rightarrow \infty } [V] \frac{VarZ(V)}{\sigma ^2}\, \end{aligned}$$

(3)

with [V] denoting the volume of V and $\sigma ^2=C(0)$ the variance of the variable at a point support.

Lantuéjoul, who refers back to Yaglom (1987), explains that this limit exists for all usual covariance functions, is non-negative, has the dimension of a length in a 1D space, an area in 2D and a volume in 3D, and is an indicator of ergodicity. This stems from the fact that, if $0<IR<+\infty$, one has

$$\begin{aligned} VarZ(V)\approx \frac{\sigma ^2}{n}\, \end{aligned}$$

(4)

with $n=[V]/IR$, as [V] is very large: the variance of Z(V) is the same as that of the average of n independent grade measurements and becomes small as n becomes large, i.e., as the region V is large in comparison with the integral range.

As the integral range of the unit-range spherical model is $\pi /5$ (Table 9), the integral range in our 2D case is:

$$\begin{aligned} IR=1.0 \times 150^2\times \frac{\pi }{5} = 14{,}137.2 \end{aligned}$$

(5)

Table 9 Integral ranges

Now the volume of our 2D deposit is $1,110 \times 90 = 99,900$ $\hbox {m}^2$, that is, 7.07 ($=99,900/14,137.2$) times the integral range. So it only “allows for 7.07 independent repetitions”, which is clearly not enough to get a spatial average with little fluctuations across the realisations. Consequently we would advise readers to compute the integral range of their covariance/variogram and compare it to the volume of the deposit being considered in order to be aware of expected statistical fluctuations.

Appendix 3: Mine Planning Optimisation

In this section we show the optimisation model used in Sect. 3.3. For that, we first introduce a base model where ore grades have no uncertainty, and later we extend the model to a formulation where the ore grades follow a scenario tree.

Importantly, we take advantage of the fact that the computational experiments are performed for a simplified mine with only three levels, so the non-anticipativity constraints mentioned in Sect. 3.3 can be implicitly imposed by taking augmented indices (reflecting the indexation of the branches in the scenario tree) on the decision variables; compare e.g. equations (11) and (12) below.

Base Optimisation Problem Without Uncertainty

We have three levels, where level $l=1$ is at the top, $l=2$ is in the middle, and $l=3$ at the bottom. Each level $l=1,2,3$ consists of a set $B_l$ of blocks, and for each block $b_l$ in $B_l$ we have to decide whether or not we extract it, represented as $x_{b_l}^{\mathrm{ext}} = 1$ or $=0$, respectively, and what fraction $x_{b_l}^{\mathrm{proc}} \in [0,1]$ of it we process; see (11). We can only process blocks that have been extracted, see (7). Each block b has a grade $g_b$, whose unit selling price is p, and the extraction and processing costs of a block are $c^{\mathrm{ext}}$ and $c^{\mathrm{proc}}$, respectively, see (6). Each level l has a total extraction capacity of $K_l$, and each extracted block takes $k_l$ of it, see (8); and analogously for the processing capacity, see (9). Lastly, there is a minimal set of precedence constraints, represented as $(b^\text {below}, b^\text {above}) \in P$, where in order to extract block $b^\text {below}$ we also have to extract block $b^\text {above}$; see (10). In this way, the base optimisation model is as follows.

$$\begin{aligned} \max \quad&\sum _{l = 1}^3 \sum _{b_l \in B_l} \left[ (p \, g_{b_l} - c^{\mathrm{proc}}) x_{b_l}^{\mathrm{proc}} - c^{\mathrm{ext}} x_{b_l}^{\mathrm{ext}} \right]&\end{aligned}$$

(6)

$$\begin{aligned} \text {s.t.} \quad&x_{b_l}^{\mathrm{proc}} \le x_{b_l}^{\mathrm{ext}}&\text {level } l=1,2,3, \text { block } b_l \in B_l \end{aligned}$$

(7)

$$\begin{aligned}&\sum _{b_l \in B_l} k_l^{\mathrm{ext}} x_{b_l}^{\mathrm{ext}} \le K_l^{\mathrm{ext}}&\text {level } l = 1, 2, 3 \end{aligned}$$

(8)

$$\begin{aligned}&\sum _{b_l \in B_l} k_l^{\mathrm{proc}} x_{b_l}^{\mathrm{proc}} \le K_l^{\mathrm{proc}}&\text {level } l = 1, 2, 3 \end{aligned}$$

(9)

$$\begin{aligned}&x_{b^\text {below}}^{\mathrm{ext}} \le x_{b^\text {above}}^{\mathrm{ext}}&\text {precedence } (b^\text {below}, b^\text {above}) \in P \end{aligned}$$

(10)

$$\begin{aligned}&x_{b_l}^{\mathrm{ext}} \in \{0, 1\}, \quad x_{b_l}^{\mathrm{proc}} \in [0,1]&\text {level } l=1,2,3. \end{aligned}$$

(11)

Optimisation Problem with Scenario Tree Uncertainty

We now extend the base optimisation model (6)–(11) to the case where there is a scenario tree for the grades of each level. Specifically, we now assume that there are $N_1$ scenarios for the grades vector in level 1; then, for each of these $N_1$ scenarios, there are $N_2$ scenarios for the vector of grades in level 2 (so there are a total of $N_1 \times N_2$ scenarios for level 2 grades); and for each of these scenarios there are $N_3$ scenarios for the vector of grades in level 3 (so there are a total of $N_1 \times N_2 \times N_3$ scenarios for level 3 grades). In this way, the scenarios for the grades of level 1 are listed by index $s_1 \in S_1 = \{1, \ldots , N_1\}$, the scenarios for grades of level 2 are listed by index $(s_1, s_2) \in S_1 \times S_2$ where $S_2 = \{1, \ldots , N_2\}$, and the scenarios for grades of level 3 are listed by index $(s_1, s_2, s_3) \in S_1 \times S_2 \times S_3$ where $S_3 = \{1, \ldots , N_3\}$.

We are assuming that the extraction of level l “sees” or “knows” the ore grade of the previous levels but not the current one, but the processing of blocks in a level sees the grade in that level (as well as in previous levels). In this way, the variables and constraints in (11) are now

$$\begin{aligned} x_{b_1}^{\mathrm{ext}}, \ x_{b_2 , s_1}^{\mathrm{ext}}, \ x_{b_3 , (s_1, s_2)}^{\mathrm{ext}} \in \{0, 1\}, x_{b_1 , s_1}^{\mathrm{proc}}, \ x_{b_2 , (s_1, s_2)}^{\mathrm{proc}}, \ x_{b_3 , (s_1, s_2, s_3)}^{\mathrm{proc}} \in [0, 1], \end{aligned}$$

(12)

for $s_1 \in S_1$, $s_2 \in S_2$ and $s_3 \in S_3$. This reflects, for instance, that for level 2 the extraction decision $x_{b_2 , s_1}^{\mathrm{ext}}$ does depend on the scenario $s_1$ for level 1 (so “it sees” the ore grades of level 1), while the processing decision $x_{b_2 , (s_1, s_2)}^{\mathrm{proc}}$ does depend on the scenario $(s_1,s_2)$ of level 2 (so “it sees” the grades of level 2).

It follows that constraints (7)–(10) are replicated according to the additional indexing of the scenarios for each variable. Indeed, constraint (7) is now

$$\begin{aligned} x_{b_1, s_1}^{\mathrm{proc}} \le x_{b_1}^{\mathrm{ext}}, \quad x_{b_2, (s_1,s_2)}^{\mathrm{proc}} \le x_{b_2, s_1}^{\mathrm{ext}}, \quad x_{b_3, (s_1,s_2,s_3)}^{\mathrm{proc}} \le x_{b_3, (s_1,s_2)}^{\mathrm{ext}} \end{aligned}$$

(13)

for all $b_1 \in B_1$, $b_2 \in B_2$ and $b_3 \in B_3$, and all $s_1 \in S_1$, $s_2 \in S_2$ and $s_3 \in S_3$; constraint (8) is now

$$\begin{aligned}&\sum _{b_1 \in B_1} k_1^{\mathrm{ext}} x_{b_1}^{\mathrm{ext}} \le K_1^{\mathrm{ext}}, \quad \sum _{b_2 \in B_2} k_2^{\mathrm{ext}} x_{b_2, s_1}^{\mathrm{ext}} \le K_2^{\mathrm{ext}}, \nonumber \\&\quad \sum _{b_3 \in B_3} k_3^{\mathrm{ext}} x_{b_3, (s_1,s_2)}^{\mathrm{ext}} \le K_3^{\mathrm{ext}} \end{aligned}$$

(14)

for all $s_1 \in S_1$ and $s_2 \in S_2$; constraint (9) is now

$$\begin{aligned}&\sum _{b_1 \in B_1} k_1^{\mathrm{proc}} x_{b_1, s_1}^{\mathrm{proc}} \le K_1^{\mathrm{proc}}, \quad \sum _{b_2 \in B_2} k_2^{\mathrm{proc}} x_{b_2, (s_1,s_2)}^{\mathrm{proc}} \le K_2^{\mathrm{proc}}, \nonumber \\&\quad \quad \sum _{b_3 \in B_3} k_3^{\mathrm{proc}} x_{b_3, (s_1,s_2,s_3)}^{\mathrm{proc}} \le K_3^{\mathrm{proc}} \end{aligned}$$

(15)

for all $s_1 \in S_1$, $s_2 \in S_2$ and $s_3 \in S_3$; and constraint (10) also changes to reflect additional scenario indexing of the variables $x_{b^\text {below}}^v$ and $x_{b^\text {above}}^{\mathrm{ext}}$ for each precedence pair $(b^\text {below}, b^\text {above}) \in P$. Finally, the objective function (6) now takes the average along each level of the scenario tree, so it is now as follows.

$$\begin{aligned}&\sum _{s_1 \in S_1} \frac{1}{N_1} \sum _{b_1 \in B_1} \left[ (p \, g_{b_1, s_1} - c^{\mathrm{proc}}) x_{b_1, s_1}^{\mathrm{proc}} - c^{\mathrm{ext}} x_{b_1}^{\mathrm{ext}} \right] \nonumber \\&\quad + \sum _{s_1 \in S_1, s_2 \in S_2} \frac{1}{N_1 N_2} \sum _{b_2 \in B_2} \left[ (p \, g_{b_2, (s_1,s_2)} - c^{\mathrm{proc}}) x_{b_2, (s_1,s_2)}^{\mathrm{proc}} - c^{\mathrm{ext}} x_{b_2, s_1}^{\mathrm{ext}} \right] \nonumber \\&\quad + \sum _{s_1 \in S_1, s_2 \in S_2, s_3 \in S_3} \frac{1}{N_1 N_2 N_3} \sum _{b_3 \in B_3} \left[ (p \, g_{b_3, (s_1,s_2,s_3)} - c^{\mathrm{proc}}) x_{b_3, (s_1,s_2,s_3)}^{\mathrm{proc}} - c^{\mathrm{ext}} x_{b_3, (s_1,s_2)}^{\mathrm{ext}} \right] . \end{aligned}$$

(16)