Phase Planning for Open Pit Coal Mines through Nested Pit Generation and Dynamic Programming

This paper presents a phase planning method specially designed for coal deposits with nearly horizontal, bedded coal seams. The geology of this type of deposit is modeled into a column model, instead of a block model, to avoid coal-rock mixing in blocks. A nested pit generation algorithm is developed for producing a series of nested, least-strip ratio pits with a column model as its input. The algorithm completely overcomes the troublesome gap problem. Taking the least-strip ratio pits as possible phase states, a dynamic programming formulation is proposed to simultaneously optimize the number of phases, the phase-pits, and the ultimate pit, with an objective of maximizing the net present value. The merits and capability of the proposed method are demonstrated through a case study on a large coal deposit.


Introduction
Large open pit mines are often mined in a number of phases, with intermediate pits for the phases referred to as phase-pits or pushbacks. e phases must be carefully planned since they provide a long-term strategic guide for the sequential development of a mine and for detailed production scheduling. In designing the phases, three elements must be determined: the number of phases, the phase-pits, and the ultimate pit. ese elements are interrelated and should be simultaneously optimized to maximize the net present value (NPV) of an open pit project.
Phase optimization, unlike production scheduling, has not been an extensively studied topic in open pit planning in recent years, and most of the research work had been done before the turn of the century. Lerchs and Grossmann first introduced the parametric analysis approach, where the block values of the block model are systematically changed and a pit optimization algorithm is executed repeatedly to obtain an optimal pit each time after the block values are changed, producing a series of nested pits [1]. ese pits can then be evaluated to choose the appropriate phase-pits. Technical parameterization of reserves is another method for nested pit generation, with an objective of finding the family of " technically optimal pits. " A technically optimal pit for a total volume, V, and ore quantity, Q, is the pit that has the highest quantity of mineral of interest among all pits of the same V and Q. A number of authors have elaborated on the mathematical formulations and solution algorithms for technical parameterization of reserves [2][3][4][5][6][7]. Dagdelen and Johnson formulated the pushback optimization problem into an integer programming (IP) model and used Lagrangian relaxation for solution [8]. All these methods have an inherent gap problem; that is, the size increment between two consecutive pits can be very large. Big gaps can inflict serious difficulties when the pits are used for phase planning or production scheduling, with the resulting solution far from the optimal or with no feasible solution at all. To overcome the gap problem, Wang and Sevim proposed a heuristic algorithm using a cone eliminating process [9]. e algorithm is capable of finding a series of nested, maximummetal pits with a pit increment almost the same as the user specified value and, thus, completely eliminates the gap problem. Meagher et al. formulated the pushback optimization problem as an IP model [10]. To facilitate the solution process and to overcome the gap problem, they solved the linear programming relaxation version of the IP model first and, then, applied a method known as " pipage rounding " to convert a fractional solution into an integral solution. e authors claimed that their approach completely overcomes the gap problem.
One of the most significant developments in the related and wider area of open pit planning, over the last two decades, is probably the incorporation of geological uncertainty in planning schemes. Conditional geostatistical simulation techniques have been developed that are capable of generating multiple, equally probable, realizations of an orebody [11][12][13][14]. ese realizations show the possible variations of the mineral content and corresponding tonnages of a deposit and provide the basis for quantifying geological uncertainty. Based on simulated orebody realizations, stochastic open pit planning techniques can be used to integrate geological uncertainty into planning processes, so as to allow some sort of geological risk management while maximizing the expected NPV. Several authors have incorporated geological uncertainty in their pushback/phase-pit optimizing approaches. Goodfellow and Dimitrakopoulos applied the simulated annealing algorithm to pushback design based on simulated orebody realizations [15]. eir goal was to modify an existing pushback design to better account for the joint local uncertainty in metal grades and material types, while remaining similar to the original design in terms of pushback tonnages and the tonnages sent to various destinations. eir case study on a copper mine showed that the approach achieved a 35-61% reduction in variability in terms of material quantities sent to the processes, leading to a reduced level of risk in the economic value of the design. Asad and Dimitrakopoulos presented a stochastic parametric maximum flow algorithm for pushback design under uncertainties in both metal content and commodity price [16]. ey used Lagrangian relaxation together with subgradient method to accommodate knapsack constraints for ore quantities in the pushbacks and addressed the gap problem by introducing a modification in the subgradient method to minimize the size difference between consecutive pushbacks. e authors applied the algorithm to a gold mine and compared the outcome with that from the conventional (deterministic) nested pit approach [17]. e case study demonstrated that the stochastic approach gave 30% more discounted cash flow, a 21% larger ultimate pit, and about 7% more metal than the conventional approach. Asad et al. later expanded the formulation to include multiple ore processing streams [18].
As mentioned above, the number of phases, the phasepits, and the ultimate pit should be optimized simultaneously. However, few approaches are capable of providing simultaneous solutions to the three elements, especially, for coal mines. Asad and Dimitrakopoulos' algorithm solves the phase-pits and ultimate pit simultaneously, but the gap problem is not completely overcome and the algorithm is for metal mines [16]. Gu et al. proposed a method capable of providing simultaneous solutions to the three elements [19]. It consists of a heuristic algorithm for generating a series of geologically optimal pits and a dynamic programming (DP) formulation for sequencing the pits, but the method is also for metal mines.
Based on the basic framework from a previous research for metal mines [19], we have developed a phase planning method particularly for open pit coal mines with nearly horizontal and bedded coal seams, with the objective of simultaneously optimizing the number of phases, the phasepits, and the ultimate pit.

