Engineering Optimization

ISSN: (Print) (Online) Journal homepage: https://www.tandfonline.com/loi/geno20

Effective solution of the long-term open pit

production planning problem using block
clustering

Reza Lotfian , Javad Gholamnejad & Yousef Mirzaeian Lardkeyvan

To cite this article: Reza Lotfian , Javad Gholamnejad & Yousef Mirzaeian Lardkeyvan (2020):
Effective solution of the long-term open pit production planning problem using block clustering,
Engineering Optimization, DOI: 10.1080/0305215X.2020.1771703

To link to this article: https://doi.org/10.1080/0305215X.2020.1771703

Published online: 07 Aug 2020.

Submit your article to this journal

Article views: 5

View related articles

View Crossmark data

Full Terms & Conditions of access and use can be found at

https://www.tandfonline.com/action/journalInformation?journalCode=geno20
ENGINEERING OPTIMIZATION
https://doi.org/10.1080/0305215X.2020.1771703

Effective solution of the long-term open pit production planning

problem using block clustering
Reza Lotfian, Javad Gholamnejad and Yousef Mirzaeian Lardkeyvan

Mining and Metallurgical Engineering, Yazd University, Yazd, Iran

ABSTRACT ARTICLE HISTORY

Production planning of an open pit mine is defined as determining the Received 14 September 2019
extraction sequence of blocks in the final pit limit with the goal of maxi- Accepted 17 May 2020
mizing Net Present Value (NPV) while considering the constraints imposed KEYWORDS
on the mining system. Mathematical programming models are used to Long-term production
optimize long-term open pit production planning problems. In real mines planning; open pit mine;
with several thousands of blocks and several scheduling periods, however, mathematical programming;
these methods lose their efficiency by increasing the number of decision block clustering
variables, and it is almost impossible to solve and obtain a solution in a rea-
sonable time using current solvers. In this paper, a new method based on
block clustering was proposed to reduce the model size. In the first step,
blocks are aggregated in clusters called “mother clusters” using the concept
of clustering as a mathematical programming problem and solved using
Genetic Algorithm (GA). In the second step, the mother clusters are con-
verted to the mining clusters considering the practicality issues. The results
of implementing the proposed method on a real deposit showed that block
clustering could decrease the number of binary variables from 476,100 to
5,652. This led to a significant reduction in solution time from 372 hours to
23.4 seconds compared to the original block model. The clustering strategy
was also able to realize more than 82% of the NPV obtained from scheduling
original blocks.

1. Introduction
Mining can be considered a complex industry with high investment costs and revenue potentials.
These huge costs require optimal planning for the extraction of ores to maximize return on invest-
ment, which ensures large revenues and the profitability of mines. Open pit mining is a well-known
and commonly used method for extracting minerals from the ground. Long-term production plan-
ning (LTPP) has been one of the most important issues in both underground and open pit mine
design and production planning. In open pit mining, the extraction of shallow and high-grade ores,
and selection of deeper and lower grade deposits, have doubled the importance of optimal production
planning. Mathematical programming models are suitable for optimizing the LTPP problem because
of their ability to model various mining constraints, such as multiple ore processing, multiple material
stockpiling, and blending constraints. Since 1960, such models have been used with deterministic and
uncertainty-based approaches. Osanloo, Gholamnejad, and Karimi (2008) investigated these models
and algorithms with deterministic and uncertainty-based approaches. The formulation of this prob-
lem is more conventional and well known in the form of integer programming (IP) or mixed integer
programming (MIP) models. In real mines with thousands of blocks and several scheduling periods,

CONTACT Javad Gholamnejad [email protected]

© 2020 Informa UK Limited, trading as Taylor & Francis Group
2 R. LOTFIAN ET AL.

however, these methods are no longer efficient owing to the increasing number of binary decision
variables and consequently an increase in the problem size. Therefore, it is almost impossible to obtain
an optimal solution in a reasonable time using current solvers.
For this reason, various approaches have been developed to solve this problem, which can be
divided into two general categories. The first group comprises methods that try to solve the prob-
lem by changing the structure of the mathematical model and are themselves classified into three
categories of decomposition, relaxation and metaheuristic techniques. The second group of methods
tries to aggregate blocks and convert them into larger extraction units (clusters) to reduce the number
of variables and constraints. The clustering technique was used in this article. In most of the previous
studies, clusters were built on each bench for simplicity so that the thickness or depth of each cluster
had to be as high as one bench. This approach does not necessarily achieve the highest net present
value (NPV). In this article, new mining units, called ‘mining clusters’, are produced in such a way
that, in addition to having the highest intracluster similarity and intercluster dissimilarity in terms of
distance and grade of blocks, the clusters are not limited to a bench and their thickness can be greater
than the height of one bench. In addition, the safe slope angle is satisfied in each cluster.
In what follows, the literature on the LTPP problem is reviewed based on aggregation of blocks.
Then, a two-step algorithm is proposed for aggregating blocks and creating mining clusters. In the
next section, the algorithm is implemented on a hypothetical two-dimensional model, and then on a
real deposit. Finally, conclusions and some suggestions for future research are presented.

2. Literature review on block aggregation approaches used to solve the LTPP

problem
The difficulty of solving the LTPP problem depends to a large extent on the number of blocks in the
ultimate pit limit, which should be planned. Therefore, if planners can group or aggregate the blocks
into larger extraction units, the number of decision variables and constraints in the optimization
model will decrease and then the model can be solved by conventional mathematical programming
methods. Various mathematical models of production planning have been suggested based on group-
ing and aggregation of blocks. Gershon (1983) presents a column model whereby, instead of using
blocks as extraction units, the blocks are aggregated in order to divide the deposit into columns or
pillars, which develop from the bottom to the top of the deposit. In his model, continuous instead
of integer variables are used in production planning. Although some integer variables are needed in
other constraints, such as for differentiating between ore and waste, the number of integer variables
is far lower than in models where original blocks are used for scheduling. With partial modification
of the column model, Gershon (1983) suggested a simpler and more generalized model, called the
directional model, in which strips of ore are defined as extraction units in a given direction. The main
criterion for selecting such a direction is ore homogeneity. Later on, Busnach, Mehrez, and Sinuany-
Stern (1985) developed an IP model for production planning of a phosphate mine to maximize NPV.
They aggregated the extraction blocks and turned them into horizontal mining layers. Klingman and
Phillips (1988) used IP to optimize the extraction of phosphate mines. They formed extraction units,
taking into account the layering characteristic of phosphates and the capabilities of mining equipment
in the mine. In their research, each mining unit has a thickness equal to that of one stratum and an
area as large as a region that can be extracted without the need to relocate the dragline.
Smith (1999) used a goal programming model and MIP to optimize short-term production
scheduling in a phosphor formation in Idaho State, USA. Each mining block in this study has a
height of 20 feet, a length of 100 feet along the seam, and a width equal to that of the seam. This
model uses two types of continuous decision variables and two types of binary decision variables.
To reduce the number of binary variables, the study used strategies such as reducing the number of
blocks per deposit via combining two or more blocks and converting them into an aggregated unit
called a ‘super-block’, removing air blocks from the model, adding blocks containing insignificant
tonnage to the next block in the lift, using the concepts of early start and late start extraction of a
 ENGINEERING OPTIMIZATION 3

