Long-term production scheduling of open pit mines using particle swarm and bat algorithms under grade uncertainty

important and integral part of the planning process for any open pit mine. It aims to define an extraction sequence of the mineralized material from the ground that produces the maximum possible discounted profit, i.e. NPV, while satisfying a set of physical and operational constraints. In the conventional approaches, the planning process usually starts with the construction of a geological block model that divides the orebody and the surrounding rock into three-dimensional arrays of regular, usually identical-sized, blocks. A set of attributes such as grade, specific gravity etc. is then assigned to each of these blocks, estimated using some form of spatial interpolation technique e.g. kriging, inverse distance method etc. and the drill-hole sample data. The blocks are then divided into two groups – waste and ore blocks. The blocks for which the prospective profit exceeds the processing cost are categorized as ore, while the rest are the waste blocks. An economic value is then assigned to each of these individual blocks by taking into account the group to which they belong, their respective estimated grade or metal content, the recovery, and the economic parameters such as metal price, mining cost, and processing costs. The ultimate pit limits (UPLs) are then determined using this economic block model, followed by the more complex production scheduling problem to define the most profitable annual extraction sequence of the blocks lying inside the ULP while satisfying different physical and operational constraints. A major drawback of this approach is its assumption that all the input parameters are known with certainty, while on the contrary a certain degree of uncertainty is almost always associated with these parameters, and ignoring this may result into unrealistic and erroneous scheduling decisions. The uncertainty in the input parameters may be caused by different economic and technical factors such as metal prices, currency exchange rates, mining and processing costs, and block grades. The graderelated uncertainty is considered to be the most important source of uncertainty for the scheduling process for open pit mines and is the main focus of study in this paper. The grade uncertainty stems from the fact that the grade values of the individual blocks are estimated using very limited drill-hole sample data, and usually a significant and variable level of uncertainty is associated with each of these estimated values. Geostatistical conditional simulation techniques provide a framework to generate multiple, equally probable realizations of the orebody’s Long-term production scheduling of open pit mines using particle swarm and bat algorithms under grade uncertainty

