**PROFESSIONAL TECHNICAL AND SCIENTIFIC PAPERS**

**Clustering-based iterative approach to stope layout optimization for sublevel stoping**

**Y.A. Sari; M. Kumral**

Department of Mining and Materials Engineering, McGill University, Canada

]]>

**SYNOPSIS**

Underground mining operations tend to have higher operating costs than surface mines. When metal prices decrease, profitability is jeopardized due to the high costs. Therefore, mining management harnesses new practices that increase operational efficiency. One way to manage this challenge is to invest in new mine planning practices. Stope layout optimization as a part of underground mine planning aims to identify a portion of the orebody in the form of production volumes (stopes) to maximize profit under roadway and stope dimension constraints. In this paper we propose a novel approach based on identifying ore-rich areas of the deposit and prioritizing their extraction through an iterative heuristic clustering approach. The proposed approach is compared with and validated by an exact method through a small mining example. The heuristics produced nearly identical results in a very short time. Finally, a case study was carried out using a larger data-set. The cluster-based iterative approach generated near-optimal stope layouts in a computationally effective manner.

**Keywords**: underground mining, iterative optimization, stope layout planning, sublevel stoping.

**Introduction**

Surface mines have long constituted the majority of worldwide mining operations. The main reasons are that

(1) The new excavation and slope monitoring technologies allow the surface mines to go deeper due to the lowered overburden removal costs

(2) Surface operations are more selective, resulting in less dilution and loss

]]> (3) Surface extraction techniques usually involve fewer safety concerns compared to underground mining.

Underground mining has complicated engineering considerations such as rock stress calculations, safety, and air distribution through appropriate ventilation. Current circumstances favour underground mining because of the following reasons:

➤ Many surface deposits have been depleted

➤ Overburden waste rock in surface mines is much less or even does not exist in underground mining; lower stripping ratio generates higher profits as it lowers the mining and waste handling costs

➤ Underground mines have less environmental impact than surface mines.

Mining projects are already very risky due to the high uncertainty related to grade distribution and volatile commodity prices (Maseko and Musingwini, 2019; Sauvageau and Kumral, 2017), thus emphasizing the importance of operational optimization. Although underground mining is considered as being governed by rock mechanics constraints, with no room for optimization, the recent prevalence of underground mining techniques increases the importance of computer-aided tools for planning and layout optimization, in order to maximize the profit and minimizie the environmental impact.

Sublevel stoping is an unsupported underground mining technique that is typically used when the orebody is massive, steep, thick, and large in size (Hartman and Mutmansky, 2002). Additionally, the rock substance strength should be medium to strong (Nicholas, 1981). The orebody is accessed through underground access roads (levels) from the shaft, and between levels rectangular extraction areas (stopes) are outlined and accessed through sublevels. When development is completed, the stope is drilled from several access points and the ore blasted. The comminuted rock collapses to the bottom of the stope, from where it is carried from the drawpoints to the shaft by haulage trucks. The minimum and maximum stope length, width, and height are determined by the geological engineer according to the rock characteristics.

Sublevel stoping comprises two main problems: stope layout planning and stope sequencing. Stope layout planning entails positioning the sublevels and stopes, as well as deciding stope dimensions, in such a way that profit is maximized. Stope sequencing aims to determine the order of mining of the stopes, taking into consideration the mine stability, equipment transportation, and net present value (NPV) maximization.

The majority of current computerized techniques approach the stope layout planning problem by partitioning the orebody into a block model, in which each block is estimated or simulated based on drill-hole samples. The grades are then converted to economic values using parameters such as ore price, and mining and processing costs. This conversion facilitates the evaluation of the blocks because the profitability of extracting a particular block can be directly recognized. The blocks are selected to be extracted such that the profit is maximized and rock stability constraints are not violated. These constraints are explained in detail in the formal problem definition section of the paper.

]]> In this paper we introduce a new clustering-based algorithm that is inspired by the practical approach implemented by mining engineers to plan the stope layout. The originality of this paper is two-fold: (1) A new formulation of the stope layout optimization problem is proposed, and (2) this problem is solved by a new three-stage approach that forms a block model, places the sublevels, then the stopes using clustering heuristics.This paper is organized as follows: in the next section, different approaches to stope layout planning are discussed. Then, the stope layout optimization problem is then defined mathematically with a new formulation. The proposed heuristic is presented in detail and the corresponding mathematical model given. The approach is tested with two case studies and a comparison with the mixed-integer programming model formulation is provided. Finally, the findings are summarized and conclusions drawn.