Coal Deposit Modeling: the Column Model
Almost all open pit optimization formulations and solution algorithms take 3D block models of deposits as their geological inputs. However, most coal deposits consist of nearly horizontal (inclination angle smaller than 15°or so), bedded coal seams. As depicted in Figure 1, blocks at the coal-rock interfaces contain a mixture of coal and rock (e.g., the blocks outlined by the dotted lines), and such blocks constitute a significant portion of all the blocks containing coal. is will cause large errors in coal quantity calculation and subsequent economic evaluation, since whole blocks are classified as coal or waste. e problem is more pronounced in cases of multiple coal seams and/or the coal seams are thin.
In this study, we model this type of coal deposit into a " column model, " where the whole deposit is divided into vertical columns instead of blocks, with all columns having a square horizontal cross section of the same size. Figure 1 illustrates a vertical cross section of a column model with the columns numbered from 1 to 23. Each column in such a model is assigned a group of attributes defining the physical and chemical properties of all the coal seams along the column. e attributes generally include floor elevation and thickness of each coal seam and the unconsolidated layer at the central line of the column, and heat value, ash content, and sulfur content of each coal seam along the column. ese attributes may be estimated based on drill hole data using a method such as Kriging or inverse distance interpolation. With a column model, the coal quantity in any given volume (e.g., a pit or a cone) is calculated based on the thickness of each seam falling inside the volume on each column. e resulting coal quantity is much more accurate than that based on whole blocks with a block model. e reduction in coal-rock mixing with a column model, as compared with a block model, depends on the geometries of the coal seams (especially, their thicknesses) and the block and column sizes. For the coal deposit model used in the case study, we estimated that the reduction is at least 50% in terms of the total amount of rock mixed in coal and coal classified as waste.

