Scenario Reduction Applied to Geostatistical Simulations

Abstract

Having a large number of geostatistical simulations of a mineral or petroleum deposit provides a better idea of its upside potential and downside risk; however, large numbers of simulated realizations of a deposit may pose computational difficulties in subsequent decision-making phases. Hence, depending on the specific case, there can be a need to select a representative subset of conditionally simulated deposit realizations. This paper examines and extends an approach developed by the stochastic optimization community based on stochastic mathematical programming with recourse and is discussed here in the context of mineral deposits while it is possibly suitable for other earth science applications. The approach is based on measuring the “distance” between simulations and the introduced distance measure between simulated realizations of a mineral deposit is based on the metal above a given set of cutoff grades while a pre-existing mine design is available. The approach is tested on 100 simulations of the Walker Lake data with promising results.

Appendices

Appendix A: Multistage Programming with Recourse

Multistage programming with recourse (Sahinidis 2004; van der Vlerk 2007; Beasley 2010) is an optimization method that allows for the consideration of different mining operational options in order to optimize profits, such as stockpiling low grade ore, changing a cutoff grade strategy, and so on. Some variables such as the commodity price and the exchange rate are outside the control of decision-makers. The former are called state variables while the latter are called control variables. As the commodity price fluctuates, the decision-makers change the values of the control variables (for example, the quantity sent to the stockpile or taken from it) so as to optimize an objective function. More formally, at the outset the decision-maker sets the initial value of the control variable x ₀ knowing the initial value ξ ₀ of the state variable. At each time period, the new value of the control variable x _i is selected once the new ξ _i is known, so as to optimize an objective function.

The procedure is

(3)

Kall and Mayer (2005) show that the associated dynamic programming equations can be written as a set of nested conditional expectations that are evaluated starting from the last time period and working backward

(4)

where f denotes the objective function. The same authors comment that even if the terms to be optimized are linear, the resulting multistage optimization problem is nonlinear if the state variables have a continuous distribution. However, if their distribution is discrete, then the possible values of the state variables (ξ ₁,ξ ₂,…,ξ _T ) can be represented by a finite set of scenarios, each with a certain probability, and the system can be solved by linear programming. This is why branching trees are used to model the variables.

Appendix B: Theory on Stability

Although many algorithms have been developed specifically for solving large scale problems by exploiting specific features, there comes a time when it is no longer possible to solve the problem numerically. This is referred to as the curse of dimensionality. It is therefore natural to ask whether it is possible to reduce the size of the problem by, for example, reducing the number of scenarios. The approach that described here is based on work by Dupacova et al. (2003) and Heitsch and Romisch (2003).

Suppose that the discretized scenario tree contains N scenarios $\xi^{k} = \{ \xi_{1}^{k}, \ldots,\xi_{T}^{k} \}$, k=1,…,N, each with a probability p _k for k=1,…,N. Let P denote this discrete distribution P(ξ=ξ ^k)=p _k . The principle behind the scenario reduction is to find the subset of scenarios of a predetermined size n (n<N) and to attribute new probabilities to these scenarios so that the new probability measure Q is as close as possible to the original one P in the sense of a certain probability distance called a canonical or ideal probability metric. To start with, a general result on stability is reviewed, similar to Dupacova et al. (2003).

2.1 B.1 Key Stability Result

We wish to minimize

$$ \mathop{\mathrm{Min}} _{{x} \in\mathcal {X}}\mathbb{E}_{{P}}f ( \xi,x ) = \int_{\varXi} f ( \xi,x ) {P} ( d\xi ). $$

(5)