**Literature review**

Sublevel stope layout planning is more difficult than surface mine planning as it involves a similar amount of decision variables but is subject to more constraints. Exact methods ideally yield optimal results, but in practise they are not able to handle large deposits as they either run out of memory or simply take a very long time to solve, which is impractical unless heuristic approximations are made. Typically a block model consists of thousands to millions of blocks. Generally, to deal with the large number of decision variables, exact methods are modified such that they are faster but non-optimal.

Jalali and Ataee-pour (2004) presented a dynamic programming approach for vein-type orebodies based on a modification of the Johnson and Sharp's (1971) dynamic programming algorithm for open pit layout optimization. Instead of taking into consideration the slope constraints as in open pit mining, a maximum variation of the elevation of both the floor and the roof from one column to the next is allowed for draw control. The algorithm provides two-dimensional (2D) solutions by combining column economic values that are perpendicular to the vein direction. Ovanic and Young (1995) introduced the branch and bound technique to optimize the start- and endpoints at each row of blocks. Two piecewise linear cumulative functions at each row, representing the physical location of the start- and end-points, are declared and the problem solved using a mixed integer programming (MIP) approach. This approach, known as SOS2 (type-two special ordered sets), also called separate programming, allows at most two adjacent ordered sets of variables to be non-zero. Although row-by-row this approach is optimal, this does not guarantee overall optimality.

Copland and Nehring (2016) proposed a MIP model to simultaneously optimize stope layout and sequence. Bai, Marcotte, and Simon (2013) followed a very different approach, where the model is defined on a cylindrical coordinate around the initial vertical raise. Blocks are converted into nodes and a source and a sink are added to the model, which is then solved by the maximum flow approach. Having a vertical raise limits the optimality in cases where the orebody is inclined, resulting in the inclusion of too much waste, particularly in cases where the deposit is larger and more than one raise will be required.

The octree division approach (Cheimanoff, Deliac, and Mallet, 1989) recursively divides the three-dimensional (3D) model into two equal parts in each dimension, resulting in eight subvolumes, and includes the subvolume in the final stope layout if it is valuable throughout. Since partial stopes are not allowed during the optimization and the algorithm works by equal division of volumes, stope locations are checked only where the minimum stope dimension is a proper divisor. The downstream geostatistical approach (Deraisme, de Fouquet, and Fraisse, 1984) uses dynamic programming on the 2D sections of the drilling and blasting data to minimize dilution.

Heuristic methods, being fast and practical, are most commonly used in the industry. In the floating stope algorithm (Alford, 1995; Alford, Brazil, and Lee, 2007) a potential stope with minimum dimensions is floated through the block model. All economic stopes are included in the final stope layout. However, a problem arises when two stopes overlap. Also, different results are obtained depending on the starting point of the floating process. Similarly, the maximum value neighbourhood algorithm (Ataee-pour, 1997, 2004) examines the neighborhood of each block in sequence and from all possibilities, the neighbourhood with maximum economic value is included in the final stope. The preference-based profit maximization approach by Topal and Sens (2010) creates a list of all possible stopes and chooses from them according to the user preference. Wang and Webber (2012) implemented a two-stage approach where the rings that do not contain ore are filtered out and the design is completed manually. Sandanayake, Topal, and Asad (2015) developed an algorithm that incorporates stope size variation by aggregating the mining blocks into a possible set of stopes, then modifying the stope attributes. Nikbin *et al. *(2019) developed two heuristic algorithms, namely the greedy and iterative enumeration algorithms, to overcome the shortcomings of exact algorithms. Later, Nikbin *et al. *(2020) combined the greedy algorithm with dynamic programming to obtain a hybrid approach. Sari and Kumral (2020) proposed a dynamic programming approach with an optional greedy heuristic that provides fast and efficient results if the stope size is fixed. Although it can potentially work with varying stope sizes, this would increase the solution time. Sari and Kumral (2019) extended this approach to polymetallic mines with pillars and introduced an ultimate stope limits formulation, which is analogous to ultimate pit limits in open pit mine planning.

Recently, metaheuristic methods have been implemented to determine stope layout. Villalba Matamoros and Kumral (2019a, 2019b) used genetic algorithms to optimize stope layout under grade uncertainty. Hou *et al. *(2019) also used genetic algorithms to optimize stope boundaries and access layout at the same time. Foroughi *et al. *(2019) proposed a non-dominated sorting genetic algorithm for the integrated optimization of stope layout design and production scheduling. This is a multi-objective optimization model where NPV and metal recovery are optimized at the same time.