given block, and reducing the number of periods by increasing the length of each period. Using the
above strategies, the number of binary variables in his model was reduced from 4452 to 198, which
can be solved in a reasonable time with relatively low computational resources.
Gershon and Murphy (1989) used assay results taken at regular intervals of depth using a core
drawn from a single drill hole, and considered the depth intervals where the ore body was entirely
extracted as ore or waste as mining units or segments. They used dynamic programming to make
decisions about the extraction of these segments. Samanta, Bhattacherjee, and Ganguli (2005) used a
layering approach to apply the genetic algorithm (GA) for production planning with the goal of grade
control during the extraction period.
Ramazan, Dagdelen, and Johnson (2005) and Ramazan (2007) introduced a new approach for the
LTPP problem with aggregation of blocks based on the concept of fundamental trees in the form of
fundamental tree algorithms. They proposed a linear programming model to determine the aggre-
gated units. After finding the fundamental trees, the LTPP problem is formulated as an MIP model,
where each fundamental tree is treated as a mining block containing an ore tonnage, metal content,
grade and quality parameters, and each fundamental tree is assigned a binary variable for each pro-
duction period. Ramazan, Dagdelen, and Johnson (2005) claimed that blocks were aggregated in the
form of fundamental trees to avoid the reduction of optimality and resolution of the model. Their
results obtained by solving the MIP model based on fundamental trees also showed a significant
improvement in solution time relative to the block-based models. Moreover, fundamental tree-based
scheduling created a higher NPV than the common commercial software.
Tabesh and Askari-Nasab (2011) and Askari-Nasab, Tabesh, and Badiozamani (2010) used the
hierarchical agglomerative clustering technique for aggregating blocks in their LTPP model. Their
algorithm aggregated blocks on each bench based on the similarity index, which is defined based on
rock type, ore grade and distances between blocks. The aggregated blocks, called ‘mining cuts’, were
then used in an MIP model to obtain the long-term plans. They implemented their algorithm on a
block model containing 19,492 blocks that had to be extracted during a 20 year time horizon. Ben-
Awuah, Askari-Nasab, and Awuah-Offei (2012), Badiozamani and Askari-Nasab (2012) and Ben-
Awuah and Askari-Nasab (2012) used the same approach to optimize production planning in Oil
Sands Mines. Askari-Nasab, Awuah-Offei, and Eivazy (2010) used two MIP models for LTPP at two
different levels. To reduce the number of binary variables, they aggregated the blocks on each bench
based on the fuzzy clustering approach and the fuzzy C-means technique, and converted the blocks
into larger units called extraction cuts.
Goodfellow and Dimitrakopoulos (2016) suggested a two-stage stochastic optimization model for
production planning in open pit mining complexes. To investigate grade and geological uncertainty
in their study, various realizations of the geological and grade model were created and the blocks in
the pit were clustered using the k-means++ algorithm based on grade criterion. Koushavand, Askari-
Nasab, and Deutsch (2014) presented an MIP model for LTPP in the presence of grade uncertainty
with the goal of maximizing NPV and minimizing uncertainty costs. They defined two separate vari-
ables in their model to control processing and extraction in order to achieve more freedom of the
optimizer to converge more quickly to higher NPV. In this way, the number of decision variables for
this model for a pit will be 3 * (T * N) + 2(T * R), where N is the number of blocks, T is the life of
the mine, R indicates the number of deposit realizations, and (T * N) are binary variables. In their
research, the blocks are aggregated on each bench based on the grade characteristic and using the C-
®
mean clustering function in MATLAB . Hoerger et al. (1999) studied the production planning process
at Newmont Mining Complex in Nevada, USA, suggesting a multi-period MIP model to optimize a
large-scale, complex operation.
Because of the variation in chemical characteristics of the deposit in the Twin Creeks mining
complex in Nevada, it was divided into several different metallurgical types. In addition, each met-
allurgical type was classified into several grade intervals. In this way, all the materials in different
pits were clustered based on metallurgical and grade similarities in the k types of increments, and
4 R. LOTFIAN ET AL.

