Integrated stochastic optimisation of stope design and long-term production scheduling at an operating underground copper mine

ABSTRACT

Conventional underground long-term mine planning is based on a deterministic stepwise framework, which is unable to effectively manage the synergies between the mine planning components or to manage the orebody risk in production schedules and forecasts. The present study enhances prior integrated optimisation approaches to a variation of the sublevel longhole open stoping mining method using backfilling, through a new two-stage stochastic integer programming formulation. A comparison to a sequential stochastic approach shows different stope layouts and extraction sequences for a copper mine with secondary elements being gold and uranium. The integrated approach shows lower horizontal development costs and 6% higher net present value.

KEYWORDS:

• Underground mining
• mine planning
• integrated stochastic optimisation
• geological uncertainty
• sublevel stoping
• mathematical programming application

1. Introduction

Long-term underground mine production scheduling conventionally consists of the definition of the stope layout, the access design, the ventilation system requirements, and the extraction sequence [1–4]. Initially, potential mineable volumes are defined based on geomechanical and geological properties. This step requires fixed input values for operational parameters, such as mining capacity, and a cut-off grade to maximise the undiscounted profit of the stopes [5–10]. Subsequently, accesses such as declines, ramps, and shafts connecting production areas to main haulage and ventilation systems are defined [11–13]. The stope layout and network designs are used as inputs for the life-of-mine (LOM) production schedule optimisation, which maximises the net present value (NPV) of the mining asset under a unit operation timeframe, while considering economic and capacity constraints [14–22]. An iterative procedure ensures that these steps are repeated for different cut-off grade values until the one that produces a schedule with the highest NPV is found and selected [3,4,6]. However, by considering separately the stope design and extraction sequence, the synergies between these two components are not captured in the optimisation process. In addition, this traditional approach is deterministic in which a single estimated (average type) representation of the orebody is considered, which does not quantify geological uncertainty and variability, representing a critical source of technical risk for a mining project [23].

Little, et al. [18] show that stope boundaries should be an outcome of the production schedule so that interdependencies among stope grades, development costs, and the time value of money are incorporated into an integrated process. Accordingly, Little, et al. [18], and Little, et al. [24] propose an integrated underground mine design and production scheduling optimisation approach, based on a mixed integer programming (MIP) formulation, aiming to maximise the discounted cashflow of stopes mined in an extraction sequence accounting for their size and location, while constrained by ore production and backfilling capacities. This approach does not account for critical costs and decisions associated with the development of main and secondary accesses. Copland and Nehring [25] address stope boundaries and scheduling with a binary integer programming formulation that considers the discounted profit from stopes and incorporates level access development decisions. The development of drifts and shafts as one-time decisions is integrated into the MIP by Hou, et al. [26]. Nevertheless, it assumes a predefined network design in terms of the position of the longitudinal ore drives and the shaft. In addition, the application of the proposed method is restricted to a stratiform two-dimensional deposit. Foroughi, et al. [27] develop a MIP that maximises the NPV and metal recovery, while jointly considering the production scheduling and stope layout. The above-mentioned method is applied to an iron ore deposit and is solved using a genetic algorithm. The development and application of the methods delineated thus far are suited to the sublevel open-stoping (SLOS) mining method. The development of ventilation systems and the internal accessibility of stopes are not thoroughly investigated in the studies previously presented, not allowing the generalisation to other mining methods or variants of the SLOS that require specific backfilling practices given different stope types [28]. In addition, these integrated methods are deterministic and are, thus, unable to account for uncertainty.

Deterministic LOM approaches lead to forecasts and production schedules that deviate from key production targets [29–32]. Dimitrakopoulos and Grieco [33] analyse how conventional methods for underground mine design are unable to capture the upside potential and/or downside risk of meeting forecasts, which are tied to the presence of uncertainty and variability of grades and material types. In order to account for geological uncertainty, stochastic simulations of the orebody are used as the main input to probabilistic and stochastic frameworks to generate a stope layout [33–36]. The stochastic optimisation of the underground mine production schedule, given the stope boundaries, is modelled for cut-and-fill, block-caving, and long-hole stoping mining methods [37–40]. Starting from Brickey [14] a resource-constrained project scheduling problem is proposed to address the selection of underground mine activities following their related precedence and duration. Nesbitt et al. [41] improve this previous work by incorporating uncertainty in the activities duration and grades, while Hill et al. [42] address computational techniques to reduce the problem size and algorithms to solve the optimisation. This activity scheduling approach, however, is suitable for shorter-term planning since the duration of each underground mining activity cannot be described on a yearly basis. Recently, a two-stage stochastic integer programming (SIP) [43] formulation that jointly optimises the stope design and production schedule for the SLOS mining method is developed by Furtado e Faria, et al. [44]. The mathematical model presented selects the stopes shapes and locations, and levels of extraction, using first-stage binary decision variables, in order to maximise NPV while considering effective vertical and horizontal development costs. The model also manages the geologic risk by minimising deviations from production targets modelled through second-stage decision variables. The related case study shows an improvement in the NPV when compared to the sequential framework, where a stope layout is given as input to the optimisation of the production schedule. Nevertheless, this method is tailored for a specific variant of SLOS that does not account explicitly for adjacency constraints when backfilling is used.

A new integrated stochastic optimisation of stope design and long-term production scheduling is proposed herein to go beyond the mining method specificities of previous approaches, in order to be applied to the SLOS mining method with backfilling. A variant of the SLOS is considered with assumptions and parameters derived from an existing operational mine. Horizontal extraction levels define the vertical boundaries of stopes; sublevel drifts and crosscuts are developed to enable longhole drilling. Primary, secondary, and tertiary stopes are aligned with extraction and backfilling procedures, and are combined with a bottom-up extraction approach to create precedence rules among stopes [45–48]. In addition, geometric parameters for shapes and sizes of stopes are defined for different mining zones according to geotechnical characteristics and requirements.

The proposed integrated optimisation is formulated as a two-stage SIP [43]. The method considers the selection of the stopes’ period of extraction and mining zone configuration that relates to the definition of stope shape parameters and types, as well as to sublevel positions. The maximisation of the NPV as per the objective function considers cumulative horizontal development costs and different haulage costs for different possible network systems. Geological uncertainty is modelled through a set of equiprobable geostatistical simulations of the mineral deposit [49–51]. The risk of not meeting production targets, which is inherent to geological uncertainty, is managed in the objective function through the incorporation of geological risk discounting (GRD) [52]. This discounting factor is applied to the objective function component minimising deviations from production targets, leading to the extraction of more uncertain materials in later periods on the LOM, when more information about the mineral deposit becomes available [52–58]. Physical and capacity constraints are included. Precedence rules of stopes defined based on geotechnical requirements enable the use of backfill after a stope is extracted.

The subsequent sections of this paper progress as follows. First, the main inputs to the integrated stochastic optimisation for the underground mine production scheduling and the correspondent mathematical programming formulation are outlined. Then, a case study at an operating underground mine is presented, including a comparison with a sequential stochastic framework. Conclusions and future work follow.

1.1. Method

A method for the integrated stochastic optimisation of stope design and mine production scheduling for the sublevel longhole open stoping (SLOS) mining method with backfilling is considered and presented below. The approach considers that a mining cycle, defined by all unit operations, such as drilling, blasting, loading, hauling, and backfilling, is completed for each stope on its mining period. A pattern of extraction, according to the type of stopes (e.g. primary, secondary, and tertiary), ensures that geotechnical constraints are met. Backfilling guarantees the stability of stope walls, removing the need for pillars. The sublevels are all aligned with the stopes’ lower bounds and can be used as extraction levels.

A set of geostatistical simulations $\mathit{s}\in \mathit{S}$ of the orebody describes the geological uncertainty. Initially, the orebody is represented in terms of blocks $\mathit{i}\in \mathit{I}$ that are, then, grouped into stopes $\mathit{j}\in \mathit{J}$. The bottom of the stopes that share the same sublevel $\mathit{l}\in \mathit{L}$ must be vertically aligned in order to define the location of developments, such as drifts $\mathit{d}\in D_l$. The requirement for the creation of possible cross-cuts $\mathit{c}\in C_l$ is that the stopes must also be laterally aligned throughout a mining direction. A set of possible configurations $\mathit{b}\in \mathit{B}$ of stopes that overlap each other is generated, allowing the stopes to have different allowable shapes, locations, and type options $\mathit{a}\in \mathit{A}$ that define the ordering of primary, secondary and tertiary stopes. The indices, sets, technical, economical, and geometrical parameters are listed in Tables 1-4, respectively. The decision variables of the proposed mathematical formulation are shown in Tables 5 and 6.

Table 1. List of indices.