Here, $\mathcal{X} \in \mathbb{R}^{m}$ is a closed non-empty convex set, Ξ is a closed subset of ℝ^s; $\mathfrak{B}$ is the σ algebra of the Borelians on Ξ. The function $f: \varXi\times R^{\mathrm{m}} \to\overline{\mathbb{R}}$ where $\overline{\mathbb{R}}$ is the right half axis is continuous relative in ξ and convex in x. P is a fixed probability measure on $( \varXi,\mathfrak{B} )$ and $\mathbb{E}_{{P}}$ denotes the expectation relative to P. A multistage problem with recourse can be rewritten in this form. In convex stochastic programming, the functions f(⋅,x), $x \in \mathcal{X}$ are not differentiable but have the following property: there exist a continuous non-decreasing function h:ℝ₊→ℝ₊ with h(0)=0 and a function $g : \mathbb{R}_{ +} \to \mathbb{R}_{ +}^{*}$ and an element ξ ₀∈ℝ^s such that

$$ \bigl| f( \xi,x ) - f( \tilde{\xi} ,x ) \bigr| \le g\bigl( \| x \| \bigr)c( \xi,\tilde{\xi} ), $$

(6)

for each $x \in \mathcal{X}$ and where the function c:Ξ×Ξ→ℝ₊ is defined by

$$ c ( \xi,\tilde{\xi} ) : = \max \bigl\{ 1,h \bigl( \Vert \xi- \xi_{0} \Vert \bigr),h \bigl( \Vert \tilde{\xi} - \xi_{0} \Vert \bigr) \bigr\} \Vert \xi- \tilde{\xi} \Vert , \quad \forall\xi,\ \tilde{\xi} \in\varXi $$

(7)

for a certain metric ∥⋅∥ on ℝ^s. Typical examples of the function h for applications are h≡1 and h(r)=r ^p−1 for r∈ℝ₊ with p≥1 (Dupacova et al. 2003 and Heitsch and Romisch 2003). The following sets are introduced

(8)

for each ε≥0. The set of solutions S(P) is by definition S(P):=S ₀(P).

For two types of optimization problem on ℝ^d, represented by two functions f and g, where the second problem can be considered as an approximation of the first, it is possible to study the quantitative estimates of the closeness of infg and argming with inff and argminf in terms of an epi-distance (Rockafellar and Wets 2004). In this case there are two probabilities: P as in Eq. (5), and an approximation to it called Q. The aim is to estimate the closeness of v(Q) and S _ε (Q), to v(P) and S _ε (P). This will be carried out using a probability metric defined on the appropriate space.

Using estimates of the approximation to optimality via an epi-distance between $\mathbb{E}_{P}f( \xi, \cdot)$ and $\mathbb{E}_{Q}f( \xi, \cdot)$ found in Rockafellar and Wets (2004) and Dupacova et al. (2003) showed that the model of Eq. (5) is stable relative the probability metric

$$ \zeta_{c} ( P,Q ) : = \sup _{f \in \mathcal{F}_{c} }\biggl \vert \int_{\varXi} f ( \xi, {x} ) P ( {d}\xi ) - \int_{\varXi} f ( \xi, {x} ) Q ( {d}\xi ) \biggr \vert , $$

(9)

where $\mathcal{F}_{c}$ is the class of continuous functions defined by

$$ \mathcal{F}_{c} = \bigl\{ f : \varXi\to \mathbb{R}\mbox{ such that }f ( \xi ) - f ( \tilde{\xi} ) \le c ( \xi,\tilde{\xi} )\ \forall\xi,\ \tilde{ \xi} \in\varXi \bigr\}. $$

(10)

The probabilities P and Q belong to the set $\mathcal{P}_{c}( \varXi)$ defined by

(11)

The probability metric ζ _c is called a Fortet–Mourier (type) metric and is defined on $\mathcal{P}_{c}( \varXi)$. Dupacova et al. (2003) give a similar result in a more general context. This research presents how to apply this result to the scenario reduction problem.