**Formal problem definition**

Ore deposits are conceptually divided into uniform rectangular grids (mining blocks). In a deposit, each potential stope has a certain economic value as it contains a unique set of blocks, each with an estimated or simulated ore grade. Stope layout planning involves placing non-overlapping three-dimensional stopes in a deposit within constrained dimensions such that the portion of the deposit that maximizes profit is selected. An additional constraint arises from the construction of the access roads below and above the stopes, known as sublevels. Stopes, should be vertically aligned.

Throughout the paper, four mathematical models are presented that use the notation given in Table I. In addition to the list of common notations, the notation unique to each mathematical model is presented following the model.

A new formulation of the stope layout planning problem is presented in the following equations.

Maximize:

]]>Subject to:

In this model, *X,Y,Z *are the number of blocks in the X, Y, and Z directions, respectively, *ij,k *are sets of starting coordinates of all valid stopes in the block model, θ is the stope height that is determined previously, *q _{x}*and

*q*represent offsets from the starting coordinates,

_{y}*y(ij,k,%,*ψ

*,Q)*is the binary decision variable that determines the extraction of a stope at the coordinate

*ij,k*with the sizes χ,ψ,θ, and

*v(ij,k)*is the economic value of the block at coordinates

*ij,k.*The objective function in Equation [1] maximizes the total economic value of the stopes that are selected to be extracted by the model. Equation [2] ensures that only one size can be accepted per stope. Equation [3] checks for overlapping stopes in the X-Y direction and allows only one of the overlapping stopes to be selected. The additional constraint for preventing stopes from overlapping vertically is given in Equation [4]. However, this constraint is different from the X-Y direction overlap constraint given in Equation [3] because the overlap should be avoided not only directly on the stope but also throughout the Z level. As explained earlier, this is important for stable sublevel formation. In this model, the sublevel construction cost is assumed to be included in the mining cost. In common practice, due to the high cost of building sublevels, the distance between sublevels is kept as large as possible and constant within each geological domain. In this case,

*z*can be set equal to

_{a}*z*to simplify the problem.

_{b}The stope layout design problem resembles a 3D container loading problem (Zhao *et al., *2016; Zhu, Balakrishnan, and Cheng, 2018), and it is most similar to a capacitated clustering problem (CCP) in operations research where a specified number of clusters are formed from a set of elements with certain weights. The total cluster weight is constrained within lower and upper limits. Each pair of elements has a predefined benefit that contributes to the objective value only if the pair is in the same cluster, and the objective is to maximize the overall benefit (Osman and Christofides, 1994). A CCP can be transformed into a stope layout design problem by representing possible stopes as elements with weights of unity and setting the maximum weight constraint as infinite and the minimum weight constraint as zero. Additionally, the number of clusters should be equal to the number of possible stopes, and the predefined benefit of a pair of elements should be assigned in proportion to the sum of their economic values or be infinitely negative if two stopes intersect or align in an overlapping fashion on the Z axis.

When the CCP is solved using these inputs, the stopes in the cluster that provide the highest benefit are to be extracted. The reasoning behind this is that stopes that can be selected together (not overlapping or intersecting) will be forced to be in the same cluster because the overall benefit increases when the number of pairs increases. If two stopes overlap or intersect, their benefit is negative. Thus, they will not be in the same cluster. As there are as many clusters as stopes, these stopes will be distributed in different clusters. Also, the way the weight capacities are set allows clusters to contain from zero to an infinite number of stopes, not limiting the number of stopes in the mine plan.

**Heuristic methodology**

When an optimization problem contains an excessive number of variables and constraints, the solution time becomes impractically long, thus heuristic methods are preferred (Zhu and Lei, 2018). As can be observed from the previous section, the vertical alignment condition increases the number of constraints drastically. However, it can be anticipated that with regard to the minimum and maximum stope height constraints, there are relatively few possible combinations of sublevels;. hence a heuristic that segregates sublevel design and stope layout design is proposed. When sublevel design has been carried out previously, the alignment constraints can be eliminated (decreasing the number of constraints) and the problem is divided into smaller problems (decreasing the number of variables). This is realized with the aid of a clustering heuristic in which the value of a block depends not only on its grade but also on the grades of neighbouring blocks. The heuristic aims to identify the ore-rich sections of the deposit and extract these volumes as stopes, taking into account the structural feasibility. With this information, sublevels are determined, followed by the decision regarding stopes in between the sublevels. The three-stage heuristic approach is summarized in Figure 1. Table I provides the notation used in this section.

]]>