Index	Definition
i	Block index
j	Stope index
b	Mining zone configuration index
d	Drift direction index
c	Crosscut index
l	Production level/sublevel index
a	Stope type option
k	Primary, secondary, and tertiary stope index
s	Index of a scenario quantifying the considered sources of uncertainty
h	Primary access system (shaft or ramp) index
$\mathit{\epsilon }$	Element (metal) index
t	Production period index

Table 2. List of sets.

Index	Definition
H	Set of primary access systems (shaft, ramp)
${\mathit{H}}_{\mathit{ramp}}$	Sub-set of primary access systems h which are ramps
${\mathit{H}}_{\mathit{shaft}}$	Sub-set of primary access systems h which are shafts
B	Set of all possible configurations b for the mining zone
${\mathit{B}}_{\mathit{h}}$	Set of possible mining zones b using the primary access system h
${\mathit{A}}_{\mathit{b}}$	Set of stope type option a for a mining zone b
${\mathit{L}}_{\mathit{b}}$	Set of sublevels l for mining zone configuration b
${\mathit{D}}_{\mathit{l}}$	Set of drift directions d in sublevel l
${\mathit{C}}_{\mathit{l}}$	Set of potential crosscuts c in sublevel l
${\mathit{J}}_{\mathit{b}}$	Set of potential stopes j in a mining zone configuration b
${\mathit{J}}_{\mathit{ba}}$	Set of potential stopes j in a mining zone configuration b, for stope type option a
${\mathit{J}}_{\mathit{bl}}$	Set of potential stopes j in a mining zone configuration b and using sublevel l
${\mathit{J}}_{\mathit{bdl}}$	Set of possible stopes in level l and drift direction d, and in mining zone configuration b
${\mathit{J}}_{\mathit{bcl}}$	Set of possible stopes in level l and cross-cut c and in mining zone configuration b
${\mathbf{\Phi }}_{\mathit{ja}}$	Set of predecessors of stope j in stope sequencing option a
${\mathit{I}}_{\mathit{b}}$	Set of blocks in a possible mining zone configuration b
${\mathit{I}}_{\mathit{jba}}$	Set of blocks i in stope j in mining zone configuration b, for stope type option a
S	Set of scenarios s describing the considered sources of uncertainty
T	Set of mining periods t
E	Set of elements $\mathit{\epsilon }$
K	Set of stope types (primary, secondary, tertiary)

Table 3. List of technical and economic parameters.

Index	Definition
${\mathit{\nu }}_{\mathit{i}}$	Volume of block i
${\mathit{\nu }}_{\mathit{jba}}$	Volume of stope j in mining zone configuration b for stope type option a ${\nu }_{\mathit{jba}}=\sum _{\mathit{i}\in {\mathit{I}}_{\mathit{jba}}}{\nu }_i$
${\mathit{w}}_{\mathit{is}}$	Tonnage of block i in scenario s
${\mathit{w}}_{\mathit{jbas}}$	Tonnage of stope j, in mining zone configuration b, for stope type option a and in geological scenario s $w_{\mathit{jbas}}=\sum _{\mathit{i}\in {\mathit{I}}_{\mathit{jb}}}{w}_{\mathit{is}}$
${\mathit{g}}_{\mathit{i\epsilon s}}$	Grade of element $\mathit{\epsilon }$ in block i, in scenario s
${\mathit{g}}_{\mathit{jba\epsilon s}}$	Average grade of element $\mathit{\epsilon }$ within stope j in mining zone b, for stope type option a and in scenario s $g_{\mathit{jba\epsilon s}}=\sum _{\mathit{i}\in {\mathit{I}}_{\mathit{jb}}}{w}_{\mathit{is}}{g}_{\mathit{i\epsilon s}}/w_{\mathit{jbas}}$
${\mathit{\pi }}_{\mathit{kjba}}^{\mathit{\hspace{0.17em}}}$	Extraction sequence type indicator (primary, secondary or tertiary) ${\pi }_{\mathit{kjba}}^{\mathit{\hspace{0.17em}}}$=1 if stope j is of type k in stope sequencing option a, in mining zone configuration b, and 0 otherwise
${\mathit{\rho }}_{\mathit{jba}}^{\mathit{\hspace{0.17em}}}$	Backfilling density ton/m^3 in stope j, in mining zone configuration b, in stope sequencing option a
${\mathit{R}}_{\mathit{\epsilon }}$	Processing recovery in percent of element $\mathit{\epsilon }$
${\mathit{P}}_{\mathit{\epsilon }}$	Metal selling price $/t of element $\mathit{\epsilon }$
${\mathit{v}}_{\mathit{jbas}}$	Economic value of stope j in mining zone b, stope sequencing option a and in scenario s
${\mathit{C}}_{\mathit{l}}^{\mathit{hor}}$	Unit horizontal development cost in sublevel l in $/km
$\mathit{C}{ }_{\mathit{k}}^{\mathit{\hspace{0.17em}}}\mathit{mine}$	Mining cost for type k stopes in $/t
${\mathit{C}}^{\mathit{proc}}$	Processing cost in $/t
${\mathit{C}}_{\mathit{hb}}^{\mathit{haul}}$	Haulage cost from mining zone b to primary access if $\mathit{h}\in H_{\mathit{ramp}}$ in $\mathrm{\$}/\left(\mathrm{tons}\cdot \mathrm{km}\right)\mathit{\hspace{0.17em}}$ and if $\in H_{\mathit{shaft}}$ in $/t
${\mathit{F}}_{\mathit{b}}$	Fixed cost for keeping the mining zone configuration b
${\mathit{f}}_{\mathit{t}}^{\mathit{EDR}}$	Economic discount factor for period t given an economic discount rate
${\mathit{f}}_{\mathit{t}}^{\mathit{GRD}}$	Geologic discount factor for period t given a geologic discount rate
${\mathit{U}}_{\mathit{ht}}^{\mathit{haul}}$	Hoisting capacity of primary access system h in period t (tons/year).
${\mathit{U}}_{\mathit{t}}^{\mathit{proc}}$	Processing capacity in period t (tons/year).
${\mathit{U}}_{\mathit{\epsilon t}}^{\mathit{\hspace{0.17em}}}$	Maximum target of element $\mathit{\epsilon }$ in period t (% element/year).
${\mathit{L}}_{\mathit{\epsilon t}}^{\mathit{\hspace{0.17em}}}$	Minimum target of element $\mathit{\epsilon }$ in period t (% element/year).
${\mathit{U}}_{\mathit{kt}}^{\mathit{bf}}$	Backfilling capacity for backfill type associated to k in period t (tons/year)
${\mathit{U}}_{\mathit{t}}^{\mathit{dev}}$	Development capacity in period t (m/year)
${\mathit{c}}_{\mathit{h}}^{\mathit{haul}}$	Penalty costs associated with the surplus deviations from the haulage tonnage capacity.
${\mathit{c}}^{\mathit{proc}}$	Penalty costs associated with the surplus deviations from the mill tonnage capacity.
${\mathit{c}}_{\mathit{\epsilon }}^{+}$	Penalty cost associated with surplus deviations from the maximum grade requirement for element $\mathit{\epsilon }$
${\mathit{c}}_{\mathit{\epsilon }}^{-}$	Penalty cost associated with shortage deviations from the minimum grade requirement for element $\mathit{\epsilon }$

Table 4. List of geometric parameters.

Index	Definition
${\mathit{\gamma }}_{\mathit{bj}}^{\mathit{x}},\mathit{\hspace{0.17em}}{\mathit{\gamma }}_{\mathit{bj}}^{\mathit{y}}$	Stope $\mathit{j}\in J_b$ sizes along direction x and y in the horizontal plane, in terms number of blocks for mining zone configuration $\mathit{b}\mathit{\hspace{0.17em}}\in \mathit{\hspace{0.17em}}\mathit{B}$.
${\mathit{\gamma }}_{\mathit{bj}}^{\mathit{z},\mathit{min}},\mathit{\hspace{0.17em}}{\mathit{\gamma }}_{\mathit{bj}}^{\mathit{z},\mathit{max}}$	Stope $\mathit{j}\in J_b$ minimum and maximum sizes along direction z (vertical plane) and in number of blocks for mining zone configuration $\mathit{b}\mathit{\hspace{0.17em}}\in \mathit{\hspace{0.17em}}\mathit{B}$.
${\mathit{\alpha }}_{\mathit{bl}}^{\mathit{z}}$	Sublevel spacing in z direction, for a mining zone configuration b and sublevel l
${\mathit{\delta }}_{\mathit{hbl}}^{\mathit{\hspace{0.17em}}}$	Length in km from the surface to the sublevel l of mining zone configuration b for primary access if $\mathit{h}\in H_{\mathit{ramp}}$ and 1 if $\mathit{h}\in H_{\mathit{shaft}}$
${\mathit{\delta }}_{\mathit{jdlb}}^{\mathit{drift}}$	Distance in a drift from a stope j to the h access point/ore pass along drift direction d in level l, in configuration b
${\mathit{\delta }}_{\mathit{jclb}}^{\mathit{crosscut}}$	Horizontal distance from stope j to the footwall along crosscut c in sublevel l under mining zone configuration b
${\mathit{\eta }}_{\mathit{b}}^{\mathit{sublevels}}$	Number of sublevels in mining zone configuration b
${\mathit{\eta }}_{\mathit{b}}^{\mathit{stopes}}$	Number of stopes in mining zone configuration b

