Abstract
Determining an optimum long term production schedule is an important part of the planning process of any open pit mine; however, the associated optimization problem is demanding and hard to deal with, as it involves large datasets and multiple hard and soft constraints which makes it a large combinatorial optimization problem. In this paper a procedure has been proposed to apply a relatively new and computationally less expensive metaheuristic technique known as particle swarm optimization (PSO) algorithm to this computationally challenging problem of the open pit mines. The performance of different variants of the PSO algorithm has been studied and the results are presented.
1. Introduction
Production scheduling of the open pit mines is a difficult and complex optimization problem. It aims to define the most profitable extraction sequence of the mineralized material from the ground that produces maximum possible discounted profit while satisfying a set of physical and operational constraints. Block model representation of the ore body commonly plays the role of the starting point for the planning and production scheduling of the open pit mines. The block model divides the mineralized body and the surrounded rock into a three-dimensional array of regular size blocks. A set of attributes, for example, grade, tonnage, density, and so forth, are then assigned to each one of these blocks estimated using some form spatial interpolation technique, for example, kriging, inverse distance weighting method, and so forth, and the exploratory drill hole sample data. The blocks are then divided into two categories, that is, ore and waste blocks. The blocks whose prospective profit exceeds their processing cost are categorized as an ore block to be sent for processing once mined while the rest are the waste blocks. An economic value is then assigned to each block by taking into account its group; that is, it is categorized as ore or waste, the commodity price, mining, processing, and marketing costs. The next step is to define the final contour of the pit, that is, the limit to which it is economically feasible to mine by solving the ultimate pit limit problem (UPL). The mathematical formulation of the UPL problem can be defined as follows: where represents the economic value of block , represents the total number of blocks in the block model, represents a binary variable corresponding to block which takes the value 1 if block is inside the ultimate pit limits and 0 otherwise, and represents the predecessor group of block . The solution to this problem identifies the collection of blocks that will result into maximum possible undiscounted profit. Graph based Lerchs and Grossmann algorithm [1] or Max-flow algorithms [2] are commonly used for defining the ultimate pit limits. Depending on the size of the individual blocks and of the ore deposit, the UPL may contain thousands to millions of blocks that may have to be scheduled over a time horizon typically ranging from 5 to 30 years while considering different physical and operational constraints which makes it a large combinatorial optimization problem.
This problem has been extensively discussed in the literature since 1960s and a wide variety of approaches have been proposed for solving it. Some of these techniques include parametric or ultimate pit limit (UPL) based approaches [3–5], integer programming [6, 7], dynamic programming [8, 9], and blocks aggregation/clustering [10]. The main limitations of these techniques are that either they have high computational cost when applied to real sized problem or they can only handle a simplified version of the original problem. In the recent years a new class of computationally less expensive metaheuristic techniques have attracted the attention of several researchers to solve this problem such as Genetic Algorithms [11–13], Simulated Annealing [14, 15], and Ant Colony Optimization [16, 17]. Though these techniques do not guarantee the optimality of the final solution that they produce, their lower computational cost makes them attractive alternative choice over the computationally expensive exact optimization algorithms. The objective of this paper is to propose a procedure to apply a similar metaheuristic technique known as particle swarm optimization (PSO) to the open pit mine production scheduling problem and to check and compare the efficiency of its different variants. In the past Ferland et al. [18] have proposed an approach to apply particle swarm optimization to the capacitated open pit mine production scheduling problem considering only one mining capacity constraint during each period of the mine life. They used a genotype representation of the solutions based on the priority value encoding and proposed a decoding scheme to generate feasible phenotype solutions from the genotype solutions. A GRASP procedure was suggested to produce the initial population of genotype solutions that evolve through the solution space using global version of particle swarm algorithm. Different properties of the procedure were studied by its application to a group of two-dimensional (2D) small problems. In this paper a totally different procedure has been proposed to apply the particle swarm optimization algorithm to the more general three-dimensional (3D) case of the open pit mine production scheduling problem with upper and lower bounds on the processing and mining capacity constraints in each period of the scheduling horizon.
The remainder of the paper is organized as follows. In Section 2 the general integer programming formulation of the production scheduling problem of the open pit mines is presented. In Section 3 the standard PSO algorithm is discussed. In Section 4 the proposed procedure for applying PSO algorithm to open pit production scheduling problem is presented in detail. The results of the numerical experiments are presented in Section 5 followed by the conclusions in Section 6.
2. Problem Formulation
The general integer programming formulation of the open pit production scheduling problem with the objective to maximize the net present value of the operation while satisfying a set of physical and operational constraints similar to the one defined by Gaupp in [19] can be described as follows.
Objective Function. Maximize the net present value of the extraction sequence: where is the discounted economic value of block if mined in period ; represents the annual discount rate; = {1, if block is mined in period ; 0, otherwise}; = total number of blocks; = total number of periods; = time period index, ; = block index, .This is subject to the following.
Reserve Constraints. A block cannot be mined more than once:
Slope Constraints. Each block can only be mined if its predecessors are already mined in or before period :
Mining Capacity Constraints. The total material mined during each period should be within the predefined upper and lower mining limits: , represent the lower and upper limits of the available mining capacity in period , respectively (constant), and represents the block tonnage.
Processing Capacity Constraints. The total ore processed during each period should be within the predefined upper and lower processing limits: , represent the lower and upper limits of the available processing capacity in period , respectively (constant). is a constant indicator; its value is equal to 1 if block is categorized as ore and 0 otherwise.
In a typical open pit mine the UPL may contain thousands to millions of blocks that have to be scheduled over a time horizon typically ranging from 5 to 30 periods or may be more; the resulting integer formulation may contain thousands to millions of integer variables and constraints, which may be extremely difficult and expensive to solve. For instance, defining production schedule for a small open pit mine containing 10000 blocks which needed to be scheduled over 10 periods may require 1,000,00 binary decision variables which is beyond the capabilities of currently available commercial software on the current hardware.
3. Particle Swarm Optimization
Particle swarm optimization (PSO) algorithm is a stochastic population based optimization technique first proposed by Kennedy and Eberhart in 1995 [20, 21]. PSO is a nature inspired algorithm. It is based on the social interaction of individuals living together in groups, for example, bird flock, fish schools, animal herds, and so forth. PSO algorithm performs the search process by using a population (swarm) of individuals (particles). Each individual (particle) is a potential solution to the optimization problem. At the start a random starting position and random velocity are assigned to each particle of the swarm. The velocity and position of the individual particles are then iteratively adjusted using the following equations for fining better positions in the search space: where , = position of the particle during current and previous iterations, respectively, , = velocity of the particle during current and previous iterations, respectively, = personal best position experienced by particle until iteration , = global best position discovered by the population until iteration , = inertia weight used to control the contribution of the particle’s previous velocity, , = acceleration coefficients used to control the influence of cognitive and social terms on the particle’s current velocity, , = vectors of uniform random numbers between 0 and 1.The flow chart of the standard PSO algorithm is shown in Figure 1.
3.1. Guaranteed Convergence PSO (GCPSO)
This modified variant of the PSO algorithm was introduced by Van den Bergh to counter the problem of premature convergence to solutions that are not guaranteed to be the local optimum. The following modified equation was proposed to update the velocity of the best particle of the population or a certain neighborhood: where represent the index of the best particle and is a vector of uniform random numbers. is a scaling factor used to perform a random search around the global best position; its value is modified during the iterative process using an adaptive strategy [22].
4. Particle Swarm Optimization Algorithm for Open Pit Mine Production Scheduling Problem
In this paper the continuous version of particle swarm algorithm and guaranteed convergence PSO algorithm (GCPSO) [23] (with a few modifications) has been applied to open pit mine production scheduling problem. There are two main reasons for using the continuous version of the said algorithms: firstly the open pit mine may contain thousands to millions of blocks inside the ultimate pit limit; making the scheduling decisions on the blocks level may be computationally expensive. Secondly the blocks in the same column of a block model are stacked on top of each other that defines a precedence relationship among these blocks which dictates that the block scheduling problem can be turned into optimum depth determination problem which makes the PSO calculations much faster and the implementation simple [24]. The flow chart of the proposed procedure for applying PSO algorithm to the open pit mine scheduling problem is shown in Figure 4 and will be discussed in more detail in the following sections.
4.1. Initial Solutions
To generate an initial population of random feasible solutions a sequential heuristic procedure has been used. The process starts by generating a list of the blocks that are available to be mined in the current period as their predecessor blocks are either already mined or due to their position in the block model they do not have one. A block is selected at random from the list and assigned to the current period and the list is updated again. The process continues for the current period until the mining and processing capacities are satisfied for it in the average sense before moving on to the next period. Mining and processing capacity constraints are handled as hard and soft constraints, respectively, during this process. The process stops when either the free blocks list is empty or the number of periods is finished for the current solution. To further diversify the generated solutions the annual mining capacity per period for each solution is chosen at random from the interval between average and maximum allowed mining capacity.
4.2. Solution Encoding
A solution encoding scheme is being used in an effort to apply the continuous version of particle swarm optimization algorithm to open pit production scheduling problem. The proposed encoding scheme determines the deepest block along a certain column to be mined in a certain period. This value (depth) is assigned to a variable corresponding to that column and period. The values of these depths variables are then updated using the standard PSO or GCPSO algorithm in each iteration in an effort to find better solutions.
4.3. Back Transform
After each iteration a back transform scheme is used to determine the period to which a block is being assigned. The procedure takes the depth variable for a certain column and period and determines all the blocks that are lying above that depth and below the depth defined by the PSO algorithm for the same column for the prior period. Obviously there is no upper limit for the first period. This process progresses from the first to the last period. During the optimization process it is ensured that depth for the successive periods along a column does not cross each other.
4.4. Constraint Handling
Particle swarm algorithm like other evolutionary algorithms does not have an explicit mechanism for constraint handling. Considering the special structure of the constraints the following two techniques have been used for constraint handling in this study.
4.4.1. Solution Normalization
Due to the one-dimensional nature of the PSO algorithm, the mining depths along each column for each period are determined independently from each other which may result in an infeasible solution in terms of the required slope angles as shown in Figure 2. To turn this back transformed infeasible solution into a feasible solution, a so-called normalization technique is being used.
4.4.2. Penalty Method for Capacity Constraint Violation
A constant penalty method [25] is used to deal with the mining and processing capacity constraints, where a constant penalty is added to the objective function for per ton violation of the capacity constraints, to decrease the quality of the infeasible solutions. The value of the penalty is problem dependent and is obtained for each problem using trial and error procedure.
5. Numerical Results
The proposed procedure for applying PSO algorithm to production scheduling problem of open pit mines has been implemented using Microsoft Visual studio 2010 (C++) programming environment. The capabilities and the efficiency of different variants of the PSO will be checked using this program. All the numerical experiments have been completed on AMD Phenom IIX 4 945 (3 GHz) and 4 GB Ram running under windows 7.
Two different data sets describing two different hypothetical copper deposits have been used to conduct the numerical experiments. The process started by first dividing the blocks into two categories, that is, ore and waste blocks, using a predefined fixed cutoff grade strategy followed by the definition of the economic block model. The ultimate pit limits for both the deposits were determined by solving the ultimate pit limit problems (as described in Section 1). The required slope angles were assumed to be 45° in all directions (Figure 3). The length of each scheduling period was assumed to be 1 year. The upper and lower limits for processing and mining capacity were set to be within ±15% of the average available quantity of ore and rock within the UPL for each period of the scheduling horizon, respectively. The discount rate was assumed to be 8% per year. The integer programming formulation of the production scheduling problem (as described in Section 2) using above-mentioned data sets and parameters was solved using commercial solver CPLEX; this solution was used as a benchmark to assess the performance of the different variants of the PSO algorithm in terms of computational time and solution quality. General information of the blocks lying inside the ultimate pit limits and about the solutions found by CPLEX is given in Table 1.
The efficiency of the following three different variants of the PSO and GCPSO algorithm has been checked during the numerical experiments.
Global (Gbest). It is the one where each particle is influenced by the best particle of the entire population.
Local (Lbest). Each particle is only influenced by the best particle of its immediate neighbours, where represents here the total number of indexed based neighbours of each particle of the population.
Multistart. Considering the fast convergence behaviour of the global variant of the PSO/GCPSO algorithm a Multistart start strategy has been applied as a diversification strategy, where the particles are reinitialized to their initial position whenever the algorithm does not show any improvement in Gbest value for a certain number of iterations (Figure 5).
(a)
(b)
(c)
During all the experiments the values of the inertia weight were set to 0.7298 and the values of the acceleration coefficients ( and ) were fixed to 1.49445 as proposed by Clerc and Kennedy in [26]. After conducting a series of experiments a population of 50 particles generated using the sequential heuristic procedure (described in Section 4.1) were found to be working well for both the problems. The performance of all the variants of the PSO algorithm was found out to be greatly affected by the numerical values of the penalties used during the iterative process. Using relatively higher values was causing the algorithm to get trapped in a local optimum without exploring better solution regions while on the other hand using very low penalties was causing excessive violation of the production capacity constraints. To deal with this problem the following two different sets of penalties have been used during the optimization process.
Algorithmic Penalties. These penalties were used during the iterative process for getting acceptable results in terms of the quality of the final solution produced. Their values are data dependent and needed to be determined using trial and error procedure. The appropriate values were found to be 7 and 9 $/ton for problems 1 and 2, respectively.
Actual Penalties. These are data independent penalties and were used to calculate fitness value of the final solution produced at the end of the iterative process.
During all the experiments, the maximum number of iterations was used as the termination criteria. The appropriate values were determined for each problem independently from each other after observing the convergence behavior of different variants of the PSO algorithm when applied to them (as shown in Figure 6) and were fixed to 30000 and 9000 for problems 1 and 2, respectively. These same values will be used in the case of GCPSO as well. The chances of getting better positions in the consecutive iteration were found out to be quite low so instead of using the adaptive strategy for updating the ρ value during the iterative process a constant value has been used. A value equal to 5 was found out to be working well for the problems under consideration. While applying the local variant of GCPSO algorithm the positions of the local and global best particle are updated using the GCPSO update equation once separated from the rest of the particles to avoid updating their position multiple times in the same iteration. Due to the stochastic nature of the PSO algorithm, each problem was solved 10 times using each variant of the PSO algorithm. The average relative % Gap between best solutions generated by the different variants of the PSO/GCPSO algorithm, that is, , and the optimal solution found by CPLEX, that is, , is calculated using the following equation and is reported in Tables 2, 3, 4, and 5: The global variant of the PSO algorithm is the one where each particle is influenced by the best particle of the entire population and is considered to be more susceptible to get trapped in a local optimum position and has shown similar kind of convergence behaviour which can be observed in Figure 6. On the other hand the local variants show relatively slower rate of convergence with different neighbourhood sizes. The Multistart strategy proved to be a better option to overcome the premature convergence problem of the global variants of both the PSO and GCPSO algorithms and on average has produced better solution with smaller relative gap as can be seen in Tables 2, 3, 4, and 5. The standard deviation of % Gap in almost all the cases has very small values showing the robustness of the procedure. Generally the average relative gap tends to increase with problem size. The computation time required by the algorithm depends on the factors defining and increasing the size of the required computations to be performed before the termination criteria are met such as population’s size, number of blocks, and number of periods. The computation time increases linearly with the number of iterations for a specific problem. The proposed procedure in all the cases can produce relatively better quality results in relatively shorter period of time in comparison to the exact optimization algorithms implied by CPLEX.
(a)
(b)
6. Conclusion
This paper presents a framework for applying a computationally efficient population based metaheuristic technique known as particle swarm optimization algorithm to the long term production scheduling problem of the open pit mines. Instead of making the scheduling decision on the block level the proposed approach turns the problem into depth finding problem and then applies the continuous version of the PSO algorithm and its different variants to find the optimum or near to optimum solution. A greedy and fast heuristic technique was used to generate the initial population of random feasible solutions which were transformed over several generations using particle swam procedure to get optimum or near to optimum solution of the production scheduling problem. By making comparison with the results obtained using CPLEX in terms of computational time and solution quality it was learned that the proposed procedure can produce better quality solutions in relatively shorter period of time with smaller % Gap and standard deviation showing the robustness of the procedure. The proposed procedure is more structured and quite flexible and can easily accommodate additional constraints, variable slopes angles, and grade and market uncertainties.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.