*Preparing the model*

The input consists of a block model with grades. Ore-rich sections

are identified through analysis of the block grades. Each block is assigned a score according to the grade of a given block and the grades of the surrounding blocks. Each set of surrounding blocks of a given shape is called a layer. If only the immediate blocks that are adjacent to the given block are considered, the depth of the surrounding layer becomes 1. The number of surrounding layers can be increased by considering the next set of blocks adjacent to the previous layer and their grades are added, optionally multiplying by a discount factor at each increasing level. This is demonstrated in Figure 2, in which three layers are framed. This way, the heuristic mimics the clustering approaches. With this heuristic score, each block contains information about its grade and strategic location. The depth of the layer surrounding the block is customizable in the program. The number of layers is closely related to stope size. If the depth is set to a feasible stope size, each block will contain the heuristic score for the stope where the block is centred. If there is more than one valuable metal in the block, secondary metals are converted in terms of the primary metal by using the equivalent grade (Equation [5]) and the grades are converted using Equation [6]. It can be observed that for the first metal, *F *= 1, resulting in *g**eq *= g + Σ g χ *F _{eq}. *where i ∈

*m*and

*i*≠1.

]]>

Block scores are calculated according to Equation [7] After the block scores have been calculated, the block scores on each level of the deposit are computed by adding the scores of the blocks on the corresponding level as per Equation [8] and saved as sublevel candidates.

*Sublevel design*

The selection of sublevels is a significant stage in the planning as it influences the succeeding decisions. In practice, ore-rich sublevéis are selected where the ore concentration is high such that the stopes will also have high average grades. Considering that sublevels are access roads, and are also extracted, sublevels are generally selected in ore-rich regions as this paves the way to more profitable potential stopes. This heuristic is inspired by this practical approach and attempts to optimize it by making use of computational tools. In the current stage, it is important that the height of the sublevels should satisfy the minimum stope height constraint.

Mathematically, sublevel design can be expressed as follows: Maximize:

where *k *+ θ __<__ Z

where *k* + 2θ ≤ Z

In Equations [9-11], *k *represents the starting block of a stope in the Z direction, θ represents the height of the stope, *y(kfi) *is the binary decision variable that selects a sublevel and the maximum heights of stopes accessed from that sublevel, and *v(ij,b) *is the economic value of the block at coordinates *ij,*δ*. *The objective function (Equation [9]) maximizes the total value of the sublevel by summing the economic values of each block within the range of the height of the sublevel. Equation [10] ensures that only one maximum stope height can be accepted below each sublevel. Equation [11] ensures that if a sublevel with a certain stope height is selected, an overlapping level cannot be selected. As a result of this stage, a combination of sublevels will be output. The best scoring solution of the candidate sublevel sets is chosen. This stage of the approach resembles the layer-building heuristics for the 3D container loading problem (Zhao *et al., *2016). As mentioned previously, due to the high cost of building sublevels an assumption may be made to simplify the problem: the distance between sublevels may be kept as large as possible and constant within each geological domain. For this reason, the mathematical model can be simplified as follows:

Maximize:

where *k *+ *z _{b}< Z*

Subject to:

The heuristic algorithm that selects the sublevels is designed as follows:

]]>➤ Only the levels that can possibly satisfy the stope height constraints are selected as candidates to speed up the search. The non-satisfying levels are eliminated.

➤ The remaining levels are ranked according to their scores.

➤ A candidate solution is created by taking the first level in the ranked list, then adding the following levels in the list as long as the candidate solution is feasible. If the addition of a level makes the candidate solution infeasible, the next level is added until the list is exhausted. The feasibility is tested by verifying that each level has at least the minimum stope height.

➤ The overall score of the candidate solution is calculated by averaging the scores of levels in the solution and multiplying by the number of levels that can be accessed through the sublevels.

➤ If this is the first calculated score, or it is the highest score so far, it is stored as the current best solution. Otherwise, the candidate solution is deleted.

➤ The first level in the ranked list is deleted. If there are no more levels in the list or the number of formed level combinations is equal to α, the algorithm is terminated. Otherwise, the algorithm returns to step 3.

The sublevel design algorithm is also illustrated in the flow chart given in Figure 3.

]]>

*Stope layout design*