Generation of Least-Strip Ratio Pits
e basic idea of phase optimization is first generating a candidate series of nested pits with a specified size increment and, then, selecting the best phase-pits (including the ultimate pit) from the series that maximize the NPV. e best candidate pits are the least-strip ratio pits, referred to as " least-SR pits " hereafter. A least-SR pit for a given coal quantity, Q, is defined as the pit that has the lowest strip ratio of all pits having the same Q. To overcome the gap problem, a cone eliminating algorithm is applied on a column model 2 Mathematical Problems in Engineering to generate a series of nested, least-SR pits. e basic logic of the algorithm is briefly described as follows.
Suppose that a series of n nested, least-SR pits, denoted as P { } n P 1 , P 2 , . . . , P n , is to be generated with a specified coal quantity increment of ΔQ, where P 1 is the smallest pit and P n the largest. Let Q i denote the coal quantity in P i . e algorithm starts with the largest pit, P n , which is created by projecting the pit walls from a closed surface boundary outline down to the lowest floor elevation of the lowest coal seam, according to the pre-determined slope angles in different directions or in different zones. e pit slope angles are determined beforehand through rock stability studies and are inputs to the algorithm. e surface boundary could be the boundary that encloses all exploration drill holes, or the property boundary for which the mining company has acquired the right of mining, or any boundary that is large enough to enclose the optimal ultimate pit within the acquired property. Within P n , the algorithm searches for a portion that contains a coal quantity of ΔQ and has the highest strip ratio, and then eliminates this portion from P n . e remaining part of P n after eliminating such a portion constitutes the next smaller pit, P n−1 , in the series, which has a coal quantity of Q n − ΔQ. To make P n−1 feasible with respect to pit slope constraints, the eliminated portion is constructed by combining upward cones (whose apexes point upwards) with shell inclination angles equal to the pit slope angles. From P n−1 , the same cone eliminating process is repeated to generate the nest smaller pit, P n−2 . e process continues until the coal quantity in the remaining part is equal to or below the coal quantity, Q 1 , specified for the smallest pit in the series, thus, resulting all the pits in the series.
With a column model, a pit is outlined by the bottoms and the tops of all the columns inside the pit, as depicted in Figure 2 for pit P i , where the bottom elevation of each column is equal to the elevation of the pit wall or pit bottom at the column center, and the top elevation of each column is equal to the elevation of ground surface at the column center. Without losing generality, suppose we have come to the point of generating P i−1 from pit P i . e cone eliminating process is outlined as follows.
Place the cone apex on the central line of a column at an elevation that is Δz higher than the bottom elevation of the column, as shown by Cone 1 in Figure 2. Calculate the quantities of coal, rock, and unconsolidated material in the cone inside pit P i . e cone's coal quantity is denoted as q c . If q c is not greater than the coal quantity increment, ΔQ, specified for the pit series, compute the cone's strip ratio and put the cone in an array. en, move the cone apex upwards along the central line of the same column by a distance of Δz, and do the same as for the previous cone. Continue the process of moving the cone upwards along the same column, until the cone's coal quantity, q c , is greater than ΔQ, or the cone's apex is above the ground surface, as depicted by the upward arrow and the dot lined cones in Figure 2. en, the cone is moved horizontally to another column inside pit P i , as shown by the horizontal arrow in Figure 2, and the entire process for the previous column is repeated. e above process continues until all the columns inside pit P i are traversed.
At the end of the above cone moving process, an array of J cones, each having a coal quantity q c ≤ΔQ, is obtained. Sort the cone array in order of descending strip ratio. A union of the first K cones in the sorted cone array is sought, such that the total coal quantity of the union is closest to ΔQ and not greater than 1.1ΔQ. In a cone union, the overlapping part between two or more cones is accounted only once. Such a union is obtained by sequentially combining the cones in the sorted array one at a time, starting with the first cone. e cones in the union are eliminated from pit P i and the remaining part is pit P i−1 , whose coal quantity is about ΔQ smaller than P i . Eliminating a cone from a pit is simply done by raising the bottom elevation of each column traversed by the cone up to the cone shell elevation at the column center. e algorithm outlined above is a heuristic one and does not guarantee that the resulting pit is the true least-SR pit. However, since the eliminated cones are the ones having the highest strip ratios, their union constitutes a volume whose strip ratio should be very close to the highest strip ratio of all volumes with the same coal quantity. us the remaining part should be very close to the true least-SR pit for its coal quantity.
In the above algorithm, each and every cone kept in the array has a coal quantity not greater than the specified increment, ΔQ, and the coal quantity of the eliminated cone union is controlled by an upper limit of 1.1ΔQ. us, the coal quantity increment between any consecutive pits in the generated pit series may be smaller than ΔQ, but cannot exceed 1.1ΔQ. e gap problem is, therefore, completely overcome. e up-moving step size, Δz, can affect the quality of the resulting pits: a smaller Δz generally gives better result (i.e., the resulting pits are closer to the true least-SR pits), but consumes more time and memory. Δz is an input parameter to the developed software and different values can be tried if necessary. From our experiments on different coal deposits using different Δz values, bench height, h (usually 10 m−15 m), is a good choice for Δz, and smaller values (e.g., h/2, h/4) make insignificant improvement on the resulting pits, but substantially increase the computing time. e column size has a similar effect. Blocks containing a mixture of coal and rock if a block model were to be used

