An adaptive large neighborhood search heuristic to optimize mineral value chains under metal and material type uncertainty

ABSTRACT

This paper addresses the optimisation of mineral value chains under metal and material type uncertainty. A mathematical model to simultaneously optimise the extraction decisions and the destination decisions is proposed. A fix-and-optimise scheme that exploits the structure of the problem and uses relaxation and decomposition techniques is introduced to obtain an initial solution, and an adaptive large neighbourhood search heuristic is developed to improve this solution. The proposed solution approach is tested on a real copper-gold mining complex. The results of these experiments show the ability of the proposed solution approach to efficiently address large instances of realistic size and provide schedules where the most valuable material is mined and processed early in the life of the mine and where the risk of not meeting production targets is successfully managed.

KEYWORDS:

• Heuristics
• adaptive large neighbourhood search
• mineral value chain
• mining complex
• metal uncertainty
• material type uncertainty
• risk management

Disclosure statement

No potential conflict of interest was reported by the author(s).

Annex A. Destroy methods

The 14 destroy methods used to select the $\beta$ blocks to be removed from the current solution, $x$, are described below. Unless otherwise specified, all destroy methods follow the general scheme summarised in Algorithm 2 and select blocks based on priority rules. More precisely, an index $p_i$, henceforth referred to as a priority value, is first calculated for each block $i$ using information derived from the current solution and/or from the history of the search. Blocks are then ranked in ascending or descending order of $p_i$, and the $\beta$ first blocks are selected. The way the priority values $p_i$ are computed differs from one method to another. Before describing how each method compute the $p_i$ values and how it ranks the blocks, some extra notation is introduced.

Let ${\alpha }_{is}=\sum _{d\in \mathcal{D}}{a}_{ids}$ be the number of destinations to which block $i$ can be sent without violating the admissibility constraints (recall that $a_{ids}=1$ if block $i$ is admissible for destination $d$ under scenario $s$ and 0 otherwise). Given the current solution $x$, denote by:

• $t_i\left(x\right)$: the period in which block $i$ is extracted in solution $x.$

• $E_i\left(x\right)=\underset{j\in \mathcal{P}\left(i\right)}{max}{t}_j\left(x\right):$ the earliest period in which block $i$ can be extracted without violating the slope constraints (recall that $\mathcal{P}\left(i\right)$ denotes the set of immediate predecessors of block $i$).

• $L_i\left(x\right)=\underset{j\in \mathcal{S}\left(i\right)}{min}{t}_j\left(x\right):$ the latest period in which block $i$ can be extracted without violating the slope constraints ($\mathcal{S}\left(i\right)$ denotes the set of immediate successors of block $i$).

• ${\mathcal{M}}_s^t\left(x\right)=\sum _{i\in \mathcal{N}}{w}_{is}{x}_i^t:$ the total tonnage, under scenario $s,$of blocks extracted in period $t$ in solution $x.$

• ${\mathcal{P}}_{ds}^t\left(x\right)=\sum _{i\in \mathcal{N}}{w}_{is}{y}_{ids}^t:$ the total tonnage, under scenario $s,$of blocks processed at destination $d$during period $t$ in solution $x.$

In what follows, in order to simplify the notation, the dependence on $x$ will be omitted whenever there is no risk of ambiguity; that is, $t_i$ will be used instead of $t_i\left(x\right)$ and so on.

A.1. Random picker (D1)

The main purpose of this method is to diversify the search. It simply selects at random $\beta$ blocks from the current solution $x$ to alleviate the risk of choosing the same blocks many times. This can be seen as associating with each block a random integer value $p_i$ chosen between 1 and $N$ and ranking the blocks in ascending order of $p_i$ (recall that $N$ denotes the number of blocks being scheduled).

A.2. Historical frequency (D2)

This method uses historical information to select the blocks. It relies on a frequency array $\mathcal{F}=({\mathcal{F}}_i)$ where each entry ${\mathcal{F}}_i$ is associated with a block $i$. This frequency array keeps track of the number of times that each block $i$ has been involved in destroying the solution since the beginning of the search process; that is, the ${\mathcal{F}}_i$ values are initially set to 0, and whenever the block $i$ is selected by a destroy method, then the value of the entry ${\mathcal{F}}_i$ is incremented. The priority values of the blocks are calculated as $p_i={\mathcal{F}}_i$, and the blocks are ranked in ascending order of $p_i$. Thus, this method selects the blocks less frequently chosen so far in order to diversify the search.