Given the sublevéis with the above approach or manually, the stopes are planned. In mixed integer linear programming terms, stope layout design between sublevels can be expressed as follows.

In this model, *ij,k *are sets of starting coordinates of all valid stopes in between sublevels, *Z _{1}*and

*Z*are beginning and ending coordinates of the stopes, respectively, in the Z direction, θ is the stope height that is determined previously,

_{2}*q*and

_{x}*q*represent offsets from the starting coordinates,

_{y}*y(ij,k,x,*Ψ

*>*θ

*)*is the binary decision variable that determines the extraction of a stope at coordinate

*ij,k*with sizes χ,ψ,θ, and

*v(ij,k)*is the economic value of the block at coordinates

*ij,k.*The objective function in Equation [14] maximizes the total economic value of the stopes that are selected to be extracted by the model. Equation [15] ensures that only one size can be accepted per stope. Equation [16] checks for overlapping stopes in the X-Y direction and allows only one of the overlapping stopes to be selected.

This is also realized through a similar but iterative and metaheuristic-like approach. Each level is considered separately and the stopes at a level are selected according to the following procedure.

➤ The height of the stopes is settled by observing the level height. If it is within acceptable limits, stope heights are set to the level height. Otherwise, the maximum stope height is set as the current stope height.

➤ A list of blocks in the level is established and sorted according to their scores s

]]> ➤ Similar to the sublevel selecting algorithm, a candidate stope combination solution is created by taking the first block of the ranked list. All stope size combinations are evaluated for their economic value and feasibility. The feasibility test consists of testing whether the stope is out of block model bounds and if the stope overlaps with another stope already selected. Within the set of stopes that are feasible, the stope that yields the highest economic value is selected. As long as the feasibility constraints are sustained, the stope is added. This process is continued until the list is exhausted._{w}, from high to low.➤ The overall economic value of the stope combination is calculated by adding the economic value of each stope in the combination.

➤ If the economic value of the combination is positive and is the highest economic value so far, it is stored as the current best combination. Otherwise, the combination is deleted.

➤ The first block of the ranked list is deleted. If there are no more blocks in the list or the number of formed stope combinations is equal to β, the algorithm advances to next step. Otherwise, the algorithm returns to step 3.

➤ This is the metaheuristic step. The iteration with the highest economic value is compared to the most recent iteration. If there is an improvement, the ordering of the most recent iteration is kept. Otherwise, the ordering is kept with a probability. This probability decreases as the number of iterations increases. If there is no improvement for

nconsecutive iterations, wherenis determined empirically, the algorithm terminates.➤ A random change is made to the ordering of the list and the algorithm returns to step 3.

The stope layout design algorithm is also illustrated in the flow chart given Figure 4. The stope combination selection process is condensed into one step as it is very similar to the sublevel design algorithm.

]]>

In this stage of design, although block scores *s _{xJZ}*are used when forming the priority list, the final decision is made according to the economic value of the design. This is preferred because the overall objective is to optimize profit. This process may be repeated with different sublevels to start with. As the number of sublevel combinations is limited due to stope height constraints, the number of repetitions will also be low and the result will be closer to optimal.

This strategy can be very effective when the ore-rich regions are in clusters and unevenly dispersed throughout the deposit as it conveniently prioritizes ore-rich areas, which is common in most deposits. If this is the case, the upper elements in the list will be dominant over the subsequent elements, increasing the possibility that the right combination will contain the upper elements. Another advantage of this approach is that given its modular structure, sublevels may be defined by an engineer, or if there is an existing development it can be defined directly in the program and the sublevel design stage above can be omitted.

**Case studies**

*Case study 1*

The heuristic approach was implemented in C++ and tested on an underground polymetallic gold-copper mine with 125 000 blocks. The data-set contains two grades for each block (gold and copper) that are the averages of multiple realizations generated by sequential Gaussian simulation. The mining plan incorporates one mineral processing plant and one waste dump. The parameters regarding the project are given in Table II.

]]>

The economic value of an extracted block is calculated according to Equation [17], where *g _{m}*is grade,

*p*is price,

_{m}*R*is the recovery of mineral m,

_{m }*C*is processing cost,

_{p}*C*is mining cost, and