Dynamic Programming Formulation for Phase Planning
Once a series of nested, least-SR pits, {P} n , is generated, the pits can then be sequenced using a DP model to optimize the phase plan. Figure 3 is a schematic illustration of the DP model, and for clarity of illustration, {P} n is assumed to contain 6 pits (in real-life instances, the number would be much larger). e horizontal axis represents the stage variable t with each stage corresponding to a phase, and the maximum stage number is the number of pits in {P} n . e vertical axis represents the state variable P with each state corresponding to a pit in {P} n , depicted by a circle in Figure 3. e states (pits) of a given stage (phase) are the possible phase-pits for that phase. An arrow represents a state transition from a pit of a phase to a pit of the next phase. Since any phase-pit of phase t is the result of expansion (through mining) of a smaller phase-pit of the preceding phase, t − 1, a state transition can only go upwards from a pit of a stage to one of the larger pits of the succeeding stage. at is why the starting (lowest) state of stage t corresponds to pit P t in {P} n (t � 1, 2,. . .,n), and the lower-right half of the diagram is void.
A path starting with the origin and ending at any pit in Figure 3 is a possible phase plan scenario. For example, path 0⟶P 2 ⟶P 4 ⟶P 6 , as shown by the thick arrows, represents such a phase plan: the number of phases is 3 (since the path ends at phase 3); the phase-pits for phases 1, 2, and 3 are pits P 2 , P 4 , and P 6 , respectively; and the ultimate pit is P 6 . e path with the highest NPV is the optimal phase plan, which can be found by economically evaluating all the paths. e following is a DP formulation for finding the best path.
In general, suppose that pit P i of stage t is being evaluated. P i can be transited from those smaller-than-P i pits of the preceding stage, t − 1. When pit P i of stage t is transited from pit P j of stage t − 1 (t −1 ≤ j ≤ i− 1), the quantities of coal, rock, and unconsolidated material mined in phase t, denoted as Q t,i (t − 1, j), W t,i (t − 1, j), and U t,i (t − 1, j), respectively, are calculated as where Q i , W i , and U i are the quantities of coal, rock, and unconsolidated material that can be mined from pit P i , respectively. ey are quantities after coal recovery and waste mixing incurred in mining operations are taken into account.
Suppose that the mining company has a coal processing plant and the salable product is clean coal. Such a transition brings a total revenue of V t,i (t − 1, j) and cost of C t,i (t − 1, j) for phase t.  Figure 2: Illustration of the cone eliminating process for generating nested, least-SR pits.  where r p is the coal recovery rate of the processing plant; p is the coal price; and c m , c p , c w , and c u are the unit costs of coal mining, coal processing, rock mining, and unconsolidated material stripping, respectively. Let y t,i (t − 1, j) denote the time (in years) required to make such a transition, and assume that the coal mining, waste removing, and coal processing capacities match one another. en, where A is the annual coal mining capacity. y t,i (t − 1, j) may not be an integer number of years, and let L t,i (t − 1, j) be the integer part of it. e average annual revenue and cost for each of the L t,i (t − 1, j) years, denoted by v t,i (t − 1, j) and c t, e revenue and cost for the remaining decimal part, denoted by a t, i (t − 1, j) and b t,i (t − 1, j), respectively, are Following the transition, the cumulative time to arrive at pit P i of stage t after finishing mining phase t, denoted by Y t, where Y t−1,j is the cumulative time to arrive at pit P j of the preceding stage, t − 1, following the best path. Y t−1,j has been calculated in evaluating the states of the preceding stage, t − 1. erefore, when pit P i of stage t is transited from pit P j of stage t − 1, the cumulative NPV realized at pit P i of stage t, after t phases of production, is given by where NPV t−1,j is the cumulative NPV at pit P j of stage t − 1, following the best path, which has been calculated in evaluating the states of the preceding stage, t − 1; ρ p and ρ c are the escalation rates of coal price and production cost, respectively; and d is the discount rate. As stated before, pit P i of stage t may be transited from all the smaller-than-P i pits of the preceding stage, t − 1. Obviously, when pit P i of stage t is transited from a different pit of stage t − 1, the quantities mined and processed in phase t will be different, and the revenue, cost, and time length will also be different. Consequently, different transitions (decisions in DP) give different cumulative NPVs at pit P i of stage t. e transition with the highest cumulative NPV is the best transition (optimal decision in DP) and, thus, the recursive objective function is When the pits of stage 1 are evaluated, all the pits can only be transited from the initial state at t � 0 (the origin in Fire 3). Initial conditions at the initial state are Starting from the first stage, the pits are evaluated forwards stage by stage, until all the pits of all stages are evaluated. e best transitions and the associated cumulative NPVs are obtained for all the pits of all stages. en, find the pit that has the highest cumulative NPV of all pits of all stages. is pit is the best ultimate pit, and the stage at which it is found indicates the best number of phases. en, starting from the best ultimate pit and tracing the best transitions backwards to the first stage, the optimum path (optimal policy in DP) is found, and the pits along this path indicate the best phase-pits of the corresponding phases.
us, the number of phases, the phase-pits, and the ultimate pit are simultaneously optimized.
is is a forward and open-ended DP formulation.
Mathematical Problems in Engineering 5