Long-term production scheduling is an important and integral part of the planning process for any open pit mine.It aims to define an extraction sequence of the mineralized material from the ground that produces the maximum possible discounted profit, i.e.NPV, while satisfying a set of physical and operational constraints.In the conventional approaches, the planning process usually starts with the construction of a geological block model that divides the orebody and the surrounding rock into three-dimensional arrays of regular, usually identical-sized, blocks.A set of attributes such as grade, specific gravity etc. is then assigned to each of these blocks, estimated using some form of spatial interpolation technique e.g.kriging, inverse distance method etc. and the drill-hole sample data.The blocks are then divided into two groups -waste and ore blocks.The blocks for which the prospective profit exceeds the processing cost are categorized as ore, while the rest are the waste blocks.An economic value is then assigned to each of these individual blocks by taking into account the group to which they belong, their respective estimated grade or metal content, the recovery, and the economic parameters such as metal price, mining cost, and processing costs.The ultimate pit limits (UPLs) are then determined using this economic block model, followed by the more complex production scheduling problem to define the most profitable annual extraction sequence of the blocks lying inside the ULP while satisfying different physical and operational constraints.
A major drawback of this approach is its assumption that all the input parameters are known with certainty, while on the contrary a certain degree of uncertainty is almost always associated with these parameters, and ignoring this may result into unrealistic and erroneous scheduling decisions.The uncertainty in the input parameters may be caused by different economic and technical factors such as metal prices, currency exchange rates, mining and processing costs, and block grades.The graderelated uncertainty is considered to be the most important source of uncertainty for the scheduling process for open pit mines and is the main focus of study in this paper.The grade uncertainty stems from the fact that the grade values of the individual blocks are estimated using very limited drill-hole sample data, and usually a significant and variable level of uncertainty is associated with each of these estimated values.Geostatistical conditional simulation techniques provide a framework to generate multiple, equally probable realizations of the orebody's Long-term production scheduling of open pit mines using particle swarm and bat algorithms under grade uncertainty Long-term production scheduling of open pit mines using particle swarm and bat algorithms attributes to quantify both the local variability and the potential uncertainty (Deutsch and Journel, 1997;Goovaerts, 1997).The availability of these techniques provides the opportunity to integrate the potential uncertainty into the scheduling process in order to generate a better production schedule in terms of achievable maximum NPV while minimizing the risk of deviations from production target.Different stochastic programming models for integrating grade uncertainty into the scheduling process have been proposed in the literature to produce better production schedules (Ramazan and Dimitrakopoulos, 2004;Golamnejad et al., 2006;Menabde et al., 2007).However, 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 and the potential uncertainty in the input data.This entails scheduling a large combinatorial optimization problem that may be extremely difficult and computationally expensive to solve using conventional techniques such as the branch and bound algorithm.In recent years a new class of computationally less expensive algorithms (metaheuristic techniques) for solving the mine design and production scheduling problem has attracted the attention of researchers, such as genetic algorithms (Denby andSchofield, 1994, 1996), simulated annealing (Kumral andDowd, 2002, 2005), and ant colony optimization (Sattarvand andNiemann-Delius, 2009, 2013).These techniques do not guarantee the optimality of the final solution that they produce, but they can produce better quality solutions with relatively less computational cost.
In this paper an extension of the framework previously presented by authors (Khan and Niemann-Delius, 2014) for long-term production scheduling for open pit mines is presented.The extension was made to integrate grade uncertainty in the optimization process.Two different metaheuristic techniques based on the particle swarm algorithm and the bat algorithm are proposed to solve a particular variant of the abovementioned production scheduling problem, i.e. two-stage stochastic programming model with recourse (Birge and Louveaux, 2011).
The paper is organized as follows.A two-stage stochastic linear programming model with recourse of the production scheduling problem of the open pit mines is presented, followed by a discussion of the standard PSO and bat algorithm.A heuristic procedure for producing the initial population of feasible random solutions is then presented, with the proposed procedure for applying the PSO and bat algorithm to the open pit production scheduling problem.A case study is also presented.
The two-stage stochastic (mixed integer) model with recourse formulation of the open pit mine scheduling problem, similar to the one proposed in Ramazan and Dimitrakopoulos (2004) can be represented as follows: [1] The objective function is composed of two parts.The first part is for the maximization of the expected NPV of the production schedule, and second part represents the minimization of the expected recourse costs whenever the stochastic constraints (Equations [6])-[9] are violated.

Subject to:
Reserve constraints: A block cannot be mined more than once [2] Slope constraints: Each block can only be mined if its predecessors are already mined in or before period t.
[3] i= 1,2 ….N ; t= 1,2 ….T; where j (set of predecessors blocks of block i).Mining capacity: The total material mined during each period should be within the predefined upper and lower limits.where BV is 0 is the undiscounted economic value of block i according to simulation s.