_{r}*t*is tonnage. The tonnage is calculated by multiplying the specific gravity by the block volume. It is assumed that all mined blocks are processed and un-extracted blocks have zero value. The resulting plan can be visualized using SGeMS software (Remy, Boucher, and Wu, 2009) as in Figure 5. Each stope is illustrated with a colour that corresponds to the average grade in the stope. Representative sublevels are shown in grey. The heuristic approach was able to successfully identify the ore-rich areas in the deposit and generate a stope layout plan. It can be observed from Figure 5 that the minimum average stope Au grade is 0.376 g/t, which can be considered as the stope cut-off grade for this operation.

Figure 6 demonstrates the effect of program-related parameters α,β and *l *(definitions given in Table I) on the overall profit of the mine. To observe this effect, the program was run multiple times with different parameter values. α and β were tested in the range 5-20 and *l *was tested in the range 3-4. Inspecting the figure, it can be inferred that in this case study, β influenced the profit the most and *l *did not have an effect. Also, increasing α and β increased the profit until α and β reach about 10. Above this value, the profit reached a plateau. This indicates that the heuristic approach is successful in ranking the more promising sublevels and stopes before unfavourable possibilities.

The average running time of the program for this case was 15 minutes 42 seconds on a MacBook Pro 2015 with 2.7 GHz Intel Core i5 processor and 16 GB memory. To decrease ore dilution and loss, the block size can be reduced. However, this would increase the number of decision variables, and hence the solution time.

To further evaluate the performance of the heuristic approach, case study 1 was re-run on a portion of the same deposit with dimensions 15 χ 15 χ 15 blocks and maximum frame size as 50, 50, 50 m. The linear programming model was formulated in a Zimpl programming environment (Koch, 2006) and the same case study was solved with the linear programming model using CPLEX to compare the outcomes. The optimal mine value obtained by the linear program was $382 037 496, and by the heuristic approach was $377 561 000 which is a difference of 0.1%. On the other hand, the runtime of CPLEX was 61 hours 23 minutes on a Dell Precision T3610 workspace with Intel Xeon E5-1620 3.70 GHz processor, whereas the heuristic approach took 2 minutes 4 seconds with the same computer used in the previous case study.

*Case study 2*

A second case study was carried out to test a slightly modified version of the heuristic algorithm. In this version, three changes were made:

(1) The distance between sublevels was kept constant

(2) Sublevel building cost was added

]]> (3) Internal waste was allowed to plan diluted stopes as opposed to typical mining stopes.

Instead of assuming the entire volume of the stope will be extracted, it is presumed that blasting can be adjusted such that internal waste can be left on the roof and floor of the stope. In addition, a small block size is used to decrease dilution/loss.

This algorithm was tested on an underground gold mine with 134 400 blocks. The dataset contains a gold grade for each block that is the averages of multiple realizations generated by sequential Gaussian simulation. The mining operation plans one mineral processing plant and one waste dump. The parameters regarding the project are given in Table III. The economic value of an extracted block is calculated in a fashion very similar to that in the previous case study. The only difference between the calculations is that in case study 2 the mining cost does not incorporate the sublevel building cost - this is extracted from the economic value of the mine in proportion to its length.

The running time of the program for this case was 96 minutes 17 seconds on a MacBook Pro 2015 with 2.7 GHz Intel Core i5 processor and 16 GB memory. The resulting plan for this case is visualized from two different perspectives using SGeMS software in Figure 7, where each stope is illustrated by the colour that corresponds to its average grade and representative sublevels are shown in grey. It can be observed that the sublevel distance constraint was respected and internal waste was allowed by the program. High sublevel building cost clearly forced the program to choose as few sublevels as possible that cover access to valuable sections in the mine. Although sublevel generation was faster due to the equal distance enforcement, the internal waste option increased the solution time because of the increased number of possibilities.

**Conclusion**

The case studies demonstrate that the approach works well without violating any constraints. Also, the program parameters have been shown to converge. In other words, the parameters can be set empirically and increased until the profit does not improve. Further investigation of the approach involved comparison of the results to an exact method. To achieve this, the problem was formulated as a mixed integer linear program model and a similar case solved both with a linear program solver and the proposed heuristic clustering approach. The results show that the mine profits generated by both approaches were very similar but the heuristic approach reached this result much faster.

The advantages of the proposed approach are that it (1) follows the engineering practices, (2) quickly generates comparable results to optimal, and (3) if sublevels were previously developed in the mine, they can be provided manually owing to the modular structure of the approach. The comparison to the linear programming model produced promising results. The focus in future will be on incorporating alternative stope shapes into the model as suggested by Jooste and Malan (2020). Also, the research will be directed at generating a robust solution through the inclusion of grade uncertainty.

**Acknowledgements**