then these clustered materials were cut into sequences. Each sequence nominally contains 1 year
extraction materials. Therefore, the planning units in this study are sequences that use in-house soft-
ware developed by the company to compute and format grade-tonnage reports for each sequence as
inputs of the MIP model. Kawahata, Schumacher, and Criss (2016) also studied the strategic plan-
ning process extraction in the same mineral complex and presented an MIP model similar to the
one presented by Hoerger et al. (1999). The process of clustering the material inside the pit into
increments is similar to that of the above-mentioned study, but the shapes of the mining sequences
are different. Kawahata, Schumacher, and Criss (2016) considered each mining bench as an extrac-
tion sequence. Kawahata, Schumacher, and Hufford (2013) used Optimizer and Minemax software to
optimize production planning in the Twin Creeks mining complex and compared the performance of
each software in generated production plans. Optimizer software considers two approaches to reduce
the model size. First, benches extracted in the base-plan schedule (which is entered in Optimizer as
a starting point for the optimization) over an identical period are aggregated in the form of a large
bench. In their study, the smallest planning units are benches, which are clustered with the use of
aggregation. The second technique used by Optimizer to reduce the model size is to limit the num-
ber of time-period options to optimize extraction benches. The ‘Sliding Window’ option in Minemax
optimization software formulates and solves a set of subproblems taking into account the limited time
frames in each subproblem of the model and continues iteratively until it covers the life of the mine.
The ‘Forward’ option is similar to the Sliding Window option, except that it formulates and solves the
model for each time period and locks the solution. Both the Sliding Window and Forward methods
reduce the solution time significantly. Goodwin et al. (2006) showed that the production planning
problem could be formulated using a control engineering language in the framework of receding
horizon control (RHC). The main idea of optimizing the RHC in control issues is the same as Win-
dow Sliding in Minemax software. The RHC approach, however, uses non-uniform discretization
of time, which is associated with changes in the mining capacity constraint. Chanda and Ricciar-
done (2002) used Minemax planning software for LTPP of a multi-element open pit mine. Given
that Minemax uses the MIP model with the branch-and-bound algorithm as the solution method,
model size and solving time still present challenges for large-scale mines. As pointed out above,
the Sliding Window and Forward strategies can be considered as two ways of reducing the model
size. Another approach in Minemax to reduce the model size is the ‘Grade Re-block’ option through
importing block model files into the software. Accordingly, when the block model is imported into
the software, the blocks fall into clusters, called ‘bin’, based on their characteristics (such as grade). In
Minemax, it is assumed that every pit is extracted bench by bench, and there is no need to consider
the sequence preferences or working slope angles of the wall or the maximum vertical development
constraint.
The NPV Scheduler software is also widely used to optimize the LTPP of open pit mines. This
software divides the blocks on each bench within the optimal pit into units called atoms, and then
solves the problem, but it does not provide details on how to aggregate the blocks and convert them
to atoms. To solve the model of the LTPP problem, Jélvez et al. (2016) suggested an approach based
on block aggregation that aggregated blocks based on the distance between them. Smith and Wicks
(2014) studied the application of an MIP model with the goal of maximizing the recovered metal for
medium-term production scheduling. The primary block model in their study included blocks with
dimensions of 15 × 15 × 12 m. They defined mining blocks with dimensions of 100 × 100 × 12 m
on each bench as a bulk of material extracted by the shovel during a scheduling period (1 month).
Accordingly, the primary blocks were then clustered in packs, called bins, based on the mineralogical
criteria, metallurgical performance and metal grade, and then a continuous variable was defined for
deciding on the amount of tonnage sent from each bin.
Overall, the inefficacy of non-layered natural deposits, inhomogeneity of mining clusters and cre-
ating clusters on mining benches are the most important disadvantages of the above research. In most
previous studies in the field of block aggregation, the thickness of each cluster was set to the height of
each bench. This does not allow the growth of a cluster in depth, and hence does not necessarily lead
 ENGINEERING OPTIMIZATION 5

to an optimal solution. In this article, therefore, a method is proposed to address the flaws of previous
methods.

2.1. Other emerging areas in the field of clustering methodology

The widespread use of clustering in different sciences has led to new methods and more developed and
efficient algorithms in the methodology of this concept. One of these methods is clustering ensemble,
which attempts to create a robust partitioning of the data by combining the different partitionings
generated by the basic (primary) clustering methods (Strehl and Ghosh 2002; Fred and Jain 2005;
Kuncheva and Hadjitodorov 2004). Clustering ensemble can provide better solutions in terms of
robustness, novelty, stability and flexibility than basic methods (Topchy, Jain, and Punch 2003; Ayad
and Kamel 2008; Fred and Lourenço 2008). In summary, clustering ensemble consists of two main
steps (Strehl and Ghosh 2002): (1) generating different clusterings as primary clustering results based
on different methods, called the diversity creation step (Yousefnezhad, Huang, and Zhang 2018); and
(2) aggregating the obtained primary partitionings to obtain the final partitioning, which is done by
a consensus function (Huang, Wang, and Lai 2018).
Along with the advantages of clustering ensemble, this method also faces some challenges, one
of the most important of which is how clusters are combined. In clustering ensemble, the results of
primary clusterings may be inconsistent. For example, two different methods have the same sam-
ple in two different clusters. As a result, the idea is not to treat the views of different clusterings
alike, but rather to weigh each one. In most previous studies, each partition is given an equal weight
in the combination process and all clusters in each partition with equal weights participate in the
final ensemble (Fred and Jain 2006). In this regard, Rashidi et al. (2019) proposed a new clustering
ensemble approach using a weighting strategy. A lack of quality in the base clusterings makes their
consensus clustering weak. Therefore, a subset of high-quality basic clusters should be selected to
create a consensus. The quality of the clusters is determined by the assessment metric of the cluster.
Nazari et al. (2019) proposed a new clustering ensemble framework based on cluster-level weight-
ing. Inspired by bagging and boosting algorithms in classification, Parvin et al. (2013) proposed and
studied the non-weighing and weighing-based sampling approaches for clustering. With the same
approach, Minaei-Bidgoli et al. (2014) also provided non-adaptive and adaptive resampling schemes
for the integration of multiple independent and dependent clusterings.
Furthermore, traditional clustering algorithms consider equally all the features of the data set
under consideration. It is obvious that some features have more information than others in a data
set. Thus, it is highly probable that some features should have lower degrees of importance during a
clustering or a classification algorithm owing to their lower information, higher variances, etc. Parvin
and Minaei-Bidgoli (2013) considered the participation of features in the clustering process, as they
weighted and used the weighted locally adaptive clustering algorithm instead of the locally adaptive
clustering algorithm. Parvin and Minaei-Bidgoli (2015) developed the fuzzy weighted locally adap-
tive clustering method by presenting the fuzzy weighted clustering algorithm. Abbasi et al. (2019)
defined the edited normalized mutual information evaluation measure rather than the normalized
mutual information (NMI) criterion used by most previous studies, and then used three types of con-
sensus function to create the final clusters. Alizadeh, Minaei-Bidgoli, and Parvin (2014) applied the
new Alizadeh–Parvin–Moshki–Minaei criterion instead of the well-known NMI measure. Bagherinia
et al. (2019) proposed a novel fuzzy clustering ensemble framework based on a new fuzzy diversity
measure and a fuzzy quality measure to find the base clusterings with the best performance. In their
research, diversity and quality are considered simultaneously and a final clustering is obtained by
applying two types of consensus function. In line with the second step of clustering ensemble men-
tioned in the first paragraph of this subsection, Mojarad, Nejatian, et al. (2019) provided an ensemble
aggregator, called the robust clustering ensemble based on sampling and cluster clustering. Mojarad,
Parvin, et al. (2019) also proposed a clustering ensemble method with a new aggregation function.
6 R. LOTFIAN ET AL.

