Abstract
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.
1. 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. The 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. These 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]. These 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–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]. The algorithm is capable of finding a series of nested, maximum-metal 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. The 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–14]. These 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]. Their 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. Their 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]. They 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. The authors applied the algorithm to a gold mine and compared the outcome with that from the conventional (deterministic) nested pit approach [17]. The 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 phase-pits, 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 phase-pits, and the ultimate pit.
2. 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. This will cause large errors in coal quantity calculation and subsequent economic evaluation, since whole blocks are classified as coal or waste. The 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. The 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. These 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. The resulting coal quantity is much more accurate than that based on whole blocks with a block model. The 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.
3. Generation of Least-Strip Ratio Pits
The 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. The 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 to generate a series of nested, least-SR pits. The basic logic of the algorithm is briefly described as follows.
Suppose that a series of n nested, least-SR pits, denoted as , is to be generated with a specified coal quantity increment of ΔQ, where P1 is the smallest pit and Pn the largest. Let Qi denote the coal quantity in Pi. The algorithm starts with the largest pit, Pn, 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. The pit slope angles are determined beforehand through rock stability studies and are inputs to the algorithm. The 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 Pn, 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 Pn. The remaining part of Pn after eliminating such a portion constitutes the next smaller pit, Pn−1, in the series, which has a coal quantity of Qn − ΔQ. To make Pn−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 Pn−1, the same cone eliminating process is repeated to generate the nest smaller pit, Pn−2. The process continues until the coal quantity in the remaining part is equal to or below the coal quantity, Q1, 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 Pi, 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 Pi−1 from pit Pi. The 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 Pi. The cone’s coal quantity is denoted as qc. If qc 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. Then, 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, qc, 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. Then, the cone is moved horizontally to another column inside pit Pi, as shown by the horizontal arrow in Figure 2, and the entire process for the previous column is repeated. The above process continues until all the columns inside pit Pi are traversed.
At the end of the above cone moving process, an array of J cones, each having a coal quantity qc ≤Δ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. The cones in the union are eliminated from pit Pi and the remaining part is pit Pi−1, whose coal quantity is about ΔQ smaller than Pi. 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.
The 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. Thus 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. Thus, the coal quantity increment between any consecutive pits in the generated pit series may be smaller than ΔQ, but cannot exceed 1.1ΔQ. The gap problem is, therefore, completely overcome.
The 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. The column size has a similar effect.
4. 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). The 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. The vertical axis represents the state variable P with each state corresponding to a pit in {P}n, depicted by a circle in Figure 3. The 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. That is why the starting (lowest) state of stage t corresponds to pit Pt 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⟶P2⟶P4⟶P6, 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 P2, P4, and P6, respectively; and the ultimate pit is P6. The path with the highest NPV is the optimal phase plan, which can be found by economically evaluating all the paths. The following is a DP formulation for finding the best path.
In general, suppose that pit Pi of stage t is being evaluated. Pi can be transited from those smaller-than-Pi pits of the preceding stage, t − 1. When pit Pi of stage t is transited from pit Pj of stage t − 1 (t −1 ≤ j ≤ i− 1), the quantities of coal, rock, and unconsolidated material mined in phase t, denoted as Qt,i (t − 1, j), Wt,i (t − 1, j), and Ut,i (t − 1, j), respectively, are calculated aswhere Qi, Wi, and Ui are the quantities of coal, rock, and unconsolidated material that can be mined from pit Pi, respectively. They 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 Vt,i (t − 1, j) and cost of Ct,i (t − 1, j) for phase t.where rp is the coal recovery rate of the processing plant; p is the coal price; and cm, cp, , and cu are the unit costs of coal mining, coal processing, rock mining, and unconsolidated material stripping, respectively.
Let yt,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. Then,where A is the annual coal mining capacity.
yt,i (t − 1, j) may not be an integer number of years, and let Lt,i (t − 1, j) be the integer part of it. The average annual revenue and cost for each of the Lt,i (t − 1, j) years, denoted by and ct,i (t − 1, j), respectively, are
The revenue and cost for the remaining decimal part, denoted by at,i (t − 1, j) and bt,i (t − 1, j), respectively, are
Following the transition, the cumulative time to arrive at pit Pi of stage t after finishing mining phase t, denoted by Yt,i (t − 1, j), iswhere Yt−1,j is the cumulative time to arrive at pit Pj of the preceding stage, t − 1, following the best path. Yt−1,j has been calculated in evaluating the states of the preceding stage, t − 1.
Therefore, when pit Pi of stage t is transited from pit Pj of stage t − 1, the cumulative NPV realized at pit Pi of stage t, after t phases of production, is given by NPVt,i (t − 1, j) aswhere NPVt−1,j is the cumulative NPV at pit Pj 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 Pi of stage t may be transited from all the smaller-than-Pi pits of the preceding stage, t − 1. Obviously, when pit Pi 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 Pi of stage t. The 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. The best transitions and the associated cumulative NPVs are obtained for all the pits of all stages. Then, find the pit that has the highest cumulative NPV of all pits of all stages. This pit is the best ultimate pit, and the stage at which it is found indicates the best number of phases. Then, 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. Thus, the number of phases, the phase-pits, and the ultimate pit are simultaneously optimized. This is a forward and open-ended DP formulation.
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. The topography of the coal field is nearly flat with a maximum relief of less than 20 m. The 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/m3. The deposit was modeled into a column model having around 39000 columns, each having a horizontal cross section of 30m × 30 m. The 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.
The 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. The coal is to be sold without processing. The time span of a single phase is generally more than 5 years to avoid complications associated with frequent transitions between phases. Therefore, in generating the least-SR pits, the coal quantity of the smallest pit, Q1, 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. The 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. This is a verification of the algorithm’s capability of completely overcoming the gap problem.
The phase plan was optimized with the generated pit series and the parameter values in Table 1 as inputs to the DP model. The 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. The 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 phase-pits superimposed on the coal seams in the column model. The 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).
(a)
(b)
(c)
(d)
The 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. The 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 over time. The optimization method presented herein and the developed software can be a handy tool for updating phase plans.
6. Conclusions
The phase planning method presented herein is specially designed for open pit coal mines with nearly horizontal and bedded coal seams. The 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. The method is capable of handling large real-life instances and produces rational results, as demonstrated by the case study. The 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.
The 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 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
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The authors acknowledge the National Natural Science Foundation of China (51674062, 51974060, and 51474049), National Science Foundation for Young Scientists of China (51604061), and Liaoning Province Key Research and Development Project (2019JH2/10300051), for their financial supports in the course of the research work presented in this paper.