A.3. Historical best (D3)

This method is also based on historical information, but it additionally uses the value of the current solution, $f\left(x\right)$. It aims to select the blocks that seem to be sent to the wrong destinations in the current solution with regard to the best-known solutions. More precisely, let ${\mathcal{X}}_{it}$ be the set of solutions found so far in which block $i$ is extracted in period $t$. We define an $N\times T$ matrix $\mathcal{Z}$. The value ${\mathcal{Z}}_{it}$ of entry $\left(i,t\right)$ in this matrix corresponds to the value of the best solution in the set ${\mathcal{X}}_{it}$ (i.e. ${\mathcal{Z}}_{it}=\underset{sol\in {\mathcal{X}}_{it}}{max}f\left(sol\right)$). All ${\mathcal{Z}}_{it}$ values are initially set to a large negative value, and they are updated each time a new solution is found. Recall that $t_i$ denotes the period in which block $i$ is extracted in the current solution. The priority value of block $i$ is calculated as $p_i={\mathcal{Z}}_{i{t}_i}-f\left(x\right)$ . Hence, a positive value of $p_i$ means that in one of the solutions found in the past $\left(sol\right)$, block $i$ is extracted in the period as it is in the current solution ($x$); however, the value of $sol$ is better than the value of $x$. This might be because in the current solution $i$ is sent to the wrong destinations, and thus removing it from these destinations might result in an improvement. The blocks are ranked in descending order of $p_i$ to favour the blocks that present the largest deviations ${\mathcal{Z}}_{i{t}_i}-f\left(x\right)$.

A.4. Greedy picker (D4)

This method selects the costliest blocks in the current solution in an attempt to extract them in other periods where they will generate more profit and/or send them to better destinations. Identifying these blocks reduces to identifying blocks whose removal increases the value of the objective function the most. Let $f\left(x-\left\{i\right\}\right)$ denote the value of the current solution $x$ if block $i$ is removed from the schedule. The priority value of $i$ is calculated as the difference between this value and the value of the current solution; i.e. $p_i=f\left(x-\left\{i\right\}\right)-f\left(x\right)$. The blocks are ranked in decreasing order of $p_i$.

A.5. Period mobility (D5)

This method selects the blocks that can be extracted in other periods without violating the slope constraints, for it is easier to modify the period in which these blocks are extracted and thus create new feasible solutions different from the current one. Moreover, when selecting such blocks, we take care to favour blocks whose removal does not increase the mining shortage cost (third term of the objective function (14)). Let $\in \hspace{0.17em}$ be a small value ($\in$ = 0.0001 in the numerical results presented in Section 5). The priority value of block $i$ is calculated as follows:

$p_i=\frac{L_i-E_i}{\sum _{s}max\left(\in ,\underset{\_}{E_{t_i}}-\left(M_s^{t_i}-w_{is}\right)\right)}.$

Note that the numerator $\left(L_i-E_i\right)$represents the number of periods in which block $i$ can be extracted in the current solution without violating the slope constraints, while the denominator $(\sum _s\mathrm{m}\mathrm{a}\mathrm{x}\left(\in ,\underset{\_}{E_{t_i}}-\left(M_s^{t_i}-w_{is}\right)\right))$ is equal to a small value if removing $i$ from its current period $t_i$ does not incur a mining shortage in $t_i$ under any scenario, and it is equal to the total shortage amount considering all scenarios otherwise. Hence, a block with a high value of $p_i$ has more feasible reinsertion possibilities (in terms of periods of extraction) and is less likely to incur a mining shortage if removed from its current period. For this reason, the blocks are ranked in descending order of $p_i$.

A.6. Destination mobility (D6)

This method has the same objective as the previous one; that is, to select blocks that can lead to a new feasible solution different from the current one. However, the blocks are ranked by the number of destinations to which they can be sent rather than the number of periods in which they can be extracted. Recall that ${\alpha }_{is}=\sum _{d\in \mathcal{D}}{a}_{ids}$ represents the number of admissible destinations for block $i$ in scenario $s$. The priority values are calculated using the formula below, and the blocks are ranked in descending order of these values:

$p_i=\left\{\begin{array}{c}\sum _{s\in S}{\alpha }_{is}ifiisextractedinthecurrentsolution\left(\mathrm{i}.\mathrm{e}.,\mathrm{i}\mathrm{f}\sum _{t\in T}{x}_{it}=1\right),\\ 0,otherwise.\end{array}\right.$

Clearly, an unmined block cannot be sent to any destination, and thus if selected, one cannot modify its destination and get a new solution different from the current one. This is why unmined blocks are given less priority (the corresponding $p_i$ values are set to 0 to avoid selecting them).

A.7. Combined mobility (D7)

This method combines the two previous ones (Period mobility and Destination mobility). It accounts not only for the periods in which each block can be extracted, but also for the destinations to which the block can be sent. More precisely, the priority values are computed using the following formula:

$p_i=\left\{\begin{array}{c}\sum _{s\in n}\left({\alpha }_{is}-1\right)+\sum _{t=E_i,t\ne t_i}^{L_i}\sum _{s\in n}{\alpha }_{is}if\sum _{t\in n}{x}_{it}=1,\\ 0otherwise.\end{array}\right.$

The term ($\sum _{s\in \mathcal{S}}\left({\alpha }_{is}-1\right)$) accounts for the number of destinations to which block $i$ can be sent if its period of extraction is not modified (its current destination under each scenario does not contribute to this term), while the term ($\sum _{t=E_i,t\ne t_i}^{L_i}\sum _{s\in \mathcal{S}}{\alpha }_{is}$) considers the other periods to which the block can be moved while satisfying the slope constraints. The blocks are ranked in descending order of $p_i$ to select those with many feasible reinsertion possibilities.

A.8. Predecessor relatedness (D8)

This method has the same objective as the last three ones: create new feasible solutions different from the current one. However, it does not rely on the number of feasible reinsertion possibilities of each block nor does it follow the general scheme outlined in Algorithm 2. It selects blocks along with their predecessors extracted in the same period, as this should allow modifying their period of extraction and create a new feasible solution (more specifically, advance their extraction together which will ensure the satisfaction of the slope constraints). This is done in three steps. In the first step, a random integer value $\tau$ is chosen between 1 and $T$ (recall that $T$ is the number of periods over which the blocks are being scheduled). Then, $\tau$ periods are selected at random. The priority values are not calculated for all blocks but only for blocks extracted in one of these $\tau$ periods. Let $t$ be such a period and $i$ be a block extracted in $t$. Denote by ${\gamma }_i$ the number of blocks in the inverted cone formed by $i$ and its predecessors extracted in $t$. Recall that $\beta$ blocks should be selected. We set $p_i=\left|{\gamma }_i-\frac{\beta }{\tau }\right|$ and among the blocks currently extracted in $t$, the block with the smallest value of $p_i$ is selected, as well as its predecessors extracted in $t$. Ties are broken up randomly. This process is repeated for each of the $\tau$ periods.

A.9. Successor relatedness (D9)

This method is similar to the previous one except that it selects blocks along with their successors. Apart from the fact that ${\gamma }_i$ represents the number of blocks in the cone formed by $i$ and its successors extracted in the same period, the procedure to select the blocks is identical to the procedure described in the previous section (A.8).

A.10. Mining reduction (D10)

The objective of this method is to select blocks whose removal can reduce the mining surplus (i.e. decrease the value of the fourth term of the objective function) or reduce the tightness of the mining capacity constraints to make room for new blocks. The priority values are computed as follows:

$p_i=\left\{\begin{array}{c}min\left(1,L_i{-}_i\right)\underset{s\in n}{max}\frac{n_s^{t_i}}{\overline{E^{t_i}}}if\sum _{t\in n}{x}_{it}=1,\\ 0otherwise.\end{array}\right.$

The blocks are ranked in descending order of$p_i$. Note that if an extracted block $i$cannot be moved to another period ($t\ne t_i$) without violating the slope constraints (i.e. if $\sum _{t\in \mathcal{T}}{x}_{it}=1$ and $E_i=L_i=t_i$), then the corresponding priority value $p_i$is equal to 0 to avoid selecting it. This is done to make it easy for the repair method to create new feasible solutions. The priority value of a block that is not extracted in the current solution (i.e. such that $\sum _{t\in \mathcal{T}}{x}_{it}=0$) is also set to $p_i=0$ to avoid selecting it, as such blocks are not extracted and thus have no influence on the mining capacity constraints (changing their periods does not affect the fourth term of the objective function (14)).

A.11. Processing reduction (D11)

This method is based on similar ideas as the previous one and aims to reduce the amount of surplus at the different destinations and/or the tightness of the processing capacity constraints. To be more precise, recall that ${\mathcal{P}}_{ds}^t$ denotes the total tonnage of blocks processed at destination $d$ during period $t$ under scenario $s$ in the current solution, and that $\overline{F_d^t}$ is the processing capacity at $d$ during $t$. Denote by ${\sigma }_{is}=\mathrm{m}\mathrm{i}\mathrm{n}\left(1,\left({\alpha }_{is}-1\right)+{\alpha }_{is}\left(L_i-E_i\right)\right)$ the number of feasible reinsertion possibilities for block $i$ under scenario $s$. Thus, if block $i$ cannot be moved to another period (i.e. if $E_i=L_i=t_i$), and if it can be sent to only one destination under scenario $s$ (i.e. and if ${\alpha }_{is}=1$), then ${\sigma }_{is}=0$. Otherwise, ${\sigma }_{is}\ge 1$. The priority values $p_i$ are calculated using the formula below, and the blocks are ranked in descending order of $p_i$:

$p_i=\left\{\begin{array}{c}\underset{s\in \mathcal{S}}{max}\left\{{\sigma }_{is}\frac{{\mathcal{P}}_{ds}^{t_i}}{\overline{F_d^{t_i}}}\right\}if\sum _{t\in \mathcal{T}}{x}_{it}=1,\\ 0otherwise.\end{array}\right.$

Note that the way the $p_i$’s are defined implies that the blocks that are not extracted in the current solution, as well as the blocks that are extracted but can neither be moved to a different period nor sent to another destination for all scenarios, have the lowest priority values (0) to avoid selecting them. As for the previous method, the idea is to prevent getting a solution similar to the current one when applying the repair method.

A.12. Shortage cautious (D12)

The purpose of this method is to find blocks that can be sent to destinations more profitable than their current destinations without incurring a shortage in their current destination. Denote by $c$ the destination to which block $i$ is sent under scenario $s$ in the current solution. Let $d^{\ast }$ be the best destination to which block $i$ can be sent under scenario s; i.e. $d^{\ast }=argma{x}_{d\in \mathcal{D}}{v}_{ids}$(recall that $v_{ids}$ represents the economic value to be generated if block $i$ is processed at destination $d$ in scenario $s$. This value is calculated as the return from selling the recovered metal minus the processing, transportation, and selling costs. $v_{ids}$ is set to a large negative value if $i$ is not admissible for destination $d$ under scenario $s$). Again, let $\in$ be a small positive value, and let $C$ be a large positive value. Recall that ${\delta }_1$ is the economic discount rate. The priority values are computed using the following formula, and the blocks are ranked in descending order of the priority values:

$p_i=\left\{\begin{array}{c}\frac{1}{{\left(1+{\delta }_1\right)}^{t_i}}\sum _{s\in n}\frac{v_{i{d}^{\ast }s}-v_{ics}}{\mathrm{m}\mathrm{a}\mathrm{x}\left(\in ,\underset{\_}{F_c^{t_i}}-\left(n_{cs}^{t_i}-w_{is}\right)\right)}if\sum _{t\in n}{x}_{it}=1,\\ -Cotherwise.\end{array}\right.$

In this formula, $\frac{1}{{\left(1+{\delta }_1\right)}^{t_i}}$is used to account for the discount factor. The numerator ($v_{i{d}^{\ast }s}-v_{ics}$) is used to favour blocks that can improve the second term of the objective function (14) the most. The denominator ($\mathrm{m}\mathrm{a}\mathrm{x}\left(\in ,\underset{\_}{F_c^{t_i}}-\left({\mathcal{P}}_{cs}^{t_i}-w_{is}\right)\right)$)) is used to favour blocks that will not incur a shortage if removed from their current destination (will not increase the fifth term of the objective function (14)). Finally, the priority values of blocks that are not extracted in the current solution are set to a large negative value to avoid selecting them, as these blocks do not contribute to any term of the objective function (14).

A.13. Empty one period (D13)

This method does not compute the priority values and does not necessarily select $\beta$ blocks. It randomly selects one period and adds all the blocks extracted in that period to the list $\mathcal{L}$ (list of selected blocks). The motivation is to allow all destinations in a given period to be empty and completely remake the destination decisions (reconstruct the destination plans) with the repair method. By doing so, more opportunities for new block combinations in the different destinations are created.

A.14. Empty waste dump (D14)

In a given period and under a given scenario, some blocks might be sent to the waste dump while they can be sent to profitable destinations where they can be processed and generate revenue. This method aims to select such blocks to improve the solution. Let ${\pi }_{is}$ be a parameter equal to 1 if block $i$ is sent to the waste dump under scenario $s$ in the current solution, and 0 otherwise. Clearly, ${\pi }_{is}=0\mathrm{\forall }s$ if $i$is not extracted in the current solution. The priority values are computed using the formula below, and the blocks are ranked in descending order of these values:

$p_i=\sum _{s\in \mathcal{S}}{\pi }_{is}\left[\left({\alpha }_{is}-1\right)\left(L_i-E_i+1\right)\right]$

The term $\left({\alpha }_{is}-1\right)$ accounts for the number of the destinations to which block $i$ can be sent under scenario $s$, excluding its current destination, while $\left(L_i-E_i+1\right)$ accounts for the periods in which $i$ can be extracted without violating the slope constraints, including its current period of extraction. The factor ${\pi }_{is}$ is used to avoid selecting blocks that are currently not in the waste dump under any scenario.

Annex B. Repair methods

Referring to Algorithms 1 and 2, the blocks in the list $\mathcal{L}$, identified by the selected destroy method, are removed from the current solution $x$, resulting in an infeasible solution $x^{-}$. This means that all the variables associated with the blocks that are not in $\mathcal{L}$ are fixed, and the remaining variables are ‘free’. One of the seven methods described below is used to reinsert each block in $\mathcal{L}$ in other feasible periods and/or destinations to obtain a new feasible solution $x{ }^{\mathrm{\prime }};$that is, to optimise the ‘free’ variables. The following notation is used. Given the solution to repair $x^{-}$, ${\mathcal{B}}^t=\left\{i\in \mathcal{N}:x_i^t=1\right\}$ denotes the set of blocks extracted in period $t$. The set of blocks processed in destination $d$ during period $t$ under scenario $s$ is denoted by ${\mathrm{\Lambda }}_{ds}^t=\left\{i\in \mathcal{N}:y_{ids}^t=1\right\}.$

B.1. Random repair (R1)

This method considers one block $i\in \mathcal{L}$ at a time, which is selected randomly, and sequentially chooses the period and the destinations in which $i$ will be scheduled. This is done as follows: The method starts by identifying the set of feasible periods $\mathcal{F}\mathcal{P}\left(i\right)$ in which $i$ can be extracted without violating the slope constraints. In doing so, the predecessors and the successors of $i$ that are in the list $\mathcal{L}$ are not accounted for. Then, one of the periods in $\mathcal{F}\mathcal{P}\left(i\right)$ is selected randomly, and $i$ is scheduled to be extracted in that period. The next step is to decide in which destination $i$ will be processed under each scenario, and again this is done randomly; i.e. under each scenario, $i$ can be processed in any destination as long as it is an admissible destination. When all scenarios are considered, $i$ is removed from $\mathcal{L}$, another block is chosen, and the process is repeated until the list $\mathcal{L}$ is empty.

B.2. Greedy repair (R2)

This method is similar to the previous one in the sense that it considers blocks $i\in \mathcal{L}$ one at a time and sequentially determines the period and the destinations for the selected block before considering another block, but, to this end, it uses selection criteria different from those used by R1 as explained below. To simplify the presentation, we denote by

(29) $g\left({\mathcal{B}}^t\right)=\sum _{i\in {\mathcal{B}}^t}\frac{E\left[c_i\right]}{{\left(1+{\delta }_1\right)}^t}+\frac{1}{S}\sum _{s\in \mathcal{S}}\left[\frac{p^{-}}{{\left(1+{\delta }_2\right)}^t}\mathrm{m}\mathrm{a}\mathrm{x}(\underset{\_}{E^t}-\sum _{i\in {\mathcal{B}}^t}{w}_{is},0)+\frac{p^{+}}{{\left(1+{\delta }_2\right)}^t}\mathrm{m}\mathrm{a}\mathrm{x}(\sum _{i\in {\mathcal{B}}^t}{w}_{is}-\overline{E^t},0)\right]$(29)

(30) $h\left({\mathrm{\Lambda }}_{ds}^t\right)=\sum _{i\in {\mathrm{\Lambda }}_{ds}^t}\frac{v_{ids}}{{\left(1+{\delta }_1\right)}^t}-\frac{q_d^{-}}{{\left(1+{\delta }_2\right)}^t}\mathrm{m}\mathrm{a}\mathrm{x}(\underset{\_}{F_d^t}-\sum _{i\in {\mathrm{\Lambda }}_{ds}^t}{w}_{is},0)-\frac{q_d^{+}}{{\left(1+{\delta }_2\right)}^t}\mathrm{m}\mathrm{a}\mathrm{x}(\sum _{i\in {\mathrm{\Lambda }}_{ds}^t}{w}_{is}-\overline{F_d^t},0).$(30)

$g\left({\mathcal{B}}^t\right)$ is used to evaluate the cost of extracting a block $i$ in period $t$ accounting for all blocks that are already extracted in this period, whereas $h\left({\mathrm{\Lambda }}_{ds}^t\right)$ is used to measure how profitable it is in scenario $s$ and period $t$ to process an additional block in destination $d$ accounting for blocks that are already processed in $d$. Again, let $\mathcal{F}\mathcal{P}\left(i\right)$ denote the set of feasible periods in which $i$ can be extracted, considering only its predecessors and successors that are not in $\mathcal{L}$. For each period $t\in \mathcal{F}\mathcal{P}\left(i\right)$, the repair method R2 first uses function (29) to compute the cost of extracting $i$ in period $t$: ${\mathrm{\Lambda }}_1\left(i,t\right)=g\left({\mathcal{B}}^t{\cup }^{\left\{i\right\}}\right)-g\left({\mathcal{B}}^t\right)$. Then, it considers the scenarios sequentially and for each scenario, it finds, following a greedy approach, the destination in which $i$ can be processed. For that, function (30) is used and the following is computed to measure how feasible and profitable it is to process $i$ in $d$: ${\mathrm{\Delta }}_2\left(i,d,s,t\right)=h\left({\mathrm{\Delta }}_{ds}^t{\cup }^{\left\{i\right\}}\right)-h\left({\mathrm{\Lambda }}_{ds}^t\right)$ if $i$ is admissible to $d$ under scenario $s$, and ${\mathrm{\Lambda }}_2\left(i,d,s,t\right)$ is set equal to a large negative value otherwise.

Let $d^{\ast }\left(i,s,t\right)=arg\underset{d\in \mathcal{D}}{max}{\mathrm{\Delta }}_2\left(i,d,s,t\right)$. Clearly, the maximum profit that one can expect if block $i$ is extracted in period $t$ is $\mathrm{\Delta }\left(i,t\right)=-{\mathrm{\Delta }}_1\left(i,t\right)+\frac{1}{S}\sum _{s\in \mathcal{S}}{\mathrm{\Delta }}_2\left(i,d^{\ast }\left(i,s,t\right),s,t\right)$. Once all periods in $\mathcal{F}\mathcal{P}\left(i\right)$ are considered, the period $t^{\ast }=argma{x}_{t\in \mathcal{F}\mathcal{P}\left(i\right)}\mathrm{\Delta }\left(i,t\right)$ is identified and block $i$ is scheduled to be extracted in $t^{\ast }$ (i.e. it is included in the set ${\mathcal{B}}^t$). Finally, for each scenario $s$, $i$ is sent to destination $d^{\ast }\left(i,s,t\right)$ (i.e. it is included in the set ${\mathrm{\Lambda }}_{d^{\ast }\left(i,s,t\right)s}^t$).

B.3. Capacity cautious repair (R3)

This method is very similar to the previous one except for the following differences:

1. $g\left({\mathcal{B}}^t\right)$ and $h\left({\mathrm{\Lambda }}_{ds}^t\right)$ are here defined by equations (31)(31) $g\left({\mathcal{B}}^t\right)=\frac{{\sum }_{s\in \mathcal{S}}{\sum }_{i\in {\mathcal{B}}^t}{w}_{is}}{\overline{E^t}}$(31) and (32(32) $h\left({\mathrm{\Lambda }}_{ds}^t\right)=\frac{{\sum }_{i\in {\mathrm{\Lambda }}_{ds}^t}{w}_{is}}{\overline{F_d^t}}.$(32) ) below rather than (29) and (30) in order to select the periods/destinations that will leave the most residual capacity for forthcoming blocks. This is done to introduce some lookahead perspective when reinserting the blocks.

(31) $g\left({\mathcal{B}}^t\right)=\frac{{\sum }_{s\in \mathcal{S}}{\sum }_{i\in {\mathcal{B}}^t}{w}_{is}}{\overline{E^t}}$(31)

(32) $h\left({\mathrm{\Lambda }}_{ds}^t\right)=\frac{{\sum }_{i\in {\mathrm{\Lambda }}_{ds}^t}{w}_{is}}{\overline{F_d^t}}.$(32)

(2) Accordingly, $d^{\ast }\left(i,s,t\right)=arg\underset{d\in \mathcal{D}}{min}{\mathrm{\Delta }}_2\left(i,d,s,t\right)$ instead of $d^{\ast }\left(i,s,t\right)=arg\underset{d\in \mathcal{D}}{max}{\mathrm{\Delta }}_2\left(i,d,s,t\right)$ (i.e. the best destination for a block in a given scenario and period is the one having the smallest capacity utilisation).

(3) $\mathrm{\Delta }\left(i,t\right)$is set equal to ${\mathrm{\Delta }}_1\left(i,t\right)+\frac{1}{S}\sum _{s\in \mathcal{S}}{\mathrm{\nabla }}_2\left(i,d^{\ast }\left(i,s,t\right),s,t\right)$ and $t^{\ast }=argmi{n}_{t\in \mathcal{F}\mathcal{P}\left(i\right)}\mathrm{\Delta }\left(i,t\right).$

B.4. MCFP repair (R4 and R5)

To repair the solution, this method combines the random repair heuristic (R1) and the MCFP heuristic described in Section 3.2. More specifically, it starts by assigning feasible periods to blocks in $\mathcal{L}$ as in R1; that is, for each block $i\in \mathcal{L}$, $\mathcal{F}\mathcal{P}\left(i\right)$ is first identified, then $i$ is included in ${\mathcal{B}}^t$ where $t$ is chosen randomly in $\mathcal{F}\mathcal{P}\left(i\right)$. Once this step is completed, the destination plans for each period $t$ that have been affected at the previous step are determined by solving the $D{P}^t$ described in 3.2. This is done by applying the MCFP heuristic on each scenario separately. When applying the MCFP heuristic, the decisions associated with the blocks that were not in $\mathcal{L}$ (i.e. blocks that were not selected by the destroy method) are fixed to their current values. Another alternative, which is more flexible and might lead to better quality solutions but at the expense of longer computational times, is to reconstruct the destination plan from scratch (i.e. none of the parts of the plan is fixed and all the decisions are to optimise). In this paper, we examine the two alternatives, which leads to two variants of the MCFP repair method. The variant that solves a partial destination problem (first alternative) is denoted by R4, while the variant that solves the full destination problem (second alternative) is denoted by R5.

B.5. MIP repair (R6 and R7)

This method is very similar to the previous one. All the extraction decisions are made first before designing the destination plans. Again, the latter are determined by considering only periods that have been affected when making the extraction decisions, considering the scenarios separately. However, rather than using MCFP, a mixed-integer programming solver is used to find the optimal values of the variables $y_{ids}^t$, $f_{ds}^{t-}$, and $f_{ds}^{t+}$ that maximise the net present economic value to be generated from processing the blocks extracted in period $t$ minus the total penalty costs of not satisfying the demands or exceeding the capacities of the different destinations during this period. Again, one can fix the binary variables $y_{ids}^t$ corresponding to the blocks $i$ that were not selected by the destroy method to their current values, which results in a variant of the MIP repair method that we will denote by R6 as one can ‘free’ all variables, which gives yield to another variant of the MIP repair method denoted by R7.

Additional information

Funding

The work in this paper was funded by NSERC CRDPJ 411270-10, NSERC Discovery Grant 239019-06, and the industry members of the COSMO Stochastic Mine Planning Laboratory: AngloGold Ashanti, Barrick, BHP, De Beers, Kinross, Newmont, and Vale. This support is gratefully acknowledged.