The rest of this article is organized as follows. In Section 3, the concept of clustering and the
definition of the k-means clustering algorithm will be described as a mathematical programming
model together with its solution by the GA. Section 4 explains the block aggregation algorithm used
in this article. An IP formulation to determine the sequence of extraction for mining clusters is pre-
sented in Section 5. Implementations on a two-dimensional model and on a real mine are given in
Section 6. A discussion is presented in Section 7, and conclusions and recommendations for future
work are presented in Section 8.

3. Data clustering
Clustering is one of the data mining algorithms. It is a general term for a number of algorithms that
group objects and put similar objects together in separate segments called clusters. Clustering in N-
dimensional Euclidean space is the process of partitioning a given set of n points into a number, say
K, of groups (or clusters) based on some similarity/dissimilarity metrics (Maulik and Bandyopad-
hyay 2000). Let the set of n points {x1 , x2 , . . . , xn } be represented by the set S and the K clusters be
represented by C1 , C2 , . . . , CK . Then:

Ci = ∅ for i = 1, . . . , K (1)
Ci ∩ Cj = ∅ for i = 1, . . . , K, j = 1, . . . , K, i = j (2)
K
∪ Ci = S (3)
i=1

The most common goal of clustering techniques is to maximize intracluster similarity and intercluster
dissimilarity. Clustering techniques are categorized into two main groups of partitional and hierar-
chical clustering. Assuming n as the input data and k as the number of final clusters, members in
partitional clustering are included in the k cluster. One of the most important techniques in parti-
tional clustering is k-means. This technique is one of the popular data clustering techniques widely
used for various applications (Jain 2010; Sinha and Jana 2016; Fahad et al. 2014). The steps of the
k-means algorithm can be described as follows (Maulik and Bandyopadhyay 2000):

• Step 1: Choose K initial cluster centres l1 , l2 , . . . , lK randomly from the n points {x1 , x2 , . . . , xn }.
• Step 2: Assign point xi , i = 1, 2, . . . , n to cluster Cj , j ∈ {1, 2, . . . , K} if xi − lj < xi − lp , p =
1, 2, . . . , K and j = p.
• Step 3: Compute new cluster centres l1∗ , l2∗ , . . . , lK∗ as follows:

1 
lj∗ = xi j = 1, 2, . . . , K
ni
xi ∈Cj

where ni is the number of elements belonging to cluster Cj .

• Step 4: If lj∗ = lj , j = 1, 2, . . . , K then terminate. Otherwise continue from Step 2.

This problem can be expressed as a mathematical optimization model:


Min Z = aij d(xi , lj ) (4)
Cj xi ∈Cj

subject to

1 xi ∈ Cj
aij
0 xi ∈
/ Cj
 ENGINEERING OPTIMIZATION 7

where j is the cluster index (j = 1,2, . . . K); i is the data index (i = 1,2, . . . n); xi is the ith data; Cj is
the jth cluster; lj is the centre of the jth cluster; d(xi ,lj ) is the distance between the ith data and centre
of the jth cluster; Z is the value of the objective function of the problem and represents the sum of the
distances of all members of each cluster from its centre; and aij is a binary variable that is equal to 1
if the ith data is a member of jth cluster, otherwise it is zero.
Solving this problem directly requires a lot of resources and time, and hence it has been proved to
be an NP-hard problem (Gonzalez 1982). Therefore, metaheuristic algorithms can be used to solve
this problem. In this study, GA is used to solve the model. The searching capability of GA has been
used in this article to search for appropriate cluster centre locations in the feature space in such a
way as to optimize a similarity metric of the resulting clusters, thereby suitably grouping the set
of n unlabelled points. The adopted clustering metric is the sum of the Euclidean distances of the
points from their respective cluster centres. Mathematically, the clustering metric M for the K clusters
C1 , C2 , . . . , CK is given by:

K 

M(C1 , C2 , . . . , CK ) = xi − lj (5)
j=1 xi ∈Cj

The task of the GA is to search for appropriate cluster centre locations in such a way as to minimize
the clustering metric M. The basic steps of the GA, which are also followed in this article, are shown
as follows (Davis 1991):

Begin
1. t = 0
2. Initialize population p(t)
3. Compute fitness p(t)
4. t = t+1
5. If termination criterion achieved go to Step 10
6. Select p(t) from p(t−1)
7. Crossover p(t)
8. Mutate p(t)
9. Go to Step 3
10. Output best and stop
End

4. Block aggregation algorithm

The block aggregation approach used in this article can be expressed in two steps. In the first step, the
blocks are aggregated in clusters, called ‘mother clusters’, using the concept of clustering. Clustering
is expressed as a mathematical programming model and solved by the GA. In the second step, the
‘mining clusters’ are created by applying an IP approach on the mother clusters. The mining clus-
ters can be extracted independently without violating slope constraints. Details of these steps are as
follows.

4.1. Step 1
In this step, blocks located in the final pit limit are aggregated in large clusters, called mother clusters,
based on the closeness of the blocks and the grade characteristic for each block. Clustering is first per-
formed by defining the clustering problem as a mathematical programming problem solved through
the GA. Input data for this model are all the blocks within the final pit limit, which are clustered
according to the distance between the blocks and their geological and metallurgical characteristics,
regardless of the overall pit slope angle.
8 R. LOTFIAN ET AL.

4.2. Step 2
After the first step, blocks located in the final pit limit are classified according to predefined charac-
teristics in the k cluster, called mother clusters. These clusters, however, are not practical and cannot
participate as planning units in the production planning process as they do not consider the slope
angle and also because of the large number of blocks in each cluster. For this purpose, a mathematical
programming model is defined and executed iteratively on each mother cluster in the second stage of
the algorithm. The purpose of this algorithm is to convert non-practical clusters into practical ones
by considering the slope angles therein. To this end, the largest possible pit algorithm is defined as
follows and applied to each cluster:
R

MaxE = yr (6)
r=1

subject to:
A

9yr − ya ≤ 0 ∀r (7)
a=1

where R is the total number of blocks in each mother cluster; A is a set of nine blocks overlying block
r; and yr is a binary variable equal to 1 if block r is inside the boundary of the largest possible pit,
otherwise it is zero.
The algorithm works in such a way that the blocks inside the largest possible pit are specified and
then removed in each run. This process is repeated as long as all blocks within the mother cluster are
assigned to a mining cluster. Therefore, all blocks in each mother cluster are classified into groups
called mining clusters, which independently respect the slope requirements of pit walls. Now, these
clusters can be used as planning units of a mathematical programming model and a decision should
be made for their extraction time.