For two probability measures P ₁ and P ₂ on a separable metric space (U,d), let $\mathcal{P}^{( P_{1},P_{2} )}$ be the space of Borel probability measures P on U×U (that is defined on σ-algebra $\mathcal{B}( U \times U )$ of Borel subsets of U×U) with fixed marginal distributions P ₁(⋅)=P(⋅×U) and P ₂(⋅)=P(U×⋅). Let $\mathcal{L}^{{( P_{1},P_{2} )}}$ be the set of finite Borel measures Q on $\mathcal{B}( U \times U )$ such that Q(A×U)−Q(U×A)=P ₁(A)−P ₂(A) for all $A \in \mathcal{B}( U )$.

Putting $\mu_{\tilde{c}}( P ) : = \int_{U \times U} \tilde{c}( x,y )P( dx,dy )$, where $\tilde{c}( x,y )$ is a positive continuous function defined on U×U, the following problems are called the Kantorovich problem and the Kantorovich–Rubinstein problem, respectively

(12)

(13)

The solution to Eq. (12) is called the Kantorovich functional, similarly the solution to Eq. (13) is called the Kantorovich–Rubinstein functional. Rachev (1991) showed that the distance ζ _c in Eq. (5) can be controlled via the Kantorovich–Rubinstein functional

$$ \zeta_{c} ( P,Q ) = \mathring{\mu} _{c} ( P,Q ) \quad \mbox{where }P,Q \in \mathcal{P}_{c} ( \varOmega ), $$

(14)

the Kantorovich functional

$$ \zeta_{c} ( P,Q ) \le\hat{\mu}_{c} ( P,Q ) \quad \mbox{where } P,Q \in \mathcal{P}_{c} ( \varOmega ). $$

(15)

The inequality in Eq. (15) becomes an equality if h≡1. So one of these functionals can be used to approximate the probability P of the original scenarios since in the case where the probability measures P and Q are discrete with a finite number of atoms, problems Eqs. (12) and (13) become finite-dimensional linear programming problems. For example, if P is a discrete probability measure with a finite number of scenarios (ξ ¹,…,ξ ^N) with probability (p ₁,…,p _N ) and if Q is another discrete probability with a finite number of scenarios $( \tilde{\xi}^{1}, \ldots,\tilde{\xi}^{M} )$ with probability (q ₁,…,q _M ), then the Kantorovich functional can be written as

$$ \hat{\mu} ( P,Q ) = \mathrm{min}\sum _{i = 1}^{N}\, \sum_{j = 1}^{M} c \bigl( \xi^{i},\tilde{\xi}^{j} \bigr) \eta_{ij}, $$

(16)

subject to the constraints

(17)

In other words, $\hat{\mu} ( P,Q )$ represents the optimal value of a linear transport problem.

2.2 B.2 Scenario Reduction

Dupacova et al. (2003) proposed using the functional $\hat{\mu}_{c}$ to reduce the number of scenarios. It is assumed that the process $\xi= \{ \xi_{t} \}_{t = 1}^{T}$ has been discretized as a scenario tree $\xi^{i} = \{ \xi_{1}^{i},\xi_{2}^{i}, \ldots,\xi_{T}^{i} \}$, i=1,…,N, where N is the total number of scenarios and p _i are the probabilities with $\sum_{i = 1}^{N} p_{i} = 1$. Let P be the distribution of this process: $P = \sum_{i = 1}^{N} p_{i}\delta_{\xi ^{i}}$. Consider a new probability measure Q whose support is composed of n scenarios ξ ^j with the probabilities q _j . Let J⊂{1,…,N} be the index set of scenarios to be eliminated. The probabilities associated with the scenarios ξ ^j, j∈J that are eliminated are attributed to the scenarios that are kept, using the redistribution rule described below. The probability measure Q can be written as $Q = \sum_{j \notin J} q_{j}\delta_{\xi^{j}}$. Let D(J;q) be the distance between the two probability measures defined by

$$ D( J;q ) : = \hat{\mu}_{c}( P,Q ) = \hat{ \mu}_{c}\Biggl(\, \sum_{i = 1}^{N} p_{i}\delta_{\xi^{i}} ,\sum_{j \notin J} q_{j}\delta_{\xi^{j}} \Biggr) , $$