This research was conducted with financial support from the Natural Sciences and Engineering Research Council of Canada (Grant No. 488262-15).

**References**

Alford, C. 1995. Optimization in underground mine design. *Proceedings of APCOM XXV: Application of Computers and Operations Research in the Minerals Industries, *Brisbane, Australia, 9-14 July 1995. Australasian Institute of Mining and Metallurgy, Melbourne. pp. 213-218. [ Links ]

Alford, C., Brazil, M., and Lee, D.H. 2007. Optimisation in underground mining. *Handbook of Operations Research in Natural Resources. *Springer. pp. 561-577. [ Links ]

Ataee-poür, M. 1997. A new heuristic algorithm to optimise stope boundaries. *Proceedings of the Second Regional APCOM Symposium on Computer Applications in the Mineral Industry, *Moscow, 24-28 August 1997, Moscow State Mining University. pp. 511-515. [ Links ]

Ataee-Pour, M. 2004. Optimisation of stope limits using a heuristic approach. *Mining Technology, *vol. 113, no. 2. pp. 123-128. [ Links ]

Bai, X., Marcotte, D., and Simon, R. 2013. Underground stope optimization with network flow method. *Computers & Geosciences, *vol. 52. pp. 361-371. [ Links ]

Cheimanoff, N.M., Deliac, E.P., and Mallet, J.L. 1989. GEOCAD: an alternative CAD and artificial intelligence tool that helps moving from geological resources to mineable reserves. *Proceedings of APCOM XXI: Application of Computers and Operations Research in the Minerals Industries, *Las Vegas, NV, 27 February - 2 March, 1989, American Institute of Mining, Metallurgical, and Petroleum Engineers. pp. 471-478. [ Links ]

Copland, T. and Nehring, M. 2016. Integrated optimization of stope boundary selection and scheduling for sublevel stoping operations. *Journal of the Southern African Institute of Mining and Metallurgy, *vol. 116, no. 12. pp.1135-1142. [ Links ]

Deraisme, J., de Fououet, C., and Fraisse, H. 1984. Geostatistical orebody model for computer optimization of profits from different underground mining methods. *Proceedings of APCOM XVIII: Application of Computers and Operations Research in the Minerals Industries, *London, UK, 26-30 March, 1984. Institution of Mining and Metakkurgy. pp. 583-590. [ Links ]

Foroughi, S., Hamidi, J.K., Monjezi, M., and Nehring, M. 2019. The integrated optimization of underground stope layout designing and production scheduling incorporating a non-dominated sorting genetic algorithm (NSGA-II). *Resources Policy, *vol. 63. https://doi.org/10.1016/j.resourpol.2019.101408 [ Links ]

Hartman, H.L. and Mutmansky, J.M. 2002. Introductory Mining Engineering: Wiley. [ Links ]

Hou, J., Xu, C., Dowd, P.a., and Li, G. 2019. Integrated optimisation of stope boundary and access layout for underground mining operations. *Mining Technology, *vol. 128, no. 4. pp. 1-13. [ Links ]

Jalali, S.E. and Ataee-pour, M. 2004. A 2D dynamic programming algorithm to optimise stope boundaries. *Proceedings of the Thirteenth International Symposium on Mine Planning and Equipment Selection, *Wroclaw, Poland, 1-3 September 2004. CRC Press. pp. 45-52. [ Links ]

Johnson, T.B. and Sharp, W.R. 1971. A three-dimensional dynamic programming method for optimal ultimate open pit design. Vol. 7553, *Report of Investigations. *Bureau of Mines, US Dep. of the Interior. [ Links ]

Jooste, Y. and Malan, D. 2020. The need for improved layout design criteria for deep tabular stopes. *Journal of the Southern African Institute of Mining and Metallurgy, *vol. 120, no. 1. pp. 23-32. [ Links ]

Koch, Τ. 2006. Zimpl user guide: Zuse Institute Berlin. https://zimpl.zib.de/download/zimpl.pdf [ Links ]

Maseko, V. and Musingwini, C. 2019. An empirical long-term commodity price range for Mineral Reserve declarations to minimize impairments in gold and platinum mines. *Journal of the Southern African Institute of Mining and Metallurgy, *vol. 119, no. 3. pp. 229-242. [ Links ]

Nicholas, D.E. 1981. Method selection - A numerical approach. *Design and Operation of Caving and Sublevel Stoping Mines, *Chap. 4. Stewart, D. (ed.). SME-AIME, New York. pp 39-53. [ Links ]