Table 5. Binary decision variables.

Index	Definition
${\mathit{z}}_{\mathit{ba}}$	Mining zone selection decision variable, equal to 1 if mining zone configuration b using stope type option a is selected, and 0 otherwise
${\mathit{y}}_{\mathit{jbat}}$	Stope selection decision variable, equal to 1 if stope j in mining zone configuration b is selected, for stope sequencing option a, in period t, and 0 otherwise

Table 6. Continuous decision variables.

Index	Definition
${\mathit{\psi }}_{\mathit{dlbt}}^{\mathbf{drift}}$	Drift’s development distance in sublevel l, in mining zone configuration b, along with drift direction $\mathit{d}\mathit{\hspace{0.17em}}\in \mathit{\hspace{0.17em}}{D}_l$ in period t
${\mathit{\psi }}_{\mathit{dlbt}}^{\mathbf{drif}{\mathbf{t}}^{\ast }}$	Effective drift’s development distance in sublevel l, in mining zone configuration b, along with drift direction $\mathit{d}\mathit{\hspace{0.17em}}\in \mathit{\hspace{0.17em}}{D}_l$ in period t
${\mathit{\psi }}_{\mathit{clbt}}^{\mathit{crosscut}}$	Crosscut’s development distance in sublevel l in mining zone configuration b, for crosscut $\mathrm{c}\in \mathit{\hspace{0.17em}}{C}_l$ in period t
${\mathit{\psi }}_{\mathit{clbt}}^{\mathit{crosscu}{\mathit{t}}^{\ast }}$	Effective crosscut’s development distance in sublevel l in mining zone configuration b, for crosscut $\mathrm{c}\in \mathit{\hspace{0.17em}}{C}_l$ in period t
${\mathit{d}}_{\mathit{hts}}^{\mathit{haul}}$	Surplus deviation of total hoisting/haulage capacity for the system h, in period t and scenario s
${\mathit{d}}_{\mathit{ts}}^{\mathit{proc}}$	Surplus deviation of total processing capacity, in period t, and scenario s
${{\mathit{d}}_{\mathit{\epsilon }\mathbf{ts}}}^{+}$	Surplus deviation from the maximum grade target of element $\mathit{\epsilon }$, in period t, and scenario s
${{\mathit{d}}_{\mathit{\epsilon }\mathbf{ts}}}^{-}$	Shortage deviation from the maximum grade target of element $\mathit{\epsilon }$, in period t, and scenario s

1.2. Input data processing

Figure 1 shows the three main steps to generate potential stopes from an orebody model by considering a given mining zone; that is, a portion of the deposit that has the same geomechanical characteristics. It is assumed that, for a given mining zone, stopes share the same geometrical parameters that define their shapes. The first step consists of dividing the orebody model into different sublevels and mining fronts. This process is done repeatedly for all allowable sublevel spacing and stope widths, producing a set of possible mining zone configurations $\mathit{b}\in \mathit{B}$. Subsequently, stope type options $A_b$ are mapped for each configuration $\mathit{b}\in \mathit{B}$. Each option $\mathit{a}\in A_b$ assigns a different type (primary, secondary, or tertiary) to sets of stopes $J_{\mathit{ba}}$. These stopes undergo an additional stage of processing that determines the stope height that yields a probabilistic highest profitability. Steps 1 and 2 are invariant to the geostatistical simulations of the orebody, while step 3 considers all orebody simulations jointly. The output of this process defines possible geometries and positions of the stopes and sublevels, which, combined with the grade and material type simulations of the orebody, are finally used as inputs to the proposed two-stage SIP.

Figure 1. Steps of the stope design and scheduling optimization.

1.2.1. Mapping configurations of shapes

First, starting from a block $\mathit{i}\in I_b$ in a mining zone, blocks are grouped according to stope dimensions, ${\mathit{\gamma }}_{\mathit{bj}}^x,{\mathit{\gamma }}_{\mathit{bj}}^y\mathit{\hspace{0.17em}}$ and the sublevel spacing ${\mathit{\alpha }}_{\mathit{bl}}^z$ defined for a given mining zone configuration $\mathit{b}\mathit{\hspace{0.17em}}\in \mathit{\hspace{0.17em}}\mathit{B}$, possible stopes $J_b$ and possible sublevels $\mathit{l}\in L_b$ within that mining zone configuration. This process is repeated for all possible configurations. Figure 2 shows two different mining zone configurations b, b`∈ B where configuration b assumes three equally spaced sublevels and configuration b` has five sublevels. It is worth noting that two different indexed mining zone configurations can have an identical structure in terms of number of mining fronts and sublevels, as they differ in the assigned haulage system option h in terms of type (i.e. shaft of ramp) and position.

Figure 2. Two mining zone configurations $\mathit{b}\in \mathit{B}$ generated in the mapping of shapes step.

It is important to note that the stope dimensions in the horizontal plane can be variable within the mining zone configuration and are directly related to the final possible stope shapes. These dimensions follow mainly geotechnical requirements related to minimum and maximum dimensions of stopes. The option of having variable shapes considering the horizontal plane is critical when there are different mining or backfilling costs for different stope types, as it enables the management of these costs according to the stope volume and type. On the other hand, for the vertical direction, the distance between sublevels, which are also potential extraction levels, is considered. The spacing ${\mathit{\alpha }}_{\mathit{bl}}^z$ only defines the sublevels l bounding the stopes and does not determine the final height of the potential stopes. A specific step to determine the actual stope height, which considers not only geometric parameters, but the economic profitability of the given mineable volume is presented in Section 2.1.3. In addition, the total number of blocks in a given direction might not be divisible by the stope dimension in that direction. Thus, configurations in which the grouping of blocks into stopes starts from different blocks in the mining zone can be generated.

As seen in Figure 3, cross-cuts c are developed parallel to a defined mining direction, in order to meet ventilation and backfilling requirements. Drifts d are developed perpendicularly to cross-cuts. The approximated dimensions of crosscuts (${\mathit{\delta }}_{\mathit{jclb}}^{\mathit{crosscut}}$) and drifts (${\mathit{\delta }}_{\mathit{jdlb}}^{\mathit{drift}-\mathit{h}}$) are calculated considering the position of the sublevels l. The drift length is always associated with the access point of a haulage system h and its respective mining zone configurations $\mathit{\hspace{0.17em}}\mathit{b}\in B_h$.