5. IP formulation for LTPP

As noted earlier, various mathematical models have been used to optimize the problem of LTPP in
open pit mines. To determine the extraction sequence of mining clusters, the model presented by
Ramazan and Dimitrakopoulos (2004) is used in this study, with some modifications.
T 
 n
Maximize cit xit (8)
t=1 i=1
n

oi xit ≥ PCmin ∀t = 1, . . . , T (9)
i=1
n

oi xit ≤ PCmax ∀t = 1, . . . , T (10)
i=1
n

(oi + wi )xit ≥ MCmin ∀t = 1, . . . , T (11)
i=1
n

(oi + wi )xit ≤ MCmax ∀t = 1, . . . , T (12)
i=1
n t

i=1 oi gi xi ≥ Gmin ∀t = 1, . . . , T (13)
n t
i=1 oi xi
 ENGINEERING OPTIMIZATION 9

Figure 1. Final pit of the hypothesized two-dimensional model.

n t

i=1 oi gi xi ≤ Gmax ∀t = 1, . . . , T (14)
n t
i=1 oi xi
t

xkt − xlr ≤ 0 ∀k = {1, 2, . . . , n}, t = {1, 2, . . . , T}, l ∈ Ck (15)
r=1
T

xit = 1 ∀i = 1, . . . , n (16)
t=1

where T is the maximum number of scheduling periods; n is the total number of mining clusters to
be scheduled; ci t is the NPV to be generated by extracting cluster i in period t; xi t is a binary variable,
equal to 1 if cluster i is to be mined in period t, 0 otherwise; oi is the ore tonnage in cluster i; wi is
the tonnage of waste material in cluster i; PCmin is the lower bound on processing capacity of ore in
period t; PCmax is the upper bound on processing capacity of ore in period t; gi is the average grade
of mining cluster i; Gmin is the lower bound on acceptable average head grade in period t; Gmax is
the upper bound on acceptable average head grade in period t; MCmin is the lower bound on mining
capacity in period t; MCmax is the upper bound on mining capacity in period t; k is the index of a
cluster considered for excavation in period t; and Ck is a set of immediate predecessor mining clusters
that must be extracted before extraction of the mining cluster k.
Equation (8) is the objective function that aims to maximize the overall NPV of the project.
According to these constraints, the total tonnage of ore processed in each period should be greater
than a lower bound (PCmin ) and less than an upper bound (PCmax ), as expressed in Equations (9) and
(10). In the same way, the total tonnage of rock (waste and ore) to be mined in each period should be
greater than a lower bound (MCmin ) and less than an upper bound (MCmax ), as expressed in Equa-
tions (11) and (12). Equations (13) and (14) state that the average grade of the mineralized material
sent to the mill has to be greater than a lower bound (Gmin ) and less than an upper bound (Gmax )
in each period. Constraints (15) ensure that cluster k is extracted only if its predecessor clusters are
extracted in the same time period or prior to k. Constraints (16) require that each cluster can be
extracted no more than once.

6. Implementation and results

6.1. Implementing the algorithm on a two-dimensional model
To investigate the efficiency of the proposed clustering algorithm, the proposed method is first imple-
mented on a two-dimensional model of a hypothetical final pit. The results are then compared with
the other two models, first the one where the original blocks in the block model are used as min-
ing units for the production planning problem (block-based model), and secondly, the model whose
clusters are formed within each bench (on-bench cluster model); in other words, the depth of clusters
is as high as one bench.
For this purpose, a vertical cross-section of a final pit is used containing 59 blocks (Figure 1), in
which the numbers in each block indicate its grade percentage.
10 R. LOTFIAN ET AL.

Figure 2. Mining clusters created using the suggested algorithm. Step 1: Formation of mother clusters, where the numbers shown
in the blocks show the number of mother clusters. Step 2: Creation of mining clusters from mother clusters, where the numbers
shown in the blocks show the number of mining clusters.

Figure 3. Production planning results of mining clusters.

Usually, a certain economic value is assigned to each block based on the recovery of the production
process, the price of the product (metal or concentrate) and the operational costs. It is worth noting
that the value of each aggregation unit is equal to the sum of the values of the blocks it contains.

6.1.1. Proposed mining clustering strategy

In this strategy, mother clusters are formed in the first step. In the second step, the mining clusters
are created from the mother clusters. Thus, blocks are aggregated in such a way that each cluster
meets operational constraints while respecting the homogeneity and proximity of blocks. Because
the aggregated units are used as mining units in the production planning model, they should meet
the slope requirements of pit walls, and have the least variability in grade. Mother and mining clusters
created with this strategy are depicted with different numbers in Figure 2. In this model, clusters can
also be developed in depth. This strategy is used to implement the production planning for a 4 year life
of a mine. Maximum and minimum production capacities are 18 and 16 blocks per year, respectively.
The calculated NPV of the project is 1066 monetary units. The production planning results are shown
in Figure 3. The numbers in the blocks indicate their extraction times.

6.1.2. On-bench clustering strategy

In this approach, block aggregation is done only on each bench based on the grade characteristic
and the distance between blocks in the two-dimensional block model shown in Figure 1, and the
results are shown in Figure 4. Then, production planning is carried out for a 4 year life of a mine
assuming maximum and minimum production capacities of 18 and 16 blocks per year, respectively.
The calculated NPV of the project is 865 monetary units. The production planning results are shown
in Figure 5. The NPV resulting from this strategy is 18.8% less than that proposed in Section 6.1.1.

6.1.3. Production planning on the block-based model

To compare the results of the two previous methods with production planning based on the original
blocks, the block scheduling is used for 59 blocks in the final pit for a 4 year life of a mine with the
same constraints of production capacity and discount rate. The NPV was 1339 monetary units. As
expected, the NPV obtained from block scheduling is higher than that derived from scheduling on
 ENGINEERING OPTIMIZATION 11

Figure 4. Aggregation of blocks on each bench based on the distance and grade characteristics, where the numbers in the blocks
show the number of mining clusters.

Figure 5. Production planning results of aggregating blocks on each bench.

Figure 6. Production planning results based on the original blocks.