Case Study
A software package has been developed based on the above model and algorithm and was used in a case study on a large coal deposit in northern China. e topography of the coal field is nearly flat with a maximum relief of less than 20 m. e surface boundary of the planning area is about 7800 m long and 4500 m wide. From drill hole information, 8 coal seams were identified, with thicknesses varying from around 1 m to around 40 m and densities between 1.28 and 1.31 t/ m 3 . e deposit was modeled into a column model having around 39000 columns, each having a horizontal cross section of 30m × 30 m. e attributes of each column, mainly the thickness and elevation of each coal seam along the column, were estimated based on drill hole data using an inverse distance interpolation method particularly designed for this study. e coal reserve within the planning boundary was estimated to be some 900 Mt. For a deposit of this scale, the annual coal mining rate was assumed to be 20 Mt of run-of-mine coal.
e coal is to be sold without processing. e time span of a single phase is generally more than 5 years to avoid complications associated with frequent transitions between phases. erefore, in generating the least-SR pits, the coal quantity of the smallest pit, Q 1 , was set to 100 Mt (5-year production), and the coal quantity increment, ΔQ, to 20 Mt. A maximum pit slope of 25°was used. With these parameter values and the column model as inputs, 40 nested pits were generated by applying the cone eliminating algorithm described above. e coal quantity increment between any two consecutive pits in the generated pit series has a very small variation (20.00 Mt to 20.17 Mt), indicating that the algorithm has produced an evenly spaced series of nested pits with increments virtually equal to the specified value.
is is a verification of the algorithm's capability of completely overcoming the gap problem. e phase plan was optimized with the generated pit series and the parameter values in Table 1 as inputs to the DP model. e best (highest-NPV) phase plan consists of 4 phases and Table 2 summarizes the major quantities to be mined in the phases. With an annual coal production of 20 Mt, each of the first three phases has a time span of 10 years and the fourth phase 11 years, giving a total mine life of 41 years. e average strip ratio increases from phase 1 to phase 4. Since maximizing NPV implies postponing waste removal as much as possible, the increasing strip ratio with time is a verification of the rationality of the proposed method. Figure 4 is a 3D view of the four phase-pits of the optimized phase plan, where phase-4 pit is also the best ultimate pit. Figure 5 is a vertical cross section of the phasepits superimposed on the coal seams in the column model. e direction of phase expansion can be clearly seen from these figures. One can also see the rationality of the optimization results from the spatial relationship between the phase-pits and the coal seams ( Figure 5).
e developed software provides an option of keeping and outputting a specified number of best phase plan scenarios. For the case study, we found five other scenarios with virtually the same NPVs as the one given in Table 2, but with different numbers of phases, and/or phase-pits (including the ultimate pit). Since it is very difficult, if not impossible, to incorporate all relevant practical considerations in any optimization model and algorithm, these phase plan scenarios, which are equally good economically, provide valuable options for further evaluation with respect to certain practical considerations to arrive at the final phase plan.
We also analyzed the effects of certain input parameters, such as coal price, production costs, and their escalation rates, on the phase planning outcome for the case study, with the above phase plan as the base case for comparison. e planning outcome was found to be sensitive with respect to these parameters. When the coal price was lowered by 20%, the size of the optimum ultimate pit decreased by 30%, and the number of phases decreased from 4 to 3. Increasing the production costs by 20% had similar effects. Increasing the coal price (or lowering the production costs) had reverse effects on the planning outcome, as expected. Setting the annual escalation rates of coal price and production costs to 2.0% and 1.5%, respectively, resulted in a 7% larger ultimate pit, the same number of phases, but larger phase-pits. Based on the outcomes of these experiments, we suggest that the phases should be updated periodically (e.g., toward the end of each phase) in a real-life operation, as the relevant economic and technical conditions change  e optimization method presented herein and the developed software can be a handy tool for updating phase plans.

Conclusions
e phase planning method presented herein is specially designed for open pit coal mines with nearly horizontal and bedded coal seams. e major merits of the method include the following: it simultaneously optimizes the number of phases, the intermediate phase-pits, and the ultimate pit; it eliminates the gap problem in generating a series of nested pits; and the column model has a clear advantage over the commonly used block model in the accuracy of coal quantity computation with nearly horizontal and bedded coal seams. e method is capable of handling large real-life instances and produces rational results, as demonstrated by the case study. e method can also be used to analyze the effects of relevant input parameters on the phase planning outcome, providing useful scenarios for decision-making in strategic planning of open pit coal mines. e proposed method in its current form has two major shortcomings. One is that the total cost of mining a phase is averaged over the years of the phase's time span while, in actuality, the quantities of rock and unconsolidated material Mathematical Problems in Engineering mined each year fluctuate within a phase, causing fluctuations in annual cost. Another shortcoming concerns the transition from one phase to the next. Transition takes place sometime before a phase is completely mined out. During the transition period, the upper benches of the next phase are mined and some working benches may traverse the boundary between the two phase-pits. Phase transitions are not incorporated in our current formulation and should be dealt with in a more detailed scheduling process. Overcoming these shortcomings will be the focus of our future research on this topic.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that they have no conflicts of interest.