Figure 3. Drifts and cross-cuts distances for a haulage system h, a potential sublevel l (in a plan view) for two potential stopes j and j`.

1.2.2. Mapping stope type options

Once the possible mining zone configurations are defined, for each configuration $\mathit{b}$, stope type options $\mathit{a}\in A_b$ are generated. A pattern of extraction is defined following geotechnical constraints. Figure 4 shows a stope-to-stope dependency that is adapted from different patterns presented by Villaescusa [28], in which the numbers define the predecessors and the colouring defines the stope type. A stope with number 1 has no predecessors, a stope tagged with number 2 has stopes with number 1 as a predecessor, and those with number 3 have all the previous ones as predecessors; this rule follows for all other stopes. Therefore, the adjacencies are mapped, such that stopes with lower numbers in Figure 4 are predecessors ($\mathit{\phi }\mathit{\hspace{0.17em}}\in \mathit{\hspace{0.17em}}{\mathrm{\Phi }}_{\mathit{ja}}$) of stopes with higher numbers ($\mathit{j}\in J_b$). It is important to note that this pattern does not predefine the extraction sequence and it is used only to define adjacencies for rock mass stability purposes.

Figure 4. Two consecutive cross-sections of the stope pattern of extraction.

The generated stopes, that are outputs from the previous step, are flagged with an indicator ${\pi }_{\mathit{kjba}}^{\mathit{\hspace{0.17em}}}$ type, where $\mathit{k}\in \mathit{K}$ defines the type of a stope j (i.e. primary, secondary, or tertiary), in a mining zone configuration b and considering a type option a that follows a dependency pattern. Thus, each type option a shows a different combination of types k and stopes j. Figure 5 represents the process of mapping the stope type options. It is observed that different combinations of mining zone configurations and stope type options are generated.

Figure 5. Stope type option for different mining zone configurations.

1.2.3. Searching variable stope heights

In the previous steps, the stopes generated come from the grouping of blocks by considering that they fully occupy the space between sublevels. In order to manage dilution, the profitability of having stopes with heights that are smaller than the distance between sublevels $\left({\mathit{\alpha }}_{\mathit{bl}}^z\right)$ and greater than the minimum stope height (${\mathit{\gamma }}_{\mathit{bj}}^{\mathit{z},\mathit{min}}$) is evaluated. This means that the alignment of the bottom of the stopes according to the sublevel that defines their lower bound is kept. However, the height might be variable within the same sublevel, and there may be waste material above certain stopes, being left as pillars. For this variant of SLOS, a sublevel can serve as an extraction level. Therefore, it is assumed that stopes above the waste portion can still be mined subsequently.

For a given mining zone configuration b, a stope type option a and a stope j, the economic value of each stope is calculated, as shown in Equation. 1, for all its possible vertical dimensions and for each simulated orebody scenario s of grades of different elements$\mathit{\epsilon }\in \mathit{E}$ and stope tonnages $w_{\mathit{jbas}}$. For each vertical dimension, given a probability threshold, the economic value that defines this quantile is calculated ($v_{\mathit{height}}^{\mathrm{\%}}$). For instance, if the P90 defines the probability threshold, there will be an economic value associated with a scenario s where 90% of all scenarios will either be equal to or will not exceed this value. Then, the economic value associated with different vertical dimensions ($v_{\mathit{height}}^{\mathrm{\%}})$ are compared, the dimension associated with the highest economic value is retained (${\mathit{\gamma }}_{\mathit{jba}}^{\mathit{z},\mathit{best}})$, defining the final possible stope shape, for that configuration b and the associated stope type option a. Figure 6 illustrates how the best stope height is chosen. It should be noted that the stope type option a definition described on the previous section (i.e. step two) does not consider dimensional parameters. However, the mining cost associated to the different stope types have an impact while searching for the best stope height. Consequently, on this third step, the stope height, volume, weight, and grade will be correlated to its possible type and indexed accordingly.

(1) $\begin{array}{c}{v}_{\mathit{jba},\mathit{s}}=w_{\mathit{jbas}}\left(\underset{\mathit{\epsilon }\in \mathit{E}}{{{\displaystyle \sum }}^{\text{}}}{g}_{\mathit{jba\epsilon s}}{R}_{\epsilon }{P}_{\epsilon }-\left(C^{\mathit{proc}}+\underset{\mathit{k}\in \mathit{K}}{{{\displaystyle \sum }}^{\text{}}}{\pi }_{\mathit{kjba}}^{}\mathit{C} k^{\mathit{mine}}\right)\right),\forall \text{b}\in \mathit{B},\mathit{a}\in A_b,\mathit{j}\in J_{\mathit{ba}},\mathit{s}\in \mathit{S}\end{array}$(1)

Figure 6. Example of how to define the stope height considering ${\mathit{\gamma }}_{\mathit{bj}}^{\mathit{z},\mathit{min}}=2$ and ${\mathit{\gamma }}_{\mathit{bj}}^{\mathit{z},\mathit{max}}=4$. The $ represents the economic value of a stope considering the $v_{\mathit{height}}^{\mathit{P}90}$. The circled stope is the best stope height ${\mathit{\gamma }}_{\mathit{jba}}^{\mathit{z},\mathit{best}}$selected.

1.3. Stochastic integer programming formulation

This section presents a mathematical programming formulation model, which is developed to optimise the underground mine production schedule and stope boundaries jointly, while considering uncertainty in material supply. Two binary decision variables are shown in Table 5. The mining zone configuration decision variables $z_{\mathit{ba}}\in \left\{0,1\right\}$ control which mining zone configuration $\mathit{b}\in \mathit{\hspace{0.17em}}\mathit{B}$ and respective stope type option is selected $\mathit{a}\in A_b$. These decision variables impact directly the selections of stopes shapes and types. A mining zone configuration is always associated with a single haulage system $\mathit{h}\in \mathit{H}$. This means that identical mining zone configuration options ($\mathit{b}$ and b`) in terms of stope shapes and sublevels can exist, but they will be associated with different available haulage systems ($\mathit{b}\in B_h,$ and b`∈ B_h`). It is assumed that vertical accesses compatible with the haulage systems are already developed, which enables the variable $z_{\mathit{ba}}$ to be time-independent. The stope selection decision variables $y_{\mathit{jba},}\in \left\{0,1\right\}$ determine if a stope $\mathit{j}\in J_b$ in a mining zone configuration $\mathit{b}\in \mathit{\hspace{0.17em}}\mathit{B}$, using stope type option $\mathit{a}\in A_b$ is mined in period t. It is assumed that all the unitary operations; that is, development of secondary accesses, drilling, blasting, hauling, and backfilling, are ready or done at the period scheduled for mining a given stope.

Two continuous decision variables ${\mathit{\psi }}_{\mathit{dlbt}}^{\mathit{drift}}$ and ${\mathit{\psi }}_{\mathit{clbt}}^{\mathit{crosscut}}$ correspond to the developed distance of a drift d or a cross-cut c, for a sublevel l, in a mining zone configuration b, in period t. To take into account the available structures developed in previous years, effective development distance ${\mathit{\psi }}_{\mathit{dlbt}}^{\mathit{drif}{\mathit{t}}^{\ast }}$and ${\mathit{\psi }}_{\mathit{clbt}}^{\mathit{crosscu}{\mathit{t}}^{\ast }}$ are used in practice, they correspond to cumulative horizontal development distances. The remaining decision variables presented in Table 6 refer to surplus deviations from haulage capacities for different haulage systems h$\left(d_{\mathit{hts}}^{\mathit{haul}}\right)$ and processing capacity ($d_{\mathit{ts}}^{\mathit{proc}}$), and deviations from lower and upper bounds for different elements $\mathit{\epsilon }\in \mathrm{E}$ requirements, ${{\mathit{d}}_{\mathit{\epsilon }\mathrm{ts}}}^{-}$ and ${{\mathit{d}}_{\mathit{\epsilon }\mathrm{ts}}}^{+}$ respectively.

1.3.1. Objective function

This section introduces the objective function of the proposed SIP and describes its main components as follows.

$max\underset{\mathit{Part}\mathit{I}:\mathit{Discounted}\mathit{\hspace{0.17em}}\mathit{Revenue}\mathit{\hspace{0.17em}}\mathit{from}\mathit{\hspace{0.17em}}\mathit{Scheduled}\mathit{\hspace{0.17em}}\mathit{Stopes}}{\underbrace{\frac{1}{\left|\mathit{S}\right|}\sum _{\mathrm{s}\in \mathrm{S}}\sum _{\mathit{t}\mathit{\hspace{0.17em}}\in \mathrm{T}}\sum _{\mathit{b}\mathit{\hspace{0.17em}}\in \mathrm{B}}\sum _{\mathrm{a}\in \mathit{\hspace{0.17em}}{A}_b}\sum _{\mathrm{j}\in \mathit{\hspace{0.17em}}{\mathit{J}}_{\mathit{ba}}}{f}_t^{\mathit{EDR}}{v}_{\mathit{jbas}}\mathit{\hspace{0.17em}}{y}_{\mathit{jbat}}}}$

$\underset{\mathit{Part}\mathit{II}:\mathit{Haulage}\mathit{\hspace{0.17em}}\mathit{Costs}}{\underbrace{-\frac{1}{\left|\mathit{S}\right|}\sum _{\mathrm{s}\in \mathrm{S}}\sum _{\mathit{t}\mathit{\hspace{0.17em}}\in \mathit{\hspace{0.17em}}T}{f}_t^{\mathit{EDR}}\left(\sum _{\mathit{h}\in \mathit{H}}\sum _{\mathit{b}\mathit{\hspace{0.17em}}\in \mathit{\hspace{0.17em}}{B}_h}\sum _{l\in L_b}\sum _{\mathit{a}\mathit{\hspace{0.17em}}\in \mathit{\hspace{0.17em}}{A}_b}\sum _{\mathit{j}\mathit{\hspace{0.17em}}\in \mathit{\hspace{0.17em}}{\mathit{J}}_{\mathit{bl}}}\mathit{\hspace{0.17em}}{y}_{\mathit{jbat}}\mathit{\hspace{0.17em}}{w}_{\mathit{jbas}}\mathit{\hspace{0.17em}}{\delta }_{\mathit{hbl}}^{\mathit{\hspace{0.17em}}}{C}_{\mathit{hb}}^{\mathit{haul}}\right)}}$