the two clustering-based models. According to the block scheduling results (Figure 6), blocks that
must be extracted each year are shown with different numbers.

6.2. Implementing the proposed algorithm on a real mine

To test the proposed aggregation algorithm, the model is applied on a copper open pit containing
23,296 blocks with dimensions of 5 × 5 × 6 m. Other model specifications are shown in Table 1. This
algorithm is coded in MATLAB release 2016a and GAMS 24.7.4 software, and is run on a system
®
with the specifications Intel CoreTM i3–6100 CPU, 3.7 GHz and 4 GB of RAM. After creating the
mining clusters, the IP model suggested in Ramazan and Dimitrakopoulos (2004) is used for the
production planning for a 4 year life of the mine. To construct the mother clusters, five characteristics,
comprising the X, Y and Z coordinates, as well as the Cu and CuO grades of blocks, are involved
in the clustering process and the blocks are grouped into seven mother clusters. The parameter K
(number of clusters) must be determined in such a way that high-grade blocks do not lose their values
in combination with low-grade blocks. In addition, by increasing the amount of K, the number of
mother clusters and subsequently that of mining clusters will increase and the planning process will
be more closely aligned with the block-based model. In this case study, it was found that the parameter
k = 7 can provide a good separation of the blocks. The clustering scheme resulting from the GA
is shown in Figure 7. Figure 8 shows the convergence of the GA, indicating the best value of the
cost function (Equation 5) in different iterations. Subsequently, by iterative running of the greatest
possible pit algorithm, the mining clusters are derived from the mother clusters and independently
respect the slope of pit walls constraint. In this way, all 23,296 blocks in the final pit limit are clustered
in 1413 mining clusters. Figure 9 shows the amounts of ore production and waste removal in each
period. Production planning for the same deposit is performed on the block-based model with GAMS
software. Table 2 summarizes the results of these two strategies.
12 R. LOTFIAN ET AL.

Figure 7. Results of clustering by the genetic algorithm for plan view 1146.

Figure 8. Convergence of the genetic algorithm.

7. Discussion
Comparison of the results of solving a hypothetical two-dimensional model shows that the highest
NPV of this project is equal to 1339, resulting from solving the problem on the original block-based
model. However, it is obvious that in large deposits, such as the case investigated in this study, the
solution time increases significantly in this method so that it is, in most cases, almost impossible to
obtain the solution in a reasonable time. In addition, in the case of obtaining the solution in many
real mines, determining the extraction sequence of the blocks by this method leads to the scattered
planning of blocks in the pit and it is hardly possible to relocate mining equipment repetitively in a
planning period. Therefore, operational production planning requires block aggregation. However,
 ENGINEERING OPTIMIZATION 13

Figure 9. Amounts of ore production and waste removal per year.

Table 1. Model specifications for a real mine.

Tonnage of rock 8,386,560
Tonnage of ore 2,733,480
Average grade of Cu (in ore) 0.51%
Maximum mining capacity (tons) 2,500,000
Minimum mining capacity (tons) 500,000
Maximum processing capacity (tons) 1,000,000
Minimum processing capacity (tons) 400,000
Maximum average grade of Cu (mill) 1.2%
Minimum average grade of Cu (mill) 0.15%

Table 2. Summary of the model solution using mining clustering and original block
model strategies by GAMS software.

Original
Mining block
clustering model
strategy strategy
Model type IP IP
Solver CPLEX CPLEX
Number of binary variables 5,652 476,100
Time solution 23.4 s 372 h
Number of iterations 10,075 1,584,423
Discount rate 10% 10%
Net present value 99,479,230,763 121,136,695,745
Relative gap 3% 0. 3%

the NPV generated by the original block-based scheduling can be considered as a good basis for
determining the quality level of planning resulting from aggregated units. Thus, the similarity of NPV
obtained from production planning based on aggregated units to that derived from production plan-
ning based on the original blocks indicates the efficiency of the proposed model. The NPV obtained
from the mining clustering strategy is at a distance of about 20% of the maximum possible NPV of
the project, while the one obtained from the on-bench clustering strategy is at a distance of 35.4%.
Comparing the NPVs of the two clustering strategies shows that access to high-grade blocks in the
early years has become possible owing to the form of clusters created by the proposed algorithm that
can grow in depth; besides, a positive cash flow can be achieved in the early years of the mine’s life.
By examining the amount of ore extraction in different years of the mine’s life, it is observed that the
mine will not produce any ore in the first year of planning in the on-bench block aggregation strategy.
14 R. LOTFIAN ET AL.

For real mines, this can lead to failure to achieve the economic goals of the project owing to a lack of
positive cash flow in the early years. This issue becomes important when the burden height is high
and the stable slope of the final pit wall is low owing to geomechanical considerations. In this case,
mine owners who apply the on-bench clustering strategy must incur the heavy costs of removing a
large amount of waste in the early years.
Using the mining clusters strategy in the real three-dimensional model significantly reduces the
number of binary variables (from 476,100 to 5652), resulting in a significant reduction in solution
time compared to the original block model (from 372 h to 23.4 s). It also resolves the issue of the
insufficiency of mathematical programming models in solving the LTPP problem. In addition, the
NPV resulting from the aggregation model is within a distance of 17.8% of the planning conducted
based on the original block model, which is acceptable and confirms the quality of the proposed
method and its practicality in real mines.

8. Conclusions and recommendations for further research

Optimization of the LTPP problem in real open pit mines using mathematical programming models
such as IP becomes a very large and intractable problem owing to the large number of blocks, followed
by the large number of binary variables. Block aggregation approaches that try to convert blocks into
larger units are used to reduce the size of the problem. In this article, mining clusters were created in
two steps. In the first step, the blocks were grouped into main clusters, called mother clusters, using
the GA based on the distance between the blocks and their geological specifications (e.g. grade). In the
second step, the greatest possible pit algorithm was run iteratively on each mother cluster to create
mining clusters that independently consider the slope requirements of pit walls. Then, one of the
commonly used IP models was used to solve the problem of production planning via aggregated
units for a real mine. In this way, 23,296 blocks in the final pit limit of a copper open pit mine were
aggregated in 1413 mining clusters. The results showed a significant reduction in the solution time
of the problem compared to that by the original block model (from 372 h to 23.4 s). The mining
clusters strategy was also able to realize more than 82% of the highest possible NPV (obtained from
original block-based scheduling). In addition, the strengths of the mining clusters approach over the
on-bench clustering approach often used in the literature were investigated and supported by studying
a hypothetical two-dimensional block model in this study.
Future research could focus on the capabilities of data mining algorithms to create mining clusters
in one step. In addition, investigating problems involving uncertainty could be another topic for future
research. Clustering ensemble methods could also be used as a new field owing to their better results
compared with the output results of the basic individual clustering methods in terms of consistency,
robustness and performance.