(18)

where $\hat{\mu}_{c}$ is the probability distance given in Eq. (15).

The principle behind optimal scenario reduction can be written as

(19)

Thus, an optimal set of indices J _∗ is sought, together with the optimal probabilities q _∗ such that D(J _∗;q _∗)=min _J,q {D(J;q)} where the minimization is carried out over all possible setsJ of a predetermined size. For a fixed set of indices, it is possible to calculate min _q {D(J;q)} and find the optimal probabilities to assign to each of the scenarios that have been kept. Dupacova et al. (2003) proved that this can be done by redistributing the probability of each of the scenarios that are eliminated to the closest scenario amongst those retained. This is called the optimal redistribution rule.

Appendix C: Testing the Heitsch and Romisch’s Sequential Selection Method

Multistage stochastic programming with recourse is widely used for optimizing projects subject to uncertainty, such as power stations. The first step is to construct a set of branching scenarios that model the evolution of the key variable(s) over time. As this research is interested in evaluating and optimizing projects, a discrete tree of the commodity price over time is constructed. Figure 1 shows a three-stage binomial (Hull 2009 and Derman et al. 1996). In finance, recombining trees are used, and so there are (N+1) nodes at the Nth stage. In multistage stochastic programming, the branches do not recombine (because earlier decisions have lasting effects) so the number of nodes increases exponentially.

Heitsch and Romisch (2003) describe the forward and backward scenario reduction algorithms in detail and give examples using data on electrical load and water inflow into dams provided by the French utility EDF. Dupacova et al. (2003) carried out numerical tests on electric load scenario trees for power management under uncertainty, and found that after 50 % reduction of the scenario tree the optimal reduced tree has still about 90 % relative accuracy. Given the quality of their work, the accuracy of their result is undoubtedly correct. But as their full tree has more than 16 million branches (2²⁴), this research sought to test the two algorithms on a much smaller example, one where all possible combinations could be exhaustively computed. To do this, a classical binomial tree for a geometric Brownian motion with 4 time periods and hence 16 terminal nodes is chosen. The parameters values (22 % volatility; 5 % risk-free rate) were taken from an example in Hull (2009) because this is the standard textbook in quantitative finance.

The decision is made to reduce the tree to one with 6 branches. As there are only 8008 combinations of 6 out of a total of 16, it is easy to compute the metric D(J;q) for all combinations and hence compare the best one with those for forward and backward algorithms. Figure 5 shows the histogram of the values of D(J;q) for all 8008 combinations. These range from 283.5 to about 2000.

Fig. 5

The histogram of the values of D(J;q) for all 8008 combinations. These range from 283.5 (optimum) to over 2000

Table 4 shows the values of D(J;q) for both algorithms and for the optimum, together with their rank out of 8008. The six scenario numbers selected are also given. The backward algorithm performed better than the forward one but its rank position is only 1516 out of 8008 and its value of D(J;q) is about 50 % greater than the best one. It is surprising to see how far from optimal their algorithms were in this case. These results suggest that an iterative random search should perform at least as well.

Table 4 Comparing the backward and forward algorithms proposed by Heitsch and Romisch with the optimum computed exhaustively from all possible combinations. The values of D(J;q) are given in the top row while their ranks are in the second row. The last row shows the scenario numbers retained in each case

Appendix D: Definitions of Recoverable Reserves

Let Q(i,k) and T(i,k) denote the metal and the tonnage above the kth cutoff grade cutoff, Z _C (k), in ith panel. The panel contains 100 small blocks. Let I _Z(i,l)>ZC(k) be the indicator function that takes the value 1 if the grade of the lth small block in the ith panel is above the cutoff Z _C (k) and 0 otherwise. These are defined as follows

Matheron (1975) defined the conventional profit as

$$B ( i,k ) = Q ( i,k ) - ZC \times T ( i,k ). $$