Nikbin, v., Ataee-pour, Μ., Shahriar, Κ., Pourrahimian, Y., and MirHassani, S.A. 2019. Stope boundary optimization: A mathematical model and efficient heuristics. *Resources Policy, *vol. 62. pp. 515-526. [ Links ]

Nikbin, V., Ataee-Pour, Μ., Shahriar, Κ., and Pourrahimian, y. 2020. A 3D approximate hybrid algorithm for stope boundary optimization. *Computers & Operations Research, *vol. 115. https://doi.org/10.1016/j.cor.2018.05.012 [ Links ]

Osman, LH. and Christofides, Ν. 1994. Capacitated clustering problems by hybrid simulated annealing and Tabu search. *International Transactions in Operational Research, *vol. 1, no. 3. pp. 317-336. doi: 10.1111/1475-3995. d01-43 [ Links ]

Ovanic, J. and Young, D.S. 1995. Economic optimisation of stope geometry using separable programming with special branch and bound techniques. *Proceedings of the Third Canadian Conference on Computer Applications in the Mineral Industry, *Montreal, Canada, 22-25 Octobe, 1995. McGill University. pp. 129-135. [ Links ]

Remy, Ν., Boucher, Α., and Wu, J. 2009. Applied Geostatistics with SGeMS: A User's Guide. Cambridge University Press. [ Links ]

Sandanayake, D.S.S., Topal, E., and Asad, M.W.A. 2015. A heuristic approach to optimal design of an underground mine stope layout. *Applied Soft Computing, *vol. 30. pp. 595-603. [ Links ]

Sari, Y.A. and Kumral, Μ. 2019. A planning approach for polymetallic mines using a sublevel stoping technique with pillars and ultimate stope limits. *Engineering Optimization, *vol. 52, no. 6. pp. 932-944. doi: 10.1080/0305215x.2019.1624739 [ Links ]

Sari, Y.A. and Kumral, Μ. 2020. Sublevel stope layout planning through a greedy heuristic approach based on dynamic programming. *Journal of the Operational Research Society, *pp. 1-10. doi: 10.1080/01605682.2019.1700179 [ Links ]

Sauvageau, Μ. and Kumral, Μ. 2017. Kalman filtering-based approach for project valuation of an iron ore mining project through spot price and long-term commitment contracts. *Natural Resources Research, *vol. 26, no. 3. pp. 303-317. doi: 10.1007/s11053-017-9329-4 [ Links ]

Topal, E. and Sens, J. 2010. A new algorithm for stope boundary optimization. *Journal of Coal Science and Engineering *(China), vol. 16, no. 2. pp. 113-119. [ Links ]

Villalba Matamoros, M.E. and Kumral, Μ. 2019b. Underground mine planning: Stope layout optimisation under grade uncertainty using genetic algorithms. *International Journal of Mining, Reclamation and Environment, *vol. 33, no. 5. pp. 353-370. [ Links ]

Villalba Matamoros, M.E. and Kumral, Μ. 2019a. A value adding approach to hard-rock underground mining operations: Balancing orebody orientation and mining direction through meta-heuristic optimization. *Journal of Central South University, *vol. 26, no. 5. pp. 3126-3139. [ Links ]

Wang, H. and Webber, Τ. 2012. Practical semiautomatic stope design and cutoff grade calculation method. *Mining Engineering, *vol. 64. no. 99. pp. 85-91. [ Links ]

Zhao, X., Bennell, J.A., Bektas, Τ., and Dowsland, Κ. 2016. A comparative review of 3D container loading algorithms. *International Transactions in Operational Research, *vol. 23, no. 1-2. pp. 287-320. doi: 10.1111/itor.12094 [ Links ]

Zhu, Τ., Balakrishnan, J., and Cheng, C.H. 2018. Recent advances in dynamic facility layout research. *INFOR: Information Systems and Operational Research, *vol. 56, no. 4. pp. 428-456. [ Links ]

Zhu, X. and Lei, D.-Y. 2018. Bi-level hybrid local search approach for three-dimensional loading problem with balancing constraints. *Journal of Central South University, *vol. 25, no. 4. pp. 903-918. doi: 10.1007/s11771-018-3793-9 ♦ [ Links ]

**Correspondence**:

Y.A. San ]]>
Email: yuksel.sari@mail.mcgill.ca

Received: 2 Jun. 2020

Revised: 14 Jun. 2020

Accepted: 20 Jan.2021

Published: March 2021