Disclosure statement
No potential conflict of interest was reported by the authors.

ORCID
Javad Gholamnejad http://orcid.org/0000-0003-3480-9644

References
Abbasi, S. O., S. Nejatian, H. Parvin, V. Rezaie, and K. Bagherifard. 2019. “Clustering Ensemble Selection Considering
Quality and Diversity.” Artificial Intelligence Review 52 (2): 1311–1340.
Alizadeh, H., B. Minaei-Bidgoli, and H. Parvin. 2014. “Cluster Ensemble Selection Based on a New Cluster Stability
Measure.” Intelligent Data Analysis 18 (3): 389–408.
Askari-Nasab, H., K. Awuah-Offei, and H. Eivazy. 2010. “Large-Scale Open Pit Production Scheduling Using Mixed
Integer Linear Programming.” International Journal of Mining and Mineral Engineering 2 (3): 185–214.
 ENGINEERING OPTIMIZATION 15

Askari-Nasab, H., M. Tabesh, and M. Badiozamani. 2010. “Creating Mining Cuts Using Hierarchical Clustering and
Tabu Search Algorithms.” In MININ 2010 Proceeding of the 4th International Conference on Mining Innovation, June,
Vol. 23, 534.
Ayad, H. G., and M. S. Kamel. 2008. “Cumulative Voting Consensus Method for Partitions with Variable Number of
Clusters.” IEEE Transactions on Pattern Analysis and Machine Intelligence 30 (1): 160–173.
Badiozamani, M. M., and H. Askari-Nasab. 2012. “Towards Integration of Oil Sands Mine Planning with Tailings and
Reclamation Plans.” In Tailings and Mine Waste Conference, Colorado State University, Keystone, Colorado, USA,
609–622.
Bagherinia, A., B. Minaei-Bidgoli, M. Hossinzadeh, and H. Parvin. 2019. “Elite Fuzzy Clustering Ensemble Based on
Clustering Diversity and Quality Measures.” Applied Intelligence 49 (5): 1724–1747.
Ben-Awuah, E., and H. Askari-Nasab. 2012. “A Mathematical Programming Model for Oil Sands Production Schedul-
ing and Waste Management.” Paper presented at the Tailings and Mine Waste Conference, Colorado State University,
Keystone, Colorado, USA, 593–608.
Ben-Awuah, E., H. Askari-Nasab, and K. Awuah-Offei. 2012. “Production Scheduling and Waste Disposal Planning for
Oil Sands Mining Using Goal Programming.” Journal of Environmental Informatics 20 (1): 20–33.
Busnach, E., A. Mehrez, and Z. Sinuany-Stern. 1985. “A Production Problem in Phosphate Mining.” Journal of the
Operational Research Society 36 (4): 285–288.
Chanda, E. K., and J. Ricciardone. 2002. “Long Term Production Scheduling Optimisation for a Surface Mining Oper-
ation: An Application of MinemaxTM Scheduling Software.” International Journal of Surface Mining, Reclamation
and Environment 16 (2): 144–158.
Davis, L. 1991. Handbook of Genetic Algorithms. New York: Van Nostrand Reinhold.
Fahad, A., N. Alshatri, Z. Tari, A. Alamri, I. Khalil, A. Y. Zomaya, S. Foufou, and A. Bouras. 2014. “A Survey of Clustering
Algorithms for Big Data: Taxonomy and Empirical Analysis.” IEEE Transactions on Emerging Topics in Computing
2 (3): 267–279.
Fred, A. L., and A. K. Jain. 2005. “Combining Multiple Clusterings Using Evidence Accumulation.” IEEE Transactions
on Pattern Analysis and Machine Intelligence 27 (6): 835–850.
Fred, A. L., and A. K. Jain. 2006. “Learning Pairwise Similarity for Data Clustering.” 18th International Conference on
Pattern Recognition (ICPR’06) 1: 925–928.
Fred, A., and A. Lourenço. 2008. “Cluster Ensemble Methods: From Single Clusterings to Combined Solutions.” In
Supervised and Unsupervised Ensemble Methods and Their Applications, edited by J. Kacprzyk, 3–30. Berlin: Springer.
Gershon, M. E. 1983. “Optimal Mine Production Scheduling: Evaluation of Large Scale Mathematical Programming
Approaches.” International Journal of Mining Engineering 1 (4): 315–329.
Gershon, M. E., and F. H. Murphy. 1989. “Optimizing Single Hole Mine Cuts by Dynamic Programming.” European
Journal of Operational Research 38 (1): 56–62.
Gonzalez, T. F. 1982. “On the Computational Complexity of Clustering and Related Problems.” In System Modeling and
Optimization, edited by R. F. Drenick and F. Kozin, 174–182. Berlin: Springer.
Goodfellow, R. C., and R. Dimitrakopoulos. 2016. “Global Optimization of Open Pit Mining Complexes with
Uncertainty.” Applied Soft Computing 40: 292–304.
Goodwin, G. C., M. M. Seron, R. H. Middleton, M. Zhang, B. F. Hennessy, P. M. Stone, and M. Menabde. 2006.
“Receding Horizon Control Applied to Optimal Mine Planning.” Automatica 42 (8): 1337–1342.
Hoerger, S., J. Bachmann, K. Criss, and E. Shortridge. 1999. “Long Term Mine and Process Scheduling at Newmont’s
Nevada Operations.” In Proceedings of the 28th Symposium on the Application of Computers and Operations Research
in the Mineral Industry, 739–748. Golden, CO: Colorado School of Mines.
Huang, D., C. D. Wang, and J. H. Lai. 2018. “Locally Weighted Ensemble Clustering.” IEEE Transactions on Cybernetics
48 (5): 1460–1473.
Jain, A. K. 2010. “Data Clustering: 50 Years Beyond K-Means.” Pattern Recognition Letters 31 (8): 651–666.
Jélvez, E., N. Morales, P. Nancel-Penard, J. Peypouquet, and P. Reyes. 2016. “Aggregation Heuristic for the Open-Pit
Block Scheduling Problem.” European Journal of Operational Research 249 (3): 1169–1177.
Kawahata, K., P. Schumacher, and K. Criss. 2016. “Large-Scale Mine Production Scheduling Optimisation with Mill
Blending Constraints at Newmont’s Twin Creeks Operation.” Mining Technology 125 (4): 249–253.
Kawahata, K., P. Schumacher, and R. Hufford. 2013. “Mine Production Scheduling Optimization at Newmont’s Twin
Creeks Mine.” Mining Engineering 65 (6): 49–54.
Klingman, D., and N. Phillips. 1988. “Integer Programming for Optimal Phosphate-Mining Strategies.” Journal of the
Operational Research Society 39 (9): 805–810.
Koushavand, B., H. Askari-Nasab, and C. V. Deutsch. 2014. “A Linear Programming Model for Long-Term Mine Plan-
ning in the Presence of Grade Uncertainty and a Stockpile.” International Journal of Mining Science and Technology
24 (4): 451–459.
Kuncheva, L. I., and S. T. Hadjitodorov. 2004. “Using Diversity in Cluster Ensembles.” In 2004 IEEE International
Conference on Systems, Man and Cybernetics (IEEE Cat. No. 04CH37583), Vol. 2, 1214–1219. IEEE, October.
Maulik, U., and S. Bandyopadhyay. 2000. “Genetic Algorithm-Based Clustering Technique.” Pattern Recognition 33 (9):
1455–1465.
16 R. LOTFIAN ET AL.