E{( NPV) i
t } is the expected discounted economic value of block i if mined in period t and is calculated as follows: [11] x it are the first-stage scenario-independent binary decision variables, which take the value unity if block i is mined in period t, and zero otherwise.
p t o -, p t o + are the discounted unit costs for shortage or surplus ore produced in period t respectively, and are calculated as: [12] p t m -, p t m + are the discounted unit costs for shortage or surplus metal produced in period t respectively and are calculated as: [13] y ts 0-, y ts 0+ are the second-stage scenario-dependent continuous variables representing the shortage or surplus amount of ore produced in period t if scenario s occurs, respectively.
y ts m-, y ts m+ are the second-stage scenario-dependent continuous variables representing the shortage or surplus amount of metal produced in period t if scenario s occurs, respectively.The second-stage (recourse) decision variables are dependent on the outcome of the orebody realizations and of the first-stage decision variables i.  Kennedy andEberhart in 1995 (Eberhart andKennedy, 1995;Kennedy and Eberhart, 1995).The PSO algorithm mimics the swarming behavior of individuals living together in groups, e.g.fish schools, bird flocks, animal herds etc.The 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.A random starting position and random velocity are assigned to each particle of the swarm.The velocity and position of each particle is then updated during the iterative process using Equations [14] and [15] respectively to find the optimum or near-optimum solution. [14] [15] where x i,t ,x i,t-1 represent the current and previous positions of particle i v i,t ,v i,t-1 represent the current and previous velocities of particle i p i,t is the personal best position experienced by particle i until iteration t g i,t is the global best position experienced by the population of particles until iteration t is the inertia weight used to control the influence/contribution of the particle's previous velocity to the current velocity, with its value fixed to 0.7298 r 1 , r 1 represent uniform random real numbers in the range of (0, 1) c 1 , c 2 are called the acceleration coefficients and are used to control the influence of cognitive and social terms on the particle's velocity, with the value commonly fixed to 1.49445 (Clerc and Kennedy, 2002).
® Initialize a population (swarm) of particles with random initial positions x i,0 and random velocities v i,0 in the search space.® Initialize each particle's personal best position p 1,0 to x i,0 its initial position ® Calculate the fitness value of each particle at its initial position x i,o and determine the initial global best position g o ® While (t < maximum number of generations) ® For all particles do ® Update the particle's velocity and position using Equations [14] and [15] respectively.® Calculate the fitness value of each particle at its current position x i,t ® If fitness (x i,t ) is better than the fitness (p i,t-1 ) )then where f i represents the current frequency of particle i, and f min , f max are the minimum and maximum allowable frequencies respectively.Initially a random frequency drawn from the uniform distribution in the interval [f min , f max ] is assigned to each bat.
[0,1] is a random number drawn from the uniform distribution in the mentioned interval, x * is the current global best position/solution, x i t+1 and x i t are the current and previous positions of particle i, v i t+1 and v i t are the current and previous velocities of particle i.For the local search part (steps 6 to 9), after selecting one solution among the current selected best solutions a new solution is generated using random walk. [19] where [-1,1] is a uniform random number and A t -is the average loudness of all the bats in the current time step.The pulse emission rate and the loudness of the bats are updated as the iteration proceed using the following equations: [20] [21] where A i t+1 and A i t are the current and previous loudness of bat i, and are constants, and r i t+1 and r i o are the current and initial pulse rates.The values of A i t and r i o can typically be [1,2] and [0,1] respectively.Yilmaz and Cengiz ( 2014) proposed a modified bat algorithm (MBA) to improve the poor exploration capabilities of the BA.They proposed to assign pulse emission rate r and loudness A to each dimension of the solution separately instead of assigning a single value to all the dimensions of a solution, and proposed the following procedure for updating the position along a certain dimension j of a solution: [22] where A j -t represents average loudness of dimension j of all solutions at time t, and u indicates a solution selected among the best solutions.The loudness and pulse rate are updated as: [23] [24] In the MBA a candidate solution is included in the population if r and < A A 'greedy' heuristic procedure has been used to generate the initial population of feasible solutions.The process starts by generating a list of the free blocks that are available to be mined in the current period as their predecessors are either already scheduled to be mined in the current or prior period or due to their position in 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 Long-term production scheduling of open pit mines using particle swarm and bat algorithms continues for the current period until the mining and processing capacities are satisfied 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, but more preference was given to satisfy the processing capacity constraint.While selecting an ore block more preference was given to the blocks with highest probability of being ore by considering the set of simulated values for that block and the predefined cut-off grade.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 is chosen at random from the range between upper and average mining capacities allowed in a period.
To apply the continuous variant of PSO and bat algorithms to the long-term production scheduling problem, the problem has been turned into optimum depth determination problem.After generating a population of initial feasible random solutions using the heuristic procedure discussed earlier, the mining depth along each column for each period is determined by using a so-called encoding scheme.These depths are then updated using the PSO or bat procedure during the iterative process.After each iteration, each solution is back-transformed to the discrete space, but due to the one-dimensional nature of the PSO and bat algorithms the resulting solutions may be infeasible in terms of the required slope angles.A normalization procedure is then used to turn the infeasible solution into a feasible one in terms of the required slope angles.A constant penalty method has been used to deal with the violation of the capacity constraints, where a constant penalty is added to the objective function for the per ton violation of the mining capacity constraints, to decrease the quality of the infeasible solutions by using Equation [25]. [25] where p t M -, p t M + are the discounted unit costs for shortage or surplus rock (ore plus waste) produced in period t respectively, and are calculated as: [26] R t Mand R t M + represent the shortage or excess amount of rock (ore plus waste) produced in period t.The value of the penalty is problem-dependent and is obtained for each problem by trial and error.A detailed description of the procedure can be found in Khan and Niemann-Delius (2014).
The pseudo codes of applying the PSO and bat algorithms to the open pit mine scheduling problem using this proposed procedure are given in algorithms 3 and 4.
® Input: block model, economic and technical parameters ® Initialize a population (swarm) of particles with random initial positions x i,0 and random velocities v i,0 in the search space (i = 1, 2 …….. N) ® Initialize each particle's personal best position p 1,0 to its initial position x i,0 ® Calculate the fitness value of each particle at its initial position x i,0 and determine the initial global best position g o ® While (t < maximum number of generations) ® Encode each particle's current position ( x i,t ), its personel best position (p i,t ), and the populations' global best position (g t ) ® for all particles do ® Update the particle's velocity and position using the equation of the PSO algorithm ® Back-transform ® Normalize the solution ® Calculate the fitness value of the particle at its current position x i,t ® If fitness (x i,t ) is better than the fitness (p i,t-1 ) then ® Input: block model, economic and technical parameters ® Initialize bat population with random initial positions x i,0 and random velocities v i,0 in the search space (i = 1, 2 ... N). ® Initialize each bat's frequency f i , pulse rate r i , and the loudness A i ® Calculate the fitness value of each particle at its initial position x i,0 and determine the initial global best position x * ® While (t < maximum number of generations) ® Encode each particle's current position (x Long-term production scheduling of open pit mines using particle swarm and bat algorithms Long-term production scheduling of open pit mines using particle swarm and bat algorithms To check the capabilities and the efficiency of the proposed procedure for generating a long-term production schedule under the condition of grade uncertainty, a set of 15 simulated realizations of a copper deposit was used.Examples of these realizations are shown in Figure 1.These realizations were generated using sequential Gaussian simulation technique (SGS) (Deutsch and Journel, 1997).
A predefined fixed cut-off grade strategy has been used to determine the destination of a particular block.While using this fixed cut-off grade strategy a block may be categorized as an ore block according to some simulations, and a waste block according to others (Figure 2).This fact has been taken into account while determining the expected economic value of the blocks.
The planning and scheduling process for open pit mines traditionally starts with the determination of the ultimate pit limits (UPLs), which define the extent to which it is economically feasible to mine.The UPLs were determined by solving the following mathematical formulation: [27] Subject to: [28] where V i represents the expected economic value of block i, N represents the total number of blocks in the block model, x i represents a binary variable corresponding to block i, and P i represents the predecessor group of block i.The required slope angles were assumed to be 45° in all directions.In the next step to define the optimum extraction sequence of the blocks lying inside the predetermined UPLs, the technical and economic parameters mentioned in Table I were used.The length of each scheduling period was assumed to be one year.The upper and lower limits for mining, processing, and metal capacity were set to be within ±20% and ±10% of the average available quantity of rock, ore, and metal available within the predefined UPL for each period of the scheduling horizon respectively.
All the numerical experiments were completed on an AMD Phenom II X4 945 (3.00 GHz and 4 GB RAM) running under Windows 7. The two-stage stochastic programming formulation of the open pit production scheduling problem (as mentioned in the section 'Problem formulation') with recourse was solved using the commercial solver CPLEX, considering the blocks within the predefined UPL and the parameters mentioned in Table I.This solution will be used as a benchmark to assess the performance of the PSO and bat algorithms in terms of computational time and solution quality.
After conducting a series of experiments, a population of 50 solutions was found out to be working well for the problem under consideration.A set 50 solutions was generated using the greedy heuristic described in the section 'Initial solutions'.These same solutions were used during all the subsequent experiments.Due to the stochastic nature of the metaheuristic algorithms the problem was solved 50 times using different variants of the PSO and bat algorithms.The maximum number of iterations was used as the termination criterion and was fixed to 2000 after observing the convergence behaviour of both the algorithms.The average relative % gap (which will be used a measure of quality of the solutions produced by the metaheuristic techniques under study) between best solutions generated by the different variants of the metaheuristic algorithms, i.e.Z approx and the optimal solution found by CPLEX, i.e.Z CPLEX are calculated using the following equation and are reported in Table II. [29] The values of the penalties used to report the violation of the capacity constraints, i.e. mining, processing, and metal capacity constraints in the objective function, were found to have a great effect on the performance of both the algorithms.When higher values were used, both the algorithms showed very few signs of improvement during the iterative process in comparison to the situations when smaller values were used.Therefore the following two sets of penalties have been used during this study.
® Algorithmic penalties: During the optimization process these penalties have been used for enabling the algorithms to efficiently explore the solution space and for generating the desired final solution with desired risk profile.These penalties are usually smaller in value than the corresponding actual penalties and are determined by trial and error.® Actual penalties: These are the actual penalties that have been used to calculate the objective value of the final solution generated by the PSO or bat algorithm at the end of the iterative process.
For the problem under consideration, both variants of the PSO algorithm performed better than the variants of the bat algorithm in terms of the achieved average % gap as mentioned in Table III.The multi-start PSO algorithm showed the best performance in comparison to all its populationbased competitors in the current study in terms of the achieved average % gap, showing the effect of diversification on the overall performance of the PSO algorithm.During the multi-start process the particle positions and velocities are reinitialized to their initial values after a fixed number of iteration as a diversification scheme.Both variants of the bat algorithm showed poor performance when the standard velocity equation, (Equation [17]) with maximum frequency value equal to unity was used to update the virtual bats' positions.
Therefore the following modified equation was used to update the velocity of a bat: [30] After conducting a series of experiments, the algorithm was found to be working well when the value of the maximum frequency was fixed to 0.6 and was chosen from the interval 0.4-0.6.The bat and MBA algorithms performed well when was fixed to 0.4, with a relatively small standard deviation showing the robustness of the technique in comparison to both variants of the PSO algorithm.In general both the metaheuristic techniques can produce better quality results with less computational cost in comparison to the exact optimization algorithm, thereby providing the opportunity of generating solutions with different risk profiles according to the needs of a certain open pit mining operation, which is quite difficult and expensive using an exact algorithm.
An extension of the framework previously formulated by the authors (Khan and Niemann-Delius, 2014)   performance in terms of average % standard deviation for the problem under study.The numerical experiments further revealed that the proposed procedure can generate sufficiently good solutions of the problem at hand in a shorter period of time.This property of the proposed framework provides the opportunity to handle grade-related uncertainties in more efficient way, which is otherwise quite difficult with currently available commercial software on the current hardware.
by A. Khan Long-term production scheduling for open pit mines is a large-scale, complex optimization problem involving large data-sets, multiple hard and soft constraints, and uncertainty in the input parameters.Uncertainty in the input parameters may be caused by different geological, economic, or technical factors.The uncertainty caused by geological factors, which is commonly termed geological or grade uncertainty, is considered to be most important source of uncertainty among these factors.It is caused by the fact that the block grade values are estimated using very sparse drill-hole sample data and the actual grade values will only be known once a block is drilled and blasted.Geostatistical conditional simulation techniques provide a framework to quantify this geological/grade uncertainty by generating multiple, equally probable simulated realizations of the orebody.Different stochastic programming models have been proposed in recent years to integrate this grade uncertainty into the optimization process, but solving these models for actual-sized open pit mines is usually extremely difficult and computationally expensive.In this paper, two different computationally efficient population-based metaheuristic techniques based on particle swarm optimization (PSO) and the bat algorithm are used to solve one particular stochastic variant of the open pit mine scheduling problem, i.e. the twostage stochastic programming model with recourse for determining the longterm production schedule of an open pit mine under the condition of grade uncertainty.open pit, production scheduling, grade uncertainty, stochastic programming, metaheuristic technique.
Stochastic constraints: The upper and lower processing capacity constraints (Equations [6] and [7] and the upper and lower metal production constraints Equations [8] and [9] represents the scenario-dependent stochastic constraints.These constraints are modelled as soft constraints where the violation is penalized in the objective function.Total number of periods t: Time period index, t = 1, 2 ….T N: Total number of blocks i: Block index, i=1,2 ..N S: Total number of simulations of the orebody s: Simulation index s = 1, 2 … S E{( NPV) i 0 } is the undiscounted expected economic value of block i and is calculated as follows: [10] the stopping criteria are met: end whileThe bat algorithm (BA) is a stochastic population-based optimization technique first proposed byYang (2010).It is based on the echolocation ability of bats.The micro-bats use their echolocation capabilities to find prey and avoid obstacles, even in complete darkness.The micro-bats emits short-duration (typically in the range of 8 to 10 ms) loudLong-term production scheduling of open pit mines using particle swarm and bat algorithms 363 VOLUME 118 L sound pulses with constant frequency in the region of 25 kHz to 150 kHz and listen for the echo that bounces back from the surrounding objects to find the food or avoid the obstacles.They usually emit 10 to 20 such sound pulses per second and can increase the pulse emission rate to about 200 pulses per second as they get close to their prey.To transform these unique properties of bats into an optimization algorithm, Yang idealized the following rules(Yang, 2010;Yilmaz, 2014) :® All bats use their echolocation capabilities to find out their distance from a certain object and can differentiate between food/prey and background barrier in some way.® Bats fly randomly with velocity v i at position x i with a frequency f min , and can vary the wavelength and loudness A 0 of their emitted sound pulses to find the food.Depending on their distance from the prey, they can adjust the rate and wavelength or frequency of their emitted pulse.® The loudness varies from a large positive value A 0 to a minimum constant value A min .Each bat's current position is a potential solution to the optimization problem.The bats use the following rules to update their positions:

®
Initialize bat population with random initial positions x i,0 and random velocities v i,0 in the search space ® Initialize each bat's frequency f i , pulse rate r i , and the loudness A i ® Calculate the fitness value of each particle at its initial position x i,0 and determine the initial global best position x * ® While (t < maximum number of generations) ® Generate a new solution by updating the frequency, velocity, and position using Equations [16], [17], and [18] respectively ® If (rand > r i ® Select a solution among the best solutions ® Generate a local solution around the selected best solution ® end if ® If (r and < A i and f(x i ) is better than f(x * ) ) ® Accept new solution ® Increase r i , reduce A i ® end if ® Rank the bats and find the current best x * ® If the stopping criteria are met: end while has been presented.This extended framework can account for grade uncertainty in the production scheduling process for open pit mines.To check the efficiency of the proposed extended framework for finding an optimum or near-optimum solution, two different population-based metaheuristic techniques, the PSO and bat algorithm, have been used.Through numerical experiments it was learnt that both variants of the PSO algorithm performed better than the variants of the bat algorithm in terms average % gap achieved, but showed poor Long-term production scheduling of open pit mines using particle swarm and bat algorithms VOLUME 118 367 L e. x it .