$-\underset{\mathit{Part}\mathit{III}:\mathit{Development}\mathit{\hspace{0.17em}}\mathit{Costs}}{\underbrace{\sum _{\mathrm{t}\in \mathrm{T}}\sum _{\mathit{b}\mathit{\hspace{0.17em}}\in \mathrm{B}}\sum _{l\in L_b}{f}_t^{\mathit{EDR}}\mathit{\hspace{0.17em}}{C}_l^{\mathit{\hspace{0.17em}}hor}\left(\sum _{\mathit{d}\mathit{\hspace{0.17em}}\in \mathit{\hspace{0.17em}}{D}_l}{\mathit{\psi }}_{\mathit{dlbt}}^{\mathit{drif}{\mathit{t}}^{\ast }}+\sum _{\mathit{c}\mathit{\hspace{0.17em}}\in \mathit{\hspace{0.17em}}{C}_l}{\mathit{\psi }}_{\mathit{clbt}}^{\mathit{crosscu}{\mathit{t}}^{\ast }}\right)}}$

$-\underset{\mathit{Part}\mathit{IV}:\mathit{Fixed}\mathit{\hspace{0.17em}}\mathit{Cost}}{\underbrace{\sum _{\mathit{b}\mathit{\hspace{0.17em}}\in \mathrm{B}}\sum _{\mathit{a}\mathit{\hspace{0.17em}}\in \mathit{\hspace{0.17em}}{A}_b}{f}_t^{\mathit{EDR}}{F}_b\mathit{\hspace{0.17em}}{z}_{\mathit{ba}}}}$

(2) $-\underset{\mathit{Part}\mathit{V}:\mathit{Geological}\mathit{\hspace{0.17em}}\mathit{Risk}\mathit{\hspace{0.17em}}\mathit{Management}}{\underbrace{\frac{1}{\left|\mathit{S}\right|}\sum _{\mathit{s}\in \mathrm{S}}\sum _{\mathit{t}\in \mathit{T}}{f}_t^{\mathit{GRD}}\left(c_h^{\mathit{haul}}{d}_{\mathit{hts}}^{\mathit{haul}}+c^{\mathit{proc}}{d}_{\mathit{ts}}^{\mathit{proc}}+\sum _{\mathit{\epsilon }\in \mathit{E}}{c}_{\epsilon }^{+}\mathit{\hspace{0.17em}}{{\mathrm{d}}_{\mathit{\epsilon }\mathrm{ts}}}^{+}+c_{\epsilon }^{-}{{\mathit{d}}_{\mathit{\epsilon ts}}}^{-}\right)}}$(2)

The objective function shown in has five parts. Part I aims to maximise the discounted revenue from scheduled stopes. Part II of the objective function minimises the haulage cost of mined stopes. When a mining zone configuration is associated with a ramp system, the haulage cost $C_{\mathit{hb}}^{\mathit{haul}}$ is expressed in $\mathrm{\$}/\left(\mathrm{t}\cdot \mathrm{km}\right)$. Thus, the length ${\delta }_{\mathit{h},\mathit{b},\mathit{l}}^{\mathit{\hspace{0.17em}}}$ of the ramp until a sublevel l must be taken into consideration, while for the shaft these values will always equal one. Part III of Equation.2 minimises the horizontal development costs. Part IV minimises particularly fixed costs associated with different mining zone configurations and patterns of extraction. Part V manages the geological risk, by minimising deviations from production targets, related to mining and processing capacities and grade blending requirements. For that purpose, penalty costs $c_h^{\mathit{haul}},$ $c^{\mathit{proc}},c_{\epsilon }^{+}\mathrm{and}{c}_{\epsilon }^{-}$are applied correspondently to the production requirements and targets, as they are discounted by a geological risk discounting factor $f_t^{\mathit{GRD}}$. Therefore, riskier stopes will be scheduled in later periods when more information regarding grades and material quality is available [52,53].

1.3.2. Constraints

The objective function is subjected to the following constraints.

(3) $\begin{array}{c}\sum _{\mathit{b}\in \mathit{B}}\sum _{\mathrm{a}\in A_b}{z}_{\mathit{ba}}\le 1\end{array}$(3)

Equation 3(3) $\begin{array}{c}\sum _{\mathit{b}\in \mathit{B}}\sum _{\mathrm{a}\in A_b}{z}_{\mathit{ba}}\le 1\end{array}$(3) imposes that a single mining zone configuration with a correspondent stope option is selected.

(4) $\begin{array}{c}\sum _{\mathit{t}\in \mathrm{T}}{y}_{\mathit{jbat}}\le z_{\mathit{ba}},\mathrm{\forall b}\in \mathit{B},\mathit{a}\in A_b,\mathit{j}\in J_{\mathit{ba}}\end{array}$(4)

Liking constraints (Equation 4(4) $\begin{array}{c}\sum _{\mathit{t}\in \mathrm{T}}{y}_{\mathit{jbat}}\le z_{\mathit{ba}},\mathrm{\forall b}\in \mathit{B},\mathit{a}\in A_b,\mathit{j}\in J_{\mathit{ba}}\end{array}$(4) ) ensure that a scheduled stope belongs to the chosen mining zone configuration.

(5) $\begin{array}{c}\sum _{\mathit{t}\in \mathit{T}}{y}_{\mathit{jbat}}\le 1,\mathrm{\forall b}\in \mathit{B},\mathit{a}\in A_b,\mathit{j}\in J_{\mathit{ba}}\end{array}$(5)

Equation 5(5) $\begin{array}{c}\sum _{\mathit{t}\in \mathit{T}}{y}_{\mathit{jbat}}\le 1,\mathrm{\forall b}\in \mathit{B},\mathit{a}\in A_b,\mathit{j}\in J_{\mathit{ba}}\end{array}$(5) defines reserve constraints, where the stope can be mined only once.

(6) $\begin{array}{c}\left|{\mathrm{\Phi }}_{\mathit{ja}}\right|y_{\mathit{jbat}}\le \left|{\mathrm{\Phi }}_{\mathit{ja}}\right|-\sum _{a\in A_b}\sum _{\mathit{\phi }\in {\mathrm{\Phi }}_{\mathit{ja}}}\sum _{t^{\prime }=\mathit{t}+1}^{\left|\mathit{T}\right|}{y}_{\phi ba,t^{\prime }},\mathrm{\forall j}\in J_b,\mathit{b}\in \mathit{B},\mathit{t}\in \mathit{T}\end{array}$(6)

Equation 6(6) $\begin{array}{c}\left|{\mathrm{\Phi }}_{\mathit{ja}}\right|y_{\mathit{jbat}}\le \left|{\mathrm{\Phi }}_{\mathit{ja}}\right|-\sum _{a\in A_b}\sum _{\mathit{\phi }\in {\mathrm{\Phi }}_{\mathit{ja}}}\sum _{t^{\prime }=\mathit{t}+1}^{\left|\mathit{T}\right|}{y}_{\phi ba,t^{\prime }},\mathrm{\forall j}\in J_b,\mathit{b}\in \mathit{B},\mathit{t}\in \mathit{T}\end{array}$(6) ensures that adjacency constraints are met. The sublevel open stope mining method allows that stopes can be mined if their predecessor are mined before or are left unmined. Differently from open-pit mine scheduling problems, the set of predecessors ${\mathrm{\Phi }}_{\mathit{ja}}$ of a stope j must contain all direct and indirect predecessors. Thus, due to this constraint, the complexity of the model is highly impacted by the size of the orebody.

(7) $\begin{array}{c}{\mathit{\psi }}_{\mathit{dlbt}}^{\mathit{drift}}\ge \underset{\mathit{j}\in {\mathit{J}}_{\mathit{bdl}}}{max}\left({\mathit{\delta }}_{\mathit{jdlb}}^{\mathit{drift}}\sum _{\mathit{a}\in \mathit{A}}{y}_{\mathit{jbat}}\right),\mathrm{\forall b}\in \mathit{B},\mathit{l}\in L_b,\mathit{d}\in D_l,\mathit{t}\in \mathit{T}\end{array}$(7)

(8) $\begin{array}{c}{\mathit{\psi }}_{\mathit{dlb}1}^{\mathit{drif}{\mathit{t}}^{\ast }}={\mathit{\psi }}_{\mathit{dlb}1}^{\mathit{drift}},\mathrm{\forall b}\in \mathit{B},\mathit{l}\in L_b,\mathit{d}\in D_l,\mathit{t}=1\end{array}$(8)

(9) $\begin{array}{c}{\mathit{\psi }}_{\mathit{dlbt}}^{\mathit{drif}{\mathit{t}}^{\ast }}\ge \left\{\left[{\mathit{\psi }}_{\mathit{dlbt}}^{\mathit{drift}}-\sum _{t^{\prime }=1}^{\mathit{t}-1}{\mathit{\psi }}_{\mathit{dlbt}\prime }^{\mathit{drif}{\mathit{t}}^{\ast }}\right]\right\},\mathrm{\forall b}\in \mathit{B},\mathit{l}\in L_b,\mathit{d}\in D_l,\mathit{t}>1\end{array}$(9)

Equations 7(7) $\begin{array}{c}{\mathit{\psi }}_{\mathit{dlbt}}^{\mathit{drift}}\ge \underset{\mathit{j}\in {\mathit{J}}_{\mathit{bdl}}}{max}\left({\mathit{\delta }}_{\mathit{jdlb}}^{\mathit{drift}}\sum _{\mathit{a}\in \mathit{A}}{y}_{\mathit{jbat}}\right),\mathrm{\forall b}\in \mathit{B},\mathit{l}\in L_b,\mathit{d}\in D_l,\mathit{t}\in \mathit{T}\end{array}$(7) -9(9) $\begin{array}{c}{\mathit{\psi }}_{\mathit{dlbt}}^{\mathit{drif}{\mathit{t}}^{\ast }}\ge \left\{\left[{\mathit{\psi }}_{\mathit{dlbt}}^{\mathit{drift}}-\sum _{t^{\prime }=1}^{\mathit{t}-1}{\mathit{\psi }}_{\mathit{dlbt}\prime }^{\mathit{drif}{\mathit{t}}^{\ast }}\right]\right\},\mathrm{\forall b}\in \mathit{B},\mathit{l}\in L_b,\mathit{d}\in D_l,\mathit{t}>1\end{array}$(9) show how the drift development costs are calculated. As shown in Figure 7 the drift development distance ${\mathit{\psi }}_{\mathit{dlbt}}^{\mathit{drift}}$ corresponds to the distance from the furthest stope mined in a year t to the access point in the sublevel l. This value is used to calculate the effective development distance ${\mathit{\psi }}_{\mathit{dlbt}}^{\mathit{drif}{\mathit{t}}^{\ast }}$, that considers only the remaining length to be developed in a given period t considering the developments done in the previous years.

Figure 7. Example of drift development costs ${\mathit{\psi }}_{\mathit{dlbt}}^{\mathit{drift}}$ and effective development costs ${\mathit{\psi }}_{\mathit{dlbt}}^{\mathit{drif}{\mathit{t}}^{\ast }}$for the extraction of three stopes in three periods.

(10) $\begin{array}{c}{\mathit{\psi }}_{\mathit{clbt}}^{\mathit{crosscut}}\ge \underset{\mathit{j}\in {\mathit{J}}_{\mathit{bcl}}}{max}\left({\mathit{\delta }}_{\mathit{jclb}}^{\mathit{crosscut}}\sum _{\mathit{a}\in \mathit{A}}{y}_{\mathit{jbat}}\right),\mathrm{\forall b}\in \mathit{B},\mathit{l}\in L_b,\mathit{c}\in C_l,\mathit{t}\in \mathit{T}\end{array}$(10)

(11) $\begin{array}{c}{\mathit{\psi }}_{\mathit{clb}1}^{\mathit{crosscu}{\mathit{t}}^{\ast }}={\mathit{\psi }}_{\mathit{clb}1}^{\mathit{crosscut}},\mathrm{\forall b}\in \mathit{B},\mathit{l}\in L_b,\mathit{c}\in C_l,\mathit{t}=1\end{array}$(11)

(12) $\begin{array}{c}{\mathit{\psi }}_{\mathit{clbt}}^{\mathit{crosscu}{\mathit{t}}^{\ast }}\ge \left\{\left[{\mathit{\psi }}_{\mathit{clbt}}^{\mathit{crosscut}}-\sum _{t^{\prime }=1}^{\mathit{t}-1}{\mathit{\psi }}_{\mathit{clbt}\prime }^{\mathit{crosscu}{\mathit{t}}^{\ast }}\right]\right\},\mathrm{\forall b}\in \mathit{B},\mathit{l}\in L_b,\mathit{c}\in C_l,\mathit{t}>1\end{array}$(12)

Equations 10(10) $\begin{array}{c}{\mathit{\psi }}_{\mathit{clbt}}^{\mathit{crosscut}}\ge \underset{\mathit{j}\in {\mathit{J}}_{\mathit{bcl}}}{max}\left({\mathit{\delta }}_{\mathit{jclb}}^{\mathit{crosscut}}\sum _{\mathit{a}\in \mathit{A}}{y}_{\mathit{jbat}}\right),\mathrm{\forall b}\in \mathit{B},\mathit{l}\in L_b,\mathit{c}\in C_l,\mathit{t}\in \mathit{T}\end{array}$(10) ,11(11) $\begin{array}{c}{\mathit{\psi }}_{\mathit{clb}1}^{\mathit{crosscu}{\mathit{t}}^{\ast }}={\mathit{\psi }}_{\mathit{clb}1}^{\mathit{crosscut}},\mathrm{\forall b}\in \mathit{B},\mathit{l}\in L_b,\mathit{c}\in C_l,\mathit{t}=1\end{array}$(11) show how the crosscut development distances are calculated. The same method is used for the crosscuts when compared to the drifts. However, by considering a mining direction that should be strictly followed, the process of deferring the crosscut development costs will also define how far from the access is worth mining.

(13) $\begin{array}{c}\sum _{b\in B_h}\sum _{a\in A_b}\sum _{\mathit{j}\in {\mathit{J}}_{\mathit{ba}}}\left(y_{\mathit{jbat}}{w}_{\mathit{jbas}}\right)-{{\mathit{d}}_{\mathit{ts}}}^{\mathit{h}}\le U_t^h,\mathrm{\forall h}\in \mathit{H},\mathit{t}\in \mathit{T},\mathit{s}\in \mathit{S}\end{array}$(13)

(14) $\begin{array}{c}\sum _{\mathit{b}\in \mathit{B}}\sum _{a\in A_b}\sum _{\mathit{j}\in {\mathit{J}}_{\mathit{ba}}}\left(y_{\mathit{jbat}}{w}_{\mathit{jbas}}\right)-{{\mathit{d}}_{\mathit{ts}}}^{\mathit{p}}\le U_t^p,\mathrm{\forall t}\in \mathit{T},\mathit{s}\in \mathit{S}\end{array}$(14)

Equation 13(13) $\begin{array}{c}\sum _{b\in B_h}\sum _{a\in A_b}\sum _{\mathit{j}\in {\mathit{J}}_{\mathit{ba}}}\left(y_{\mathit{jbat}}{w}_{\mathit{jbas}}\right)-{{\mathit{d}}_{\mathit{ts}}}^{\mathit{h}}\le U_t^h,\mathrm{\forall h}\in \mathit{H},\mathit{t}\in \mathit{T},\mathit{s}\in \mathit{S}\end{array}$(13) defines an upper bound for the extracted material considering the hoisting or haulage capacity of the haulage system associated to the mining zone configuration. Also, a constraint that limits the production of ore mined is considered (Equation 14(14) $\begin{array}{c}\sum _{\mathit{b}\in \mathit{B}}\sum _{a\in A_b}\sum _{\mathit{j}\in {\mathit{J}}_{\mathit{ba}}}\left(y_{\mathit{jbat}}{w}_{\mathit{jbas}}\right)-{{\mathit{d}}_{\mathit{ts}}}^{\mathit{p}}\le U_t^p,\mathrm{\forall t}\in \mathit{T},\mathit{s}\in \mathit{S}\end{array}$(14) ). These two constraints are modelled as soft constraints by allowing deviations from the defined boundaries considering a general case where block tonnages, consequently, stope tonnages, are variable for different geological scenarios $\mathit{s}\in \mathit{S}$.

(15) $\begin{array}{c}\sum _{\mathit{b}\in \mathit{B}}\sum _{a\in A_b}\sum _{\mathit{j}\in {\mathit{J}}_{\mathit{ba}}}\left(y_{\mathit{j},\mathit{b},\mathit{a},\mathit{t}}{\rho }_{\mathit{jba}}^{\mathit{\hspace{0.17em}}}{\nu }_{\mathit{jba}}\right)\le U_{\mathit{kt}}^{\mathit{bf}},\mathrm{\forall k}\in \mathit{K},\mathit{t}\in \mathit{T}\end{array}$(15)

Backfilling capacity constraints (Equation 15(15) $\begin{array}{c}\sum _{\mathit{b}\in \mathit{B}}\sum _{a\in A_b}\sum _{\mathit{j}\in {\mathit{J}}_{\mathit{ba}}}\left(y_{\mathit{j},\mathit{b},\mathit{a},\mathit{t}}{\rho }_{\mathit{jba}}^{\mathit{\hspace{0.17em}}}{\nu }_{\mathit{jba}}\right)\le U_{\mathit{kt}}^{\mathit{bf}},\mathrm{\forall k}\in \mathit{K},\mathit{t}\in \mathit{T}\end{array}$(15) ) are considered in a way that different upper bounds $U_{\mathit{kt}}^{\mathit{bf}}$ can be chosen for different types of backfilling$\mathit{k}\in \mathit{K}$. In this case, deviations are not allowed since no uncertainty is associated with the stope volume ${\nu }_{\mathit{jb}}$.

(16) $\begin{array}{c}\sum _{\mathit{b}\in \mathit{B}}\left(\sum _{l\in L_b}\sum _{d\in D_l}{\mathit{\psi }}_{\mathit{dlbt}}^{\mathit{drif}{\mathit{t}}^{\ast }}+\sum _{l\in L_b}\sum _{C\in C_l}{\mathit{\psi }}_{\mathit{clbt}}^{\mathit{crosscu}{\mathit{t}}^{\ast }}\right)\le U_t^{\mathit{develop}},\mathrm{\forall t}\in \mathit{T}\end{array}$(16)

Horizontal development capacities are added in terms of the maximum length of the total drifts and crosscuts (Equation 16(16) $\begin{array}{c}\sum _{\mathit{b}\in \mathit{B}}\left(\sum _{l\in L_b}\sum _{d\in D_l}{\mathit{\psi }}_{\mathit{dlbt}}^{\mathit{drif}{\mathit{t}}^{\ast }}+\sum _{l\in L_b}\sum _{C\in C_l}{\mathit{\psi }}_{\mathit{clbt}}^{\mathit{crosscu}{\mathit{t}}^{\ast }}\right)\le U_t^{\mathit{develop}},\mathrm{\forall t}\in \mathit{T}\end{array}$(16) ).

(17) $\begin{array}{c}\sum _{\mathit{b}\in \mathit{B}}\sum _{\mathrm{a}\in A_b}\sum _{j\in J_b}\left(y_{\mathit{jbat}}{w}_{\mathit{jbsa}}\right)\left(g_{\mathit{jba\epsilon s}}-U_{\mathit{\epsilon t}}^{\mathit{\hspace{0.17em}}}\right)-{{\mathit{d}}_{\mathit{\epsilon }\mathrm{ts}}}^{+}\le 0,\mathrm{\forall }\mathit{\epsilon }\in \mathit{E},\mathit{t}\in \mathit{T},\mathit{s}\in \mathit{S}\end{array}$(17)

(18) $\begin{array}{c}\sum _{\mathit{b}\in \mathit{B}}\sum _{a\in A_b}\sum _{j\in J_b}\left(y_{\mathit{jbat}}{w}_{\mathit{jbsa}}\right)\left(g_{\mathit{jba\epsilon s}}-L_{\mathit{\epsilon t}}^{\mathit{\hspace{0.17em}}}\right)+{{\mathit{d}}_{\mathit{\epsilon }\mathrm{ts}}}^{-}\ge 0,\mathrm{\forall }\mathit{\epsilon }\in \mathit{E},\mathit{t}\in \mathit{T},\mathit{s}\in \mathit{S}\end{array}$(18)

In order to guarantee that grade requirements for different elements $\mathit{\epsilon }\in \mathit{E}$ are achieved, grade blending constraints are added to the formulation (Equations 17 and 18). These constraints allow deviations ${{\mathit{d}}_{\mathit{\epsilon }\mathrm{ts}}}^{+}$and ${{\mathit{d}}_{\mathit{\epsilon }\mathrm{ts}}}^{-}$ from upper and lower bounds respectively for each scenario$\mathit{s}\in \mathit{S}$.

(19) $\begin{array}{c}{y}_{\mathit{jbat}}\in \left\{0,1\right\},\mathrm{\forall b}\in \mathit{B},\mathit{a}\in A_b,\mathit{j}\in J_{\mathit{ba}},\mathit{t}\in \mathit{T}\end{array}$(19)

(20) $\begin{array}{c}{z}_{\mathit{ba}}\in \left\{0,1\right\},\mathrm{\forall b}\in \mathit{B},\mathit{a}\in A_b\end{array}$(20)

(21) $\begin{array}{c}{\mathit{\psi }}_{\mathit{dlbt}}^{\mathit{drift}},{\mathit{\psi }}_{\mathit{dlbt}}^{\mathit{drif}{\mathit{t}}^{\ast }}\ge 0,\mathrm{\forall b}\in \mathit{B},\mathit{l}\in L_b,\mathit{d}\in D_l,\mathit{t}\in \mathit{T}\end{array}$(21)

(22) $\begin{array}{c}{\mathit{\psi }}_{\mathit{clbt}}^{\mathit{crosscut}},{\mathit{\psi }}_{\mathit{clbt}}^{\mathit{crosscu}{\mathit{t}}^{\ast }}\ge 0,\mathrm{\forall c}\in C_l\mathit{b}\in \mathit{B},\mathit{l}\in L_b,\mathit{c}\in C_l,\mathit{t}\in \mathit{T}\end{array}$(22)

(23) $\begin{array}{c}{{\mathit{d}}_{\mathit{hts}}}^{\mathit{haul}}\ge 0,\mathrm{\forall h}\in \mathit{H},\mathit{t}\in \mathit{T},\mathit{s}\in \mathit{S}\end{array}$(23)

(24) $\begin{array}{c}{{\mathit{d}}_{\mathit{ts}}}^{\mathit{proc}}\ge 0,\mathrm{\forall t}\in \mathit{T},\mathit{s}\in \mathit{S}\end{array}$(24)

(25) $\begin{array}{c}{{\mathit{d}}_{\mathit{\epsilon }\mathrm{ts}}}^{+},{{\mathit{d}}_{\mathit{\epsilon }\mathrm{ts}}}^{-}\ge 0,\mathrm{\forall }\mathit{\epsilon }\in \mathit{E},\mathit{t}\in \mathit{T},\in \mathit{S}\end{array}$(25)

Equations 19(19) $\begin{array}{c}{y}_{\mathit{jbat}}\in \left\{0,1\right\},\mathrm{\forall b}\in \mathit{B},\mathit{a}\in A_b,\mathit{j}\in J_{\mathit{ba}},\mathit{t}\in \mathit{T}\end{array}$(19) –25(25) $\begin{array}{c}{{\mathit{d}}_{\mathit{\epsilon }\mathrm{ts}}}^{+},{{\mathit{d}}_{\mathit{\epsilon }\mathrm{ts}}}^{-}\ge 0,\mathrm{\forall }\mathit{\epsilon }\in \mathit{E},\mathit{t}\in \mathit{T},\in \mathit{S}\end{array}$(25) refer to integrality and non-negativity constraints.

2. Case study – application at an operating underground copper mine

In this section, an application of the proposed method at an operating underground copper (Cu) mine is presented, where gold (Au) and uranium (U₃o₈) are secondary elements. First, the results for the underlined integrated stochastic framework are analysed. Then, these results are compared to a stochastic sequential framework in which a stope design is given as an input to the method, fixing the available stopes to be scheduled, as well as its types.

As the main input, 10 geostatistical simulations of Cu, Au, and U₃O₈ in a grid of size 5 m × 5 m x 5 m, for a mining zone of the considered mineral deposit are used. Previous studies show through a sensitivity analysis of the stochastic optimisation of the long-term mine planning, that 10–15 simulations are sufficient to produce stable results. This conclusion is attributed to the support-scale effect, once a large number of blocks are grouped to generate the schedule of a given production period [59–61]. Figure 8 shows two realisations for copper, gold and uranium grades. Mine accesses (i.e. ramp and decline) and ventilation systems are also fixed inputs. Figure 9 shows the available infrastructure in the mine and the defined mining direction. In addition, geometrical parameters for the shapes of the stopes were provided by the mining company operating the mine. The maximum and minimum dimensions of stopes are considered given geotechnical and drilling equipment requirements. Two possible mining zone configurations associated with these allowable shapes are considered (Table 7). These configurations account for the position of the access point at each sublevel according to the available ramp design. Three option types are used as an input, meaning that all stopes will have the possibility of being primary, secondary, or tertiary.

Figure 8. Realizations of (a) copper, (b) gold, and (c) uranium grades in a grid of 5m x 5m x 5m.

Figure 9. Plan view of the developed infrastructure.

Table 7. Stope geometrical parameters.

Parameter	Value/Description
Minimum dimensions	10m x 30m x 30m
Maximum dimensions	30m x 30m x 150m
Configuration 1	15m x 30m x 40m (3,240 potential stopes)
Configuration 2	30m x 30m x 80m (810 potential stopes)

For this case study the stopes are considered to be blasted and extracted bottom-up, with the pattern of extraction represented in Figure 5, and are subsequently backfilled with cemented aggregate fill (CAF), which is not considered a limiting feature in mine production. Thus, the same mining cost is used for all stope types. A horizontal development capacity is considered in terms of the maximum length that can be developed. A single haulage system is available with its maximum capacity constraining ore production. Once the sources of uncertainty considered do not directly affect the ore tonnage production, the mining capacity constraints are modelled as hard constraints. On the other hand, uncertainty in terms of copper, gold, and uranium grades is taken into consideration, thus, penalty costs for deviations from the minimum and maximum grades for these three elements’ requirements are applied and are discounted throughout the years in order to manage the geologic risk. Table 8 displays the technical and economic parameters used in the optimisation of the copper mine. It is worth noting that a fixed unit mining cost of 50 $/t is considered. This cost incorporates drilling, blasting, mucking and other fixed yearly costs such as ventilation costs and backfilling that were scaled in terms of the yearly production rate.

Table 8. Technical and economic parameters used as input in the optimisation.

Parameter	Value
Cu price ($/t)	8,500
Au price ($/ozt)	1,200
U3O8 price ($/t)	71,904
Economic discount rate	10%
Geologic discount rate	10%
Processing recovery Cu (%)	94%
Processing recovery Au (%)	70%
Processing recovery U3O8 (%)	70%
Mining cost ($/t)	50
Processing cost ($/t)	13.5
Haulage cost ($/t*km)	5
Drifts development cost ($/m)	12,000
Density (t/m3)	3.2
Block tonnage (t)	400
Mining capacity (Mt/y)	3
Drift development capacity (m/y)	5,000
Minimum copper mill head grade (%)	1.8
Minimum gold mill head grade (g/t)	1.0
Maximum gold mill head grade (g/t)	2.0
Minimum uranium mill head grade (g/t)	420
Maximum uranium mill head grade (g/t)	600
Penalty cost for deviations below minimum Cu mill head grade ($/unit)	100
Penalty cost for deviations below minimum and above maximum Au mill head grade ($/unit)	100
Penalty cost for deviations below minimum and above maximum U3O8 mill head grade ($/unit)	10

2.1. Results of the integrated stochastic optimization

The outputs of the proposed integrated stochastic approach are presented next. Figure 10 shows the optimal stope shapes, types, and sequence of extraction throughout a 12-year life-of-mine. The selection of the mining zone configuration represents the trade-off between selecting grades and development costs. The optimal configuration 2 (Table 7) guarantees a lower development cost, while deviations from grade production targets are well managed. The green curves in Figures 11 and 12 show the risk profiles of the integrated stochastic stope design and schedule in terms 10%, 50%, and 90% probabilities (i.e. P10, P50 and P90). An NPV of 3.53 B$ (Figure 11a) considering the P50 and a cumulative development cost of 157 M$ (Figure 11b) are observed. Also, the ore production follows the maximum mining capacity (Figure 11c) and the Cu, Au, and U₃O₈ grades have small deviations from the defined bounds and tend to decrease through the years showing the effect of geological risk management (Figure 12). These results are compared to a stochastic sequential framework in the next subsection.

Figure 10. Integrated stochastic optimization outputs from left to right: the stope types option selected and the extraction sequence.

Figure 11. Risk profiles of the integrated (green curves) and sequential (black curves) stochastic frameworks: a) NPV; b) cumulative horizontal development cost; c) ore tonnage.

Figure 12. Risk profiles of the integrated (green curves) and sequential (black curves) stochastic frameworks: a) cumulative Cu content; b) mill Cu head grade; c) cumulative Au content; d) mill Au head grade; e) cumulative U₃O₈ content; f) mill U₃O₈ grade.

The formulation was implemented in C++ on Visual Studio 15 and solved with CPLEX v.12.8.0. The present application is compounded by 117,684 binary decision variables and 2,172,580 constraints. Using a standard personal computer with six cores and 32 GB RAM, the preprocessing and optimisation steps took approximately 24 hours and were constrained by memory allocation limitations.

2.2. Comparison between the integrated and the sequential stochastic approaches

To show the importance of jointly optimising the stope design and extraction sequence, the proposed method is compared to the stochastic sequential approach. In this sequential approach, the preprocessing steps generate the possible mining zone configurations with the respective allowable stope shapes. From these configurations and possible stopes, a stope design that maximises the undiscounted cashflow, regardless of adjacency constraints and the horizontal development costs, is first generated. The selected stope locations and shapes are used as the inputs to generate the extraction sequence using the SIP formulation presented in Section 2.

The black curves in Figures 11 and 12 show the results of this sequential approach. The same technical, geometric, and economic parameters presented in Tables 7 and 8 are used. It is seen in Figure 11a, that the NPV of the integrated approach is 3.54 B$ while in the sequential approach it is 3.35 B$, considering the P50 in both cases. Therefore, the integrated approach is expected to produce a NPV 190 M$ higher than the sequential approach. Although less stopes are mined when the sequential approach is used, the development cost for the sequential approach is almost three times higher than that in the integrated approach (Figure 11b). Therefore, many stopes that have a positive impact economic value when the cumulative development costs are not considered become uneconomical and inaccessible when this information is actually taken into account. Figure 13 shows the comparison of the extraction sequence of the integrated and sequential approaches. In Figure 13b, the black wireframe shows the initial stope design generated with step-wise approach, that is noted to be physically different to the one produced by the integrated approach. In addition, it is seen that several stopes included in the initial design in Figure 13b are not included in the schedule, showing that, in fact, the effective development costs are critical in the economic value of stopes. It can be also noted that the stopes chosen in the sequential approach are smaller (i.e. configuration 1 in Table 7 – Stope geometrical parameters), generating more selectivity in terms of grades, but producing comparable metal contents for the three elements, Cu, Au, and U₃O₈ and meeting grade blending requirements (Figure 12). The integrated approach, however, shows a higher metal production in early periods, which has a greater positive impact on the NPV, due to the time value of money.

Figure 13. Comparison of the extraction sequence of the a) integrated approach and b) sequential approach (the wireframe corresponds to the stopes selected in the stope design).

3. Conclusions

A new mathematical programming formulation for the integrated stochastic optimisation of stope design and long-term mine production scheduling is presented, along with its application at an operating underground copper mine. The method is developed based on the sublevel longhole open stoping (SLOS) with backfilling underground mining method and overcome the limitations of previously proposed approaches that are tailored to specificities of other mining methods. The SLOS variation considered in the present study follows the assumptions and parameters derived from an existing operational mine. The proposed method generates jointly the stope boundaries and the extraction sequence assuming that all stopes follow the same geometrical parameters. In addition, a pattern of extraction and a mining direction define adjacencies among stopes and a mining cycle can be completed in one year, which defines the period of extraction.

The proposed two-stage stochastic integer programming (SIP) maximises the NPV, as well as considers metal prices for different elements, mining costs for different types of stopes, horizontal development costs and haulage costs for the different systems available, while also minimising the risk of not meeting production targets. The output of the optimisation is an operational selection of mining zone configuration that defines the stope shapes and respective types, as well as the extraction sequence of stopes that respects the optimal adjacencies and defined mining direction, generating a risk resilient production schedule.

The proposed method is applied to an underground mine that has copper as the main element, as well as gold and uranium as the secondary elements. The mining zone studied has an available ramp and a ventilation system that defines a mining direction. It follows a pattern of extraction of primary, secondary, and tertiary stopes that defines the adjacencies. The integrated framework is compared to a sequential stochastic framework, in which a design that defines the mining zone configuration and stopes that maximise the undiscounted cashflow is used as a fixed input to the proposed SIP. Firstly, it is seen that physically different stope designs and, consequently, different extraction sequences are produced for these two different frameworks. Furthermore, the integrated approach shows an NPV that is 6% higher than the step-wise approach. The horizontal development costs of the sequential approach are shown to be clearly higher than in the integrated approach. This difference can be attributed mainly to the fact that the actual development costs cannot be considered when generating the initial stope design, which, thereby, limits the decisions in the extraction sequence for the step-wise framework. Therefore, the importance of having an integrated method that is capable of exploiting the relationships between the optimisation components is validated.

A case study that accounts for multiple sources of uncertainty simultaneously is a topic for future research, once the presented SIP can directly accommodate commodity price uncertainty. The application and the several aspects of the method are built based on the assumption of a large underground mine divided into mining zones. In the presented case study, only one mining zone schedule was optimised due to the complexity associated with the orebody size and computational limitations. Considering that multiple mining zones can be mined simultaneously and can contribute to the production of ore that feeds the processing plant, an application that optimises multiple mining zones simultaneously could be considered. In addition, the simultaneous optimisation for a mining complex that assumes the existence of multiple mines, stockpiles, and processing streams, with critical considerations in terms of vertical development costs is an extension for future work. Furthermore, computational efficiency limitations are seen when a commercial solver is used, which restricts the model’s application to larger problems. Thus, a metaheuristic solver should be implemented in future developments.

Acknowledgments

This work is funded by the National Science and Engineering Research Council of Canada (NSERC) CRD Grant 500414-16, NSERC Discovery, Canada Grant 239019, the COSMO Stochastic Mine Planning Laboratory and mining industry consortium (AngloGold Ashanti, Anglo American, BHP, De Beers, IAMGOLD, Kinross, Newmont Mining, and Vale), and the Canada Research Chairs Program.

Disclosure statement

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