Minaei-Bidgoli, B., H. Parvin, H. Alinejad-Rokny, H. Alizadeh, and W. F. Punch. 2014. “Effects of Resampling Method
and Adaptation on Clustering Ensemble Efficacy.” Artificial Intelligence Review 41 (1): 27–48.
Mojarad, M., S. Nejatian, H. Parvin, and M. Mohammadpoor. 2019. “A Fuzzy Clustering Ensemble Based on Cluster
Clustering and Iterative Fusion of Base Clusters.” Applied Intelligence 49 (7): 2567–2581.
Mojarad, M., H. Parvin, S. Nejatian, and V. Rezaie. 2019. “Consensus Function Based on Clusters Clustering and Iter-
ative Fusion of Base Clusters.” International Journal of Uncertainty, Fuzziness and Knowledge-Based Systems 27 (1):
97–120.
Nazari, A., A. Dehghan, S. Nejatian, V. Rezaie, and H. Parvin. 2019. “A Comprehensive Study of Clustering Ensemble
Weighting Based on Cluster Quality and Diversity.” Pattern Analysis and Applications 22 (1): 133–145.
Osanloo, M., J. Gholamnejad, and B. Karimi. 2008. “Long-Term Open Pit Mine Production Planning: A Review of
Models and Algorithms.” International Journal of Mining, Reclamation and Environment 22 (1): 3–35.
Parvin, H., and B. Minaei-Bidgoli. 2013. “A Clustering Ensemble Framework Based on Elite Selection of Weighted
Clusters.” Advances in Data Analysis and Classification 7 (2): 181–208.
Parvin, H., and B. Minaei-Bidgoli. 2015. “A Clustering Ensemble Framework Based on Selection of Fuzzy Weighted
Clusters in a Locally Adaptive Clustering Algorithm.” Pattern Analysis and Applications 18 (1): 87–112.
Parvin, H., B. Minaei-Bidgoli, H. Alinejad-Rokny, and W. F. Punch. 2013. “Data Weighing Mechanisms for Clustering
Ensembles.” Computers & Electrical Engineering 39 (5): 1433–1450.
Ramazan, S. 2007. “The New Fundamental Tree Algorithm for Production Scheduling of Open Pit Mines.” European
Journal of Operational Research 177 (2): 1153–1166.
Ramazan, S., K. Dagdelen, and T. B. Johnson. 2005. “Fundamental Tree Algorithm in Optimising Production Schedul-
ing for Open Pit Mine Design.” Mining Technology 114 (1): 45–54.
Ramazan, S., and R. Dimitrakopoulos. 2004. “Recent Applications of Operations Research and Efficient MIP Formu-
lations in Open Pit Mining.” Society of Mining, Metallurgy, and Exploration, Inc. Transactions 316: 73–78.
Rashidi, F., S. Nejatian, H. Parvin, and V. Rezaie. 2019. “Diversity Based Cluster Weighting in Cluster Ensemble: An
Information Theory Approach.” Artificial Intelligence Review 52: 1341–1368.
Samanta, B., A. Bhattacherjee, and R. Ganguli. 2005. “A Genetic Algorithms Approach for Grade Control Planning in a
Bauxite Deposit.” In Proceedings of the 32nd International Symposium on Applications of Computers and Operations
Research in the Mineral Industry, 337–342.
Sinha, A., and P. K. Jana. 2016. “Clustering Algorithms for Big Data: A Survey.” In The Human Element of Big Data,
edited by G. S. Tomar, N. S. Chaudhari, R. S. Bhadoria, and G. C. Deka, 155–174. Portland, OR: Chapman and
Hall/CRC.
Smith, M. L. 1999. “Optimizing Inventory Stockpiles and Mine Production: An Application of Separable and Goal
Programming to Phosphate Mining Using AMPL/CPLEX.” CIM Bulletin 92 (1030): 61–64.
Smith, M. L., and S. J. Wicks. 2014. “Medium-Term Production Scheduling of the Lumwana Mining Complex.”
Interfaces 44 (2): 176–194.
Strehl, A., and J. Ghosh. 2002. “Cluster Ensembles—A Knowledge Reuse Framework for Combining Multiple Parti-
tions.” Journal of Machine Learning Research 3: 583–617.
Tabesh, M., and H. Askari-Nasab. 2011. “Two-Stage Clustering Algorithm for Block Aggregation in Open Pit Mines.”
Mining Technology 120 (3): 158–169.
Topchy, A., A. K. Jain, and W. Punch. 2003. “Combining Multiple Weak Clusterings.” In Third IEEE International
Conference on Data Mining, 331–338. IEEE, November.
Yousefnezhad, M., S. J. Huang, and D. Zhang. 2018. “WoCE: A Framework for Clustering Ensemble by Exploiting the
Wisdom of Crowds Theory.” IEEE Transactions on Cybernetics 48 (2): 486–499.