Essays on Sequence Optimization in Block Cave Mining and Inventory Policies with Two Delivery Sizes by Anita Frances Parkinson B.Sc.Eng., Queen’s University, 1991 M.S.C.E.P., Massachusetts Institute of Technology, 1994 M.Sc.B.A., The University of British Columbia, 2002 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF Doctor of Philosophy in The Faculty of Graduate Studies (Business Administration) The University Of British Columbia August, 2012 c Anita Frances Parkinson 2012 Abstract Chapter 1 is an introductory chapter for this thesis work. It sets the scene by describing the motivation and industrial setting for each project. In Chapter 2, “Optimal Inventory Replenishment with Two Delivery Sizes”, we consider a periodic review inventory system where a retailer can order in multiples of a fixed quantity Q1 , or multiples of Q2 = 2Q1 , where the per unit material cost is less for ordering Q2 . We extend results of Veinott, and of Fangruo Chen, to show that an optimal replenishment policy has a reorder point R, as well as a second parameter controlling when the last order should be for Q1 instead of Q2 under a linear cost structure. In Chapter 3, “Sequence Optimization in Block Cave Mining” we investigate the use of integer programming models to aid the practitioner in the early planning stages of a Block Cave Mine. Given the footprint of the ore body divided into draw points or grid squares, sequence optimization determines which draw points to open in which period to meet the physical constraints of the mining process and maximize the total net present value of the mine. Traditionally done by trial and error by experts in the field, this is a first attempt to use modelling techniques to automate and optimize the process. We develop three integer programming models and discuss the challenges of formulating the problem in this framework. Two additional models are developed for comparison, one using the Column Generation technique and one using a greedy or myopic algorithm. All models are run on two data sets provided by our industrial partner, and the performance and results are compared. This work demonstrates that integer programming models can generate opening sequences but, like many “real life” problems, this one is complicated. ii Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . v List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . x Dedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xi 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 2 Optimal Inventory Replenishment with Two Delivery Sizes 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1.1 The Model and Overview of Results . . . . . . . . . . 2.2 Analytical Tools in the One Delivery Size Setting . . . . . . 2.2.1 Markov Decision Process . . . . . . . . . . . . . . . . 2.2.2 Finding the Optimal Average Cost Policy . . . . . . . 2.2.3 Previous Work . . . . . . . . . . . . . . . . . . . . . . 2.3 Extensions to Two Delivery Sizes . . . . . . . . . . . . . . . 2.3.1 Formal Model . . . . . . . . . . . . . . . . . . . . . . 2.3.2 Intuitive Starting Point . . . . . . . . . . . . . . . . . 2.3.3 Structure of Optimal Policy . . . . . . . . . . . . . . 2.3.4 Stationary Distribution of a Two Delivery Size Policy 2.3.5 Counterexamples to Optimality of Alpha-Policy . . . 2.3.6 An Algorithm for Computing Exact Solutions . . . . 2.3.7 Upper and Lower Bounds . . . . . . . . . . . . . . . . 2.4 Numerical Study . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 Further Extensions . . . . . . . . . . . . . . . . . . . . . . . 2.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 5 6 8 8 9 13 14 14 15 16 17 22 31 36 37 42 42 iii Table of Contents 2.7 Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . 43 3 Sequence Optimization in Block Cave Mining . . . . . . . . 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.1 Current Practice . . . . . . . . . . . . . . . . . . . . . 3.1.2 Tunnel Details . . . . . . . . . . . . . . . . . . . . . . 3.1.3 Cave Shape . . . . . . . . . . . . . . . . . . . . . . . 3.1.4 Other Considerations . . . . . . . . . . . . . . . . . . 3.1.5 Assumptions . . . . . . . . . . . . . . . . . . . . . . . 3.1.6 Previous work . . . . . . . . . . . . . . . . . . . . . . 3.2 Model Framework . . . . . . . . . . . . . . . . . . . . . . . . 3.2.1 Data . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.2 Decision Variables . . . . . . . . . . . . . . . . . . . . 3.2.3 Objective . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.4 Constraints . . . . . . . . . . . . . . . . . . . . . . . . 3.2.5 Unconstrained Sequence Optimization as a Draw Point Scheduling Model . . . . . . . . . . . . . . . . . . . . 3.3 Single Tunnel . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.1 Alternating Constraints Formulation . . . . . . . . . 3.3.2 Extended Formulation . . . . . . . . . . . . . . . . . 3.3.3 Computational Results . . . . . . . . . . . . . . . . . 3.4 Multiple Tunnels . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.1 Introduction/Overview . . . . . . . . . . . . . . . . . 3.4.2 Computation and Data Sets . . . . . . . . . . . . . . 3.4.3 Adapting Single Tunnel Formulations . . . . . . . . . 3.4.4 Malkin and Wolsey’s 2D Integral Formulation . . . . 3.4.5 Formulations Based on 4 Vertices . . . . . . . . . . . 3.4.6 Column Generation . . . . . . . . . . . . . . . . . . . 3.4.7 Greedy/Myopic Algorithm . . . . . . . . . . . . . . . 3.5 Application of Models to Test Data Sets . . . . . . . . . . . 3.5.1 Data Sets . . . . . . . . . . . . . . . . . . . . . . . . . 3.5.2 Models . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . 3.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 45 47 47 52 54 54 54 56 56 56 57 57 References 58 64 65 68 75 76 76 76 78 84 107 115 124 126 126 135 139 179 179 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 186 iv List of Tables 2.1 Possible Transitions from After-Ordering Inventory Position R+k . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Transition Probability Matrix for Only Ordering Large Delivery Sizes . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Transition Probability Matrix for Ordering a Small Delivery Size at Only One Inventory Position . . . . . . . . . . . . . . 2.4 Candidates for Optimal Policy: Counterexample 1 . . . . . . 2.5 Alpha Policy Candidates: Counterexample 1 . . . . . . . . . 2.6 Stationary Distribution and Average Cost of Candidate Policies: Counterexample 1 . . . . . . . . . . . . . . . . . . . . . 2.7 Candidate Policies: Counterexample 2 . . . . . . . . . . . . . 2.8 Stationary Distribution and Average Costs: Counterexample 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.9 Runs with Uniform Distribution . . . . . . . . . . . . . . . . 2.10 Runs with Normal Distribution . . . . . . . . . . . . . . . . . 2.11 Details of Policy Iteration . . . . . . . . . . . . . . . . . . . . 2.12 Runs with Two delivery sizes, Q2 6= 2Q1 . . . . . . . . . . . Table of Single Tunnel Comparison of Alternating Constraints Formulation and Extended Formulation . . . . . . . . . . . . 3.2 Size of Data Sets Provided for Model Development . . . . . . 3.3 Table of Running Times of Malkin Model from Data Sets P2 and P4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Partial Table of Running Times for Data Set P2 . . . . . . . 3.5 Comparison of Results of Wandering Axes Variation . . . . . 3.6 Available Data for Data Set RT1 . . . . . . . . . . . . . . . . 3.7 Maximum Active Draw Points for Data Set RT1 . . . . . . . 3.8 Available Data for Data Set RT2 . . . . . . . . . . . . . . . . 3.9 Maximum Active Draw Points for Data Set RT2 . . . . . . . 3.10 Number of Variables and Each Type of Constraints for Basic Model Runs on Data Set RT1 and RT2 . . . . . . . . . . . . 11 19 21 25 25 26 29 29 39 40 41 44 3.1 75 78 95 96 103 128 129 132 134 136 v List of Tables 3.11 Number of Variables and Each Type of Constraints for Malkin Model Runs on Data Set RT1 and RT2 . . . . . . . . . . . . 3.12 Number of Variables and Each Type of Constraints for 2Cone Model Runs on Data Set RT1 and RT2 . . . . . . . . . . . . 3.13 Results for Data Set RT1 . . . . . . . . . . . . . . . . . . . . 3.14 Cave Development Results for Data Set RT1 . . . . . . . . . 3.15 Solution Times: . . . . . . . . . . . . . . . . . . . . . . . . . . 3.16 Results of Sensitivity Runs: . . . . . . . . . . . . . . . . . . . 3.17 Cave Sizes in Sensitivity Runs: . . . . . . . . . . . . . . . . . 3.18 Results for Data Set RT2: . . . . . . . . . . . . . . . . . . . 3.19 Cave Development Results for Data Set RT2: . . . . . . . . . 3.20 Results from Selected Runs of Malkin Model on RT2 . . . . . 3.21 Integer Programming Formulation Comparison . . . . . . . . 137 138 139 140 140 161 161 168 170 173 180 vi List of Figures 2.1 2.2 2.3 3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 3.9 3.10 3.11 3.12 3.13 3.14 3.15 3.16 3.17 Two Examples of Set of Q Q=6 . . . . . . . . . . . . Graphical Representation of ple 1 . . . . . . . . . . . . . Graphical Representation of ple 2 . . . . . . . . . . . . . Integer Minimizers of G(·) for . . . . . . . . . . . . . . . . . . . Optimal Policy: Counterexam. . . . . . . . . . . . . . . . . . . Optimal Policy: Counterexam. . . . . . . . . . . . . . . . . . . A Cross Section of the Underground Workings at the Dutoitspan Mine . . . . . . . . . . . . . . . . . . . . . . . . . . . Cross Section Picture of Block Cave Mine . . . . . . . . . . . Plan of Extraction Level of Block Cave Mining Operation . . Steps in Undercut Development. Adapted from Resolution Copper Mining (2009) . . . . . . . . . . . . . . . . . . . . . . Example of Conventional Undercut Layout Showing Blasting Pattern. Barber (2000) . . . . . . . . . . . . . . . . . . . . . . Three Examples of Single Tunnel Opening Sequences . . . . . Plan of Undercut Level of Block Cave Mining Operation . . . Typical Diamond Pattern of Opened Draw Points . . . . . . . Network Representation of Single Tunnel in One Period . . . Simple Network . . . . . . . . . . . . . . . . . . . . . . . . . . Footprints of Data Sets Provided for Model Development . . Example of Tunnel Numbering . . . . . . . . . . . . . . . . . Opening Patterns from Data Set P1 using Extended Formulation Within-Tunnels and Across-Tunnels . . . . . . . . . . . Opening Patterns from Data Set P2 using Extended Formulation Within-Tunnels and Across-Tunnels . . . . . . . . . . . Example of Centre Point 4,21 . . . . . . . . . . . . . . . . . . Example of Centre Point 4,21 and Direction of the Column Connected Constraint((3.24)-(3.25)) . . . . . . . . . . . . . . Example of Centre Point 4,21 and Direction of the Row Connected Constraint((3.26)-(3.29)) . . . . . . . . . . . . . . . . . 12 27 30 46 48 49 50 50 51 52 53 70 71 77 79 82 83 86 88 89 vii List of Figures 3.18 3.19 3.20 3.21 3.22 3.23 3.24 3.25 3.26 3.27 3.28 3.29 3.30 3.31 3.32 3.33 3.34 3.35 3.36 3.37 3.38 3.39 3.40 3.41 3.42 3.43 3.44 3.45 3.46 3.47 3.48 3.49 Example of Centre Point 4,21 and Arrows for All Constraints Data Set P2, Period 1 Opening . . . . . . . . . . . . . . . . . Data Set P4, Centre Point (14,65) . . . . . . . . . . . . . . . Horizontal Axis Bent Around Hole in Footprint, Centre Point (14,65) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Histogram of Total NPV . . . . . . . . . . . . . . . . . . . . . Surface Plot of Total NPV for Data Set P2 . . . . . . . . . . Projections of the Surface Plots for Draw Point Values of wpii , wi exp(−λpi ) and Total NPV for Data Set P2 . . . . . . . . . . . Table of Neighbour Values for Data Set P2 . . . . . . . . . . Results from the Wandering Axes Model Variation on Data Set P2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Results from the Malkin Model Centered at Tunnel 6, Draw Point 12 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Same Figures Side by Side . . . . . . . . . . . . . . . . . . . . Example of Diamonds with Vertices Having Only One Open Neighbour . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Example of Diamond with a (Top) Vertex Having Two Open Neighbours . . . . . . . . . . . . . . . . . . . . . . . . . . . . Three Examples of Pairs of Cones and Their Intersections . . An Example of the Extended Footprint for the 2Cone Model Data Set RT1 . . . . . . . . . . . . . . . . . . . . . . . . . . . Histogram of Draw Point Durations at 120 tons/day for Data Set RT1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Distribution of Draw Point Values over the Mining Footprint for Data Set RT1 . . . . . . . . . . . . . . . . . . . . . . . . . Data Set RT2 . . . . . . . . . . . . . . . . . . . . . . . . . . . Histogram of Durations for Data Set RT2 . . . . . . . . . . . Results of Basic Model on Data Set RT1 - Period 1 . . . . . . Results of Basic Model on Data Set RT1 - Period 2 . . . . . . Results of Basic Model on Data Set RT1 - Period 3 . . . . . . Results of Basic Model on Data Set RT1 - Period 4 . . . . . . Results of Basic Model on Data Set RT1 - Period 10 . . . . . Best Results of Malkin Model for Data Set RT1 - Period 1 . . Best Results of Malkin Model for Data Set RT1 - Period 2 . . Best Results of Malkin Model for Data Set RT1 - Period 3 . . Best Results of Malkin Model for Data Set RT1 - Period 4 . . Best Results of Malkin Model for Data Set RT1 - Period 10 . Results of 2Cone Model for Data Set RT1 - Period 1 . . . . . Results of 2Cone Model for Data Set RT1 - Period 2 . . . . . 90 94 97 98 99 100 101 103 104 105 105 107 109 110 112 127 129 130 131 133 141 142 142 143 143 144 145 145 146 146 147 148 viii List of Figures 3.50 3.51 3.52 3.53 3.54 3.55 3.56 3.57 3.58 3.59 3.60 3.61 3.62 3.63 3.64 3.65 3.66 3.67 3.68 3.69 3.70 3.71 3.72 3.73 3.74 3.75 3.76 3.77 3.78 3.79 Results of 2Cone Model for Data Set RT1 - Period 3 . . . . . Results of 2Cone Model for Data Set RT1 - Period 4 . . . . . Results of 2Cone Model for Data Set RT1 - Period 10 . . . . Results of ColGen Model for Data Set RT1 - Period 1 . . . . Results of ColGen Model for Data Set RT1 - Period 2 . . . . Results of ColGen Model for Data Set RT1 - Period 3 . . . . Results of ColGen Model for Data Set RT1 - Period 4 . . . . Results of ColGen Model for Data Set RT1 - Period 10 . . . . Results of Greedy Model for Data Set RT1 - Period 1 . . . . Results of Greedy Model for Data Set RT1 - Period 2 . . . . Results of Greedy Model for Data Set RT1 - Period 3 . . . . Results of Greedy Model for Data Set RT1 - Period 4 . . . . Results of Greedy Model for Data Set RT1 - Period 10 . . . . Opening Patterns for First Two Periods for Basic Model . . . Opening Patterns for First Two Periods for Malkin Model . . Opening Patterns for First Two Periods for 2Cone Model . . Opening Patterns for First Two Periods for ColGen Model . . Opening Patterns for First Two Periods for Greedy Model . . Published Sequence Optimization for Sample Data Set . . . . Results of 2Cone Model Starting at Tunnel 5 . . . . . . . . . Results of 2Cone Model Starting on Edge Tunnel . . . . . . . Results of 2Cone Model Starting on Any Edge . . . . . . . . Results of 2Cone Model Minimized . . . . . . . . . . . . . . . Results of Minimum Malkin . . . . . . . . . . . . . . . . . . . Overlay of Published Results and First 3 periods of Generated Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Results of 2Cone Model Starting at Centre . . . . . . . . . . Results from Basic Model on Data Set RT2 . . . . . . . . . . Results from Malkin Model (29,80) for Data Set RT2 . . . . . Results from Malkin Model (39,50) . . . . . . . . . . . . . . . Results from Greedy Algorithm on Data Set RT2 . . . . . . . 148 149 149 150 151 151 152 152 153 154 154 155 155 156 157 157 158 158 160 162 162 163 163 164 165 167 172 175 176 178 ix Acknowledgements Thank you to my thesis advisors Tom McCormick and Maurice Queyranne who guided and stuck with me through all the ups, downs and gaps in my thesis process. A big thank you to Tony Diering of Gemcom who was willing to take on this project with me and taught me all I know about block cave mining. Thanks also to Brian Graham, who taught me many, many things and kept me smiling. As is often said, it is not the arrival that is important, but how you get there. x Dedication This thesis is dedicated to all my fellow students who have struggled with the Ph.D. experience, and to all the family, friends and humourists who help us through it. If you read this let me know firstname.lastname at gmail “I see my light come shining From the west unto the east Any day now, any day now I shall be released” Bob Dylan xi Chapter 1 Introduction This thesis contains two chapters of work that are essentially unrelated to each other. In choosing to work on both of them I was motivated by real industrial problems. In a previous life, I worked as a chemical engineer, the perfect blend of theory in the office, and actual implementation in the chemical plant. There is nothing like staring at a three storey tall evaporator to realize that getting the temperature value that is essential for the computer model, will be problematic. Theory and models ARE great, but the art is making it work, or extracting the understanding that you can use to make it work. When I returned to school to study Operations Research I chose a master’s program that was focused on an industrial project. With this background, it is not surprising that I sought out problems with an industrial basis for this thesis. Behind the inventory policy work (Chapter 2 Optimal Inventory Replenishment with Two Delivery Sizes), is a shipping problem. A company produces a product in the southern hemisphere and transports it by using a fleet of tanker ships to customers in Europe and North America; how should they schedule their tankers to meet expected demand? The classic Economic Order Quantity (EOQ) model sets a target inventory position and prescribes ordering (or delivering) enough to meet that target in each period. The target trades off between inventory holding costs and fixed delivery charges. Unless demand is constant, this results in different order quantities in each time period. In the setting of a fleet of tanker ships, the size of the available ships may not match the desired order quantity. Also, it would be expected that the delivery costs would vary with the size of the ship, typically showing economies of scale. The chapter investigates an ordering policy for a very simple fleet — one with an unlimited number of ships of two sizes, one twice the size of the other. The second project (Chapter 3 Sequence Optimization in Block Cave Mining) addresses a problem from the mining industry; given that you are going to mine a given ore body, where do you start and what direction do you proceed? The current practice is a trial and error process where options are generated using expert knowledge and evaluated by a set of criteria. 1 Chapter 1. Introduction This is a classic opportunity where the tools of operations research can be applied to try to capture some of the experience and art of decision making and proved to be an inviting project. The planning of open pit mines has long been aided by integer programming models, but these models do not apply to block cave mines. In block cave mining, the ore body is approached from underneath. Blasting is done to start the break up of the rock in a given area called a draw point. When the blasted rock is removed, a cave is created. Over time the rock continues to break off from the top of the cave and is removed for processing. Sequence optimization plans where to do the initial blasting in each period to create a single cave that matches processing capabilities, and maximizes the total Net Present Value of the mine. By assuming a fixed removal rate from each draw point, the three dimensional ore body is reduced to a two dimensional mining footprint. In this work we develop three integer programming models. The main challenge that emerges is the constraint that the open draw points form a single contiguous cave. This proved difficult in the integer programming framework, though two of the models were able to meet it. In addition, the column generation technique is used to move the problem of feasible cave generation out of the integer programming framework. All of the models share three basic constraints. Two are straightforward to frame in the integer programming framework. The Start Once constraint ensures that each draw point is opened once and only once, and the Global Capacity constraint ensures the number of active draw points does not exceed the downstream processing capacity. The third constraint, that the opened draw points must form a single, contiguous group, or cave, is the source of the model variations. In the first model, the Basic model, contiguity is enforced along a two dimensional grid. Two formulations to achieve this are discussed and one is selected for the model. Unfortunately, two dimensions are not sufficient to guarantee a single cave, as is shown on the sample data. This model does prove useful on our sample data as it provides an upper bound on the Total Net Present Value of the mine. The model runs to completion in a reasonable time period on all the sample data sets. The Malkin model comes from Malkin and Wolsey, whose contributions are much appreciated. Their model also incorporates the preference of planning experts for caves with a diamond shape, which is based on the rock dynamics. A centre point is chosen, which establishes two perpendicular axes in a two dimensional grid. Contiguity is achieved by only allowing a draw point to open if those draw points between it and each axis are open. 2 Chapter 1. Introduction This model runs well on all sample data sets, producing the single, contiguous cave in each period as desired. The challenge comes in choosing the centre point. For the smaller sample data set it is feasible to have a run with each draw point as the centre point. This is not feasible for the larger sample data set. Achieving a diamond shape cave is also the driving force behind the 2Cone model in which a cave is defined by the intersection of two cone or triangle shaped groups of draw points. Each cone is defined by a vertex, and grows by two draw points on each subsequent row of a two-dimensional grid. The location of each vertex is a decision variable so is chosen by the model to maximize the Total Net Present Value of the mine. This model does not complete on all the sample data sets, but when it does, it produces a single contiguous cave in each period. The difficulties in achieving a single, contiguous cave in each period within the integer programming framework gave rise to the use of the Column Generation technique to create the ColGen model. Briefly, this technique iterates between generating a set of feasible caves for each period and using integer programming techniques to select a cave in each period that maximizes the objective and meets the constraints. The appeal of this technique is that the feasible caves can be generated outside of the integer programming framework. Unfortunately, the technique fails to generate valuable opening sequences. Additional algorithms are generated which ensure that each feasible cave added by the technique has caves in preceding and following periods that overlap to form a feasible sequence. This model did not complete on all data sets, but did produce a feasible solution on the smaller sets. In this chapter we compare the performance of the models on two data sets. All models perform well on the smaller data set of 236 draw points. The Basic model returns a sequence with the highest value, though two caves are generated in the early periods. The Malkin model returns a value that is 98% of the Basic, and the 2Cone and ColGen models come in at 95% of the the Basic model value. On this data set, the Basic model completes in under 10 seconds while the others take over an hour. We suggest that this speed shows a usefulness of the Basic model in setting an upper bound on the value, and indicating high value draw points for initial periods. The second data set is much larger at 6510 draw points. On this data set, the ColGen and 2Cone models failed to complete in under a week. After 19 hours the Basic model returned a sequence which again had multiple caves in the early periods. It was not feasible to run the Malkin model on all draw points, but the highest valued of those selected reached 98.8% of the Basic model value. 3 Chapter 1. Introduction The selection of individual Malkin models runs varied from 1 to 9 hours to complete. While these running times may seem prohibitively long to many researchers, in the workplace it may be feasible to run a model over night, or over a weekend if a single run produces useful results. Again, we suggest that the Basic model is useful in producing an upper bound against which the selection of Malkin runs can be evaluated. 4 Chapter 2 Optimal Inventory Replenishment with Two Delivery Sizes 2.1 Introduction Consider a retailer facing uncertain demand who periodically replenishes her stock from a manufacturer (or distributor, or . . . ). Classic models assume that the quantity ordered can be any number of units. For simple version of such models, often one can find replenishment policies that are optimal, in the sense of minimizing long-run average cost. In the absence of fixed costs for ordering, such policies are often of the form of having an order “up to level”, or reorder point, R. This means that if the retailer observes that her stock is below R at an ordering epoch, she orders up to R; otherwise, she orders nothing. In reality, the order size is often limited to, e.g., multiples of a full truckload. Suppose that a full truckload is Q units, so that all orders must be integral multiples of Q. In this case, Veinott (1965) showed that, for an appropriate R, a variant of an order up to level policy is optimal: In each period, the retailer should order the smallest multiple of Q necessary to bring her inventory position to at least R. Later, Chen (2000) extended this result to a multi-echelon setting. In these cases the per-unit cost of items is always the same, so material cost can be ignored. But in many practical situations, there is a variety of order sizes possible. For example, there could be both large and small trucks. In a consulting project we did with a bulk chemical supplier, there were several different sizes of ocean tankers available for shipping their product. This naturally raises the question of what optimal ordering policies would be in such cases. This paper considers a stylized version of such a problem. We assume that there are two delivery sizes available, Q1 (small truck), and Q2 (big 5 Chapter 2. Optimal Inventory Replenishment with Two Delivery Sizes truck), with Q2 = 2Q1 . We defend this simple model by pointing out that even this simplest generalization turns out to be quite hard to solve. Denote the per-unit cost of one item on a big truck as c2 , and on a small truck as c1 , so that the total cost of ordering a big truckload is c2 Q2 . We make the reasonable assumption that c2 ≤ c1 , due to quantity discounts and/or economies of scale in transportation. This implies that our material ordering cost is no longer independent of our replenishment policy, since it depends on the mix of Q2 and Q1 that we order. Since it is cheaper on a per-unit basis to order big truckloads, it will never pay to order more than one small truckload. This means that if we order a small truckload, we effectively pay a penalty of (c1 − c2 )Q1 for that delivery. Note that if c1 < c2 there is never any reason to order the large delivery size and this becomes Chen and Veinotts’s setting. It is natural to conjecture that some more elaborate version of a reorderpoint policy will be optimal here. Thus there should be some value R such that the retailer should order enough Q2 ’s to bring her inventory just below R, and then choose whether her last order should be a Q2 or a Q1 . This leads to the basic tradeoff in this model: If the last order is a Q2 , then we save on material cost, but at the expense of having a larger inventory, and hence larger holding costs. On the other hand, if the last order is a Q1 , then we pay more for the material, but pay lower holding costs as a result of having less inventory. Our analysis extends that of Chen (2000). 2.1.1 The Model and Overview of Results We assume that our retailer operates a single location using periodic review. Demand during each period is an integer number of units, and follows the stationary distribution function f . Demands in different periods are independent and identically distributed. If demand during a period is greater than inventory, then we backlog the unserved demand (at a cost). At the beginning of the period, the retailer assesses her inventory position and places an order. For simplicity we assume zero lead time, and so the order arrives instantly. Over the period demand is realized and at the end of the period, one-period holding or backorder costs, independent of unit purchase costs, are assessed. We assume that the selling price of the good is fixed, and does not vary with demand. As a result of this assumption, price does not enter the model as all demand is met, and price does not affect demand. The retailer can order any quantity of the form aQ2 + bQ1 , where a is 6 Chapter 2. Optimal Inventory Replenishment with Two Delivery Sizes a non-negative integer and where we can assess, without loss of generality (see Section 2.1), that b is zero or one. We assume an infinite horizon and define an “optimal” policy as one that minimizes long-run average costs. We follow Chen in deriving our results for general cost functions, assuming only some convexity properties of the one-period expected cost. As with Chen, we will be able to show that, e.g., linear holding and backorder costs lead to cost functions satisfying these properties. Let’s make an intuitive argument for what an optimal policy should look like. It certainly should have a reorder point R. As already mentioned, the retailer should order enough Q2 ’s to bring her inventory position into the interval (R − Q2 , R], and then decide on her last order. If the all but last order brings inventory into (R − Q2 , R − Q1 ], then the last order must be for a large truck to bring the inventory above R. Otherwise inventory is (temporarily) in (R − Q1 , R] and the last order can be either Q1 or Q2 . If we define oi as the amount to order from inventory position i, then a policy can be described as [oR−Q1 +1 , oR−Q1 +2 , . . . , oR ] with each oi = Q1 or Q2 . Let G(y) denote the expected one period cost given the inventory position is y after ordering. Inventory position is the sum of items on hand plus items on order. Chen proves that with one delivery size Q available and reorder point R, the stationary distribution of the Markov chain on inventoryP position is uniform on (R, R + Q] and hence the long-run average R+Q 1 cost is Q y=R+1 G(y). In our case the stationary distribution of inventory position under a two delivery size policy can be more complicated (see Section 2.3.4), so we must analyze the problem differently than Chen. We first describe the single delivery size case as a Markov Decision Process and show that the optimal policy has a set of recurrent or target states containing the Q smallest integer minimizers of G(y). Next we use similar techniques to show that the set of recurrent states for the two delivery size policy will include the integer minimizer of G(y), plus the minimizer’s Q1 neighbours which are the next ordered integer minimizers. The set of recurrent states may also include the next Q1 integer minimizers of G(y). We further provide some guidance on how R and the optimal policy can be computed in practice. Two approximations are provided as well as an algorithm for finding the optimal policy based on the Markov Decision process technique of Policy Iteration. It turns out, in the examples run, that the reorder points associated with the policies of always ordering just Q2 ’s, or always ordering just Q1 ’s, provide quite good starting approximations for 7 Chapter 2. Optimal Inventory Replenishment with Two Delivery Sizes our optimal R. 2.2 2.2.1 Analytical Tools in the One Delivery Size Setting Markov Decision Process We will begin the analysis by reviewing the standard model of Veinott (1965) and the work of Chen (2000) to set the stage for our work. We will consider the setting as a periodic review inventory model at a single location. At the beginning of the period, inventory position s is observed and an order is placed, bringing the inventory position to y. To fix ideas, let us consider the simple case of one delivery size Q available for ordering. Demand in each period D, is assumed to be independent and identically distributed taking only integer values. Unmet demand in any period is backlogged. At the end of the period a holding charge is assessed on any excess inventory which is then held to the next period, and a penalty is paid for any backlogged demand. Let G(y) be the expected one period holding and backorder cost given that the inventory position immediately after ordering is y. Since all delivery sizes have the same cost, and unmet demand is backlogged, item purchasing costs can be ignored. The optimal policy will minimize the long run average cost under an infinite time horizon. In Chen’s work, he can prove that G(y) is quasiconvex in y, but we assume that it is to make our analysis work. See Newman (1987) for a definition of a quasiconvex function. Recognizing that inventory position is a Markov chain we follow Puterman’s (1994) Markov decision process formulation by identifying the following: Decision Epochs The beginning of the period when the ordering decision must be made States The inventory position at the time the ordering decision is made. Define the set of inventory positions at the time the ordering decision is made by S, with lower bound s and upper bound s, S = {s, s + 1, . . . , −1, 0, 1, . . . , s − 1, s}. Under reasonable demand distribution, i.e., one that is bounded below by zero, and above by some finite value, and a reasonable ordering policy, i.e., one that keeps inventory position bounded, the inventory position at the beginning of the period belongs to a finite set. 8 Chapter 2. Optimal Inventory Replenishment with Two Delivery Sizes Actions Order some integer multiple a of delivery size Q. Define the set of actions for a given state s (the inventory position at the time of ordering), As , as As = {0, 1, 2, . . . , a}. Following the assumptions that lead to a bounded state space, we assume there is an upper limit, a, to the size of the order placed. Expected Rewards At the end of the period, after demand has been realized, holding and backorder costs are assessed. Define the cost for action a taken from inventory position s as r(s, a) = G(s + aQ) Transition Probabilities Once the order has been made, the transition probabilities depend on the demand distribution. Denote the probability of demand equal to d as P {D = d} and define the transition probability to inventory position j (before the next order is made), after taking action a from inventory position s as P {D = s + aQ − j}, j ≤ s + aQ; p(j|s, a) = 0, j > s + aQ. 2.2.2 Finding the Optimal Average Cost Policy We seek an ordering policy π ∗ , consisting of actions as , one for each s ∈ S, which minimizes the long run average cost. For a model such as ours with stationary and bounded rewards, stationary transition probabilities and finite, discrete state space S, the resulting infinite-horizon Markov chain {Yt : t = 1, 2, . . .}, has average reward or gain of a policy π given as N 1 π X g (y) = lim Ey { r(Yt )} N →∞ N π t=1 The policy’s transition probability matrix Pπ = Pπ (s, j) = (p(j|s, as )) has a limiting matrix Pπ∗ . Denote the stationary distribution P as q, the soT T lution to the system of equations q = q Pπ subject to j∈S q(s) = 1. It is easy to show that Pπ∗ = eq T Puterman (1994) where e is the unit vector. Under the restrictions of our model, the average reward converges to g π (s) = Pπ∗ r(s). The average cost in this form can be described as a weighted average of the expected costs of the states that are reached. The weights are the long run probability of spending a period in the state. Clearly, the average 9 Chapter 2. Optimal Inventory Replenishment with Two Delivery Sizes cost is minimized under a policy that orders to the inventory positions that minimize G(y), the expected one period costs. Let y ∗ be an integer that minimizes G(y). When G(y) is quasiconvex, then the optimal policy is to order to y ∗ in each period. If Q = 1 this would be possible from all s ∈ S, and results in the class of order-up-to policies. For Q > 1 it is not always possible to reach y ∗ . From an inventory position s < y ∗ , if it is not possible to reach y ∗ by ordering multiples of Q, i.e., (y ∗ − s) mod Q 6= 0, then cost is minimized by ordering to arg min{G(s + ∗ ∗ b y Q−s cQ), G(s + d y Q−s eQ)}. That is, ordering either to the closest position above or below y ∗ , whichever has a lower expected one period cost. From inventory s ≥ y ∗ , no orders should be placed. Thus, from any inventory position s < y ∗ ordering enough deliveries to reach y ∈ (y ∗ − Q, y ∗ ], then for each inventory position y ∈ [y ∗ − Q, y ∗ ), ordering an additional delivery size if G(y + Q) < G(y), is the optimal policy. For ease of future notation, let y1 = y ∗ be the integer that minimizes G(y) . Similarly, let y2 be the ‘second’ integer minimizer of G(y), i.e., y2 ∈ Z, G(y1 ) ≤ G(y2 ) ≤ G(y) ∀ y ∈ Z, y 6= y1 . Since G(y) is quasiconvex, y2 = y1 − 1 or y2 = y1 + 1, i.e., y2 is a neighbour of y1 Definition: Let y1 , y2 , . . . , yQ be such that G(y1 ) ≤ G(y2 ) ≤ . . . ≤ G(yQ ) ≤ G(y) ∀y ∈ Z 6= yi i = 1, . . . , Q. Define MG,Q ≡ {yi }, i = 1, . . . , Q, the ordered set of Q integer minimizers of G(y). If G(y) is quasiconvex then yi = maxj<i {yj } + 1 or yi = minj<i {yj } − 1 If G(y) is quasiconvex then the optimal policy is to order into the band of target inventory positions (R, R + Q] containing MG,Q , the Q integer minimizers of G(y). The bottom of the band is R = minMG,Q −1. Under our optimal policy structure, the inventory position after ordering, y, forms a Markov chain {Yt : t = 1, 2, . . .} over a finite set of inventory positions Y . As shown in the previous paragraph, Y contains Q inventory positions which form a contiguous band. In order to evaluate the average cost, we need the transition probability matrix Pπ . As before, the reward function r(y) = G(y). The elements of the transition probability matrix pij are the probabilities of moving from Yt = i to Yt+1 = j under the policy π, and are solely a function of the demand faced between orders. Before constructing Pπ for our optimal policy structure we illustrate the possible transitions that can be made from target after-ordering inventory position R + k for integer k, 1 ≤ k ≤ Q. Each line of Table 2.1 below shows the demand faced, the inventory position (IP) after the specified demand has be realized, the order placed to return to the set of target states Y , and the end state that is reached. 10 Chapter 2. Optimal Inventory Replenishment with Two Delivery Sizes Starting State R+k Demand 0 1 .. . IP after demand R+k R+k−1 .. . Order 0 0 .. . End State R+k R+k−1 .. . k−1 k k+1 .. . R+1 R R−1 .. . 0 Q Q .. . R+1 R+Q R+Q−1 .. . Q−1 Q Q+1 .. . R+k−Q+1 R+k−Q R+k−Q−1 .. . Q Q Q .. . R+k+1 R+k R+k−1 .. . Q+k−1 Q+k Q+k+1 .. . R−Q+1 R−Q R−Q−1 .. . Q 2Q 2Q .. . R+1 R+Q R+Q−1 .. . Table 2.1: Possible Transitions from After-Ordering Inventory Position R+k From this it is easy to construct the probability transition matrix To State → From State ↓ R+1 ∞ X R+1 P {D = jQ} ... j=0 Pπ = R+Q−1 ... ∞ X P {D = jQ + 2} j=0 R+Q ∞ X P {D = jQ + 1} j=0 . . . R+Q−1 R+Q ∞ X j=0 ∞ X j=0 P {D = jQ + Q − 2} ... ∞ X P {D = jQ} j=0 P {D = jQ + Q − 1} ... ∞ X P {D = jQ + 1} j=0 ∞ X P {D = jQ + 2} j=0 ∞ X P {D = jQ} j=0 Inspection reveals that each row of the probability matrix has the same entries as the row above, with the values each shifted one column to the right. In addition to Pπ being stochastic, the columns of Pπ sum to one. As a result of this symmetry, the solution to q T = q T Pπ is q(R + Q) = q(R + Q − 1) = . . . = q(R + 1) = 1 Q 11 Chapter 2. Optimal Inventory Replenishment with Two Delivery Sizes and the long run average cost is Q 1 X G(R + 1 + i) Q i=1 as shown by Chen (2000). Again, P we note that if G(·) is quasiconvex, the long run average cost g = Q1 Q i=1 G(R + 1 + i) is minimized by MG,Q . The sequential order of the integer minimizers will vary with G(·), but the quasiconvexity of G(·) ensures that they form a contiguous band of inventory positions as shown in Figure 2.1. Figure 2.1: Two Examples of Set of Q Integer Minimizers of G(·) for Q = 6 These figures also illustrate that for Q = i, Q0 = j with i < j the target inventory positions for the setting in which multiples of delivery size Q are ordered are a subset of the target inventory positions for the setting in which multiples of delivery size Q0 are ordered. This is because MG,Q ⊂ MG,Q0 . In this section we have shown that the optimal policy under delivery size ordering of size Q is to order into the band (R, R + Q], commonly known as the (R, nQ) policy. The long run average cost of this policy is 1 PQ i=1 G(R + 1 + i). This result is a replication of the earlier work of Q Veinott and Chen which is described below. We also demonstrated that the band (R, R + Q] contains the set of Q integer minimizers of G(x), the one period holding and backorder costs. 12 Chapter 2. Optimal Inventory Replenishment with Two Delivery Sizes 2.2.3 Previous Work Two approaches to this single delivery size problem have been presented in the literature. In 1965 Veinott showed the optimality of the (R, nQ) policy, which he called a (k, Q) policy, for the n-period and infinite period models. In his setting, at the beginning of a period inventory is reviewed and an order, consisting of integer multiples of delivery size Q, placed which will be delivered λ periods later. Over the period previous orders are delivered and new demand arrives. At the end of the period holding and penalty costs are assessed on the inventory on hand. Ordering costs of unit cost c per item, are assessed when the order is made, but are discounted to the end of the period with one period discount factor α, 0 ≤ α ≤ 1. Under the (k, Q) policy, if the inventory position at the beginning of the period is less than k, order the smallest integral multiple of Q that will bring the inventory position to at least k, otherwise, do not order. Given the inventory position after ordering is y the cost function is the sum of the ordering cost, and L(y), the expected future holding and penalty costs. G(y) = (1 − α)cy + L(y). G(y) is assumed to be unimodal and minimized at y. The policy is optimal if it minimizes the expected discounted cost. The minimizer k is chosen from all values l such that l ≤ y ≤ l + Q. The policy is first proved optimal for the n-period model by comparing the period-by-period inventory position under the (k, Q) policy to that under any other policy. Since unmet demand is backlogged, the difference in inventory positions after ordering under the two policies will be an integral multiple of Q in each period. In any given period i, the after-ordering inventory position under the (k, Q) policy, yi0 , will satisfy k ≤ yi0 ≤ k + Q and the after-ordering inventory position under the alternate policy yi will satisfy either yi = yi0 or |yi − yi0 | ≥ Q. In both cases G(yi ) ≥ G(yi0 ). Under the (k, Q) policy, yi0 ≥ k + Q can only occur if y10 ≥ k + Q and by period i there has not been enough demand to bring yi0 below k + Q. In this situation, yi ≥ yi0 i = 1 . . . j where j is the first period where yi0 < k + Q, hence G(yi ) ≥ G(yi0 ) for these periods. Since this policy minimizes G(yi ) it minimizes the expected cost in the n-period model and also the infinite period model. Chen looks at the infinite horizon, non-discounted problem where G(y) is the expected cost (holding and backorder) in a period given the afterordering inventory position is y. Chen shows that if band [R + 1, R + Q] contains the integer Pminimizers of G(y + nQ) for any y and integer n, then this R minimizes Q x=1 G(y + x). Since the demands in each period are independent, the after-ordering inventory position under the (R, nQ) policy 13 Chapter 2. Optimal Inventory Replenishment with Two Delivery Sizes will follow a Markov chain. Under some mild conditions on the demand distribution, the steady-state distribution of the P Markov chain is uniform. Thus the long run average single period cost is Q1 R+Q y=R+1 G(y). This result P is more general than Veinott’s as Chen assumes that Q x=1 G(y + x) is quasiconvex whereas Veinott assumes G(·) is quasiconvex. Also, Chen extends to a multi-echelon setting. 2.3 Extensions to Two Delivery Sizes In this section, the setting is extended to one where two delivery sizes are available for ordering. We show that this complicates the policy determination over the single delivery size, illustrate a method of determining the optimal policy and finally show computationally easy upper and lower bounds on the long run average cost. 2.3.1 Formal Model We start with the single delivery size model and extend it to multiple delivery sizes. We remind the reader of the formal model in which at the beginning of the period, inventory position s is determined and an order is placed. Demand is realized over the period with any unmet demand backlogged to a later period. Assume that demands in each period are independent and identically distributed taking only integer values, and that demand is bounded. At the end of the period a penalty is paid for any backlogged demand and any inventory held over to the next period. Define G(y) as the single period expected cost given the inventory position after ordering is y. Now we extend the model to allow for orders made up from two delivery sizes. Consider the specific case of a small delivery size Q1 and a larger delivery size Q2 = 2Q1 . Any order can consist of an integer multiple of each delivery size. There are some economies of scale so the unit cost of the small delivery size is higher than that of the larger delivery size. If ci is the unit cost of delivery size Qi , then c1 > c2 . In this case, it is never economical to order more than one of the small delivery sizes and orders will be of the form aQ2 + bQ1 for integers a ≥ 0, b = 0, 1. Since all demand is met in the model, all reasonable policies will order the same number of items over the long run. We can consider the material or purchase cost at the lower unit cost c2 a sunk cost and we can account for differences in unit cost by charging a cost of Q1 (c1 − c2 ) for each small delivery size ordered. 14 Chapter 2. Optimal Inventory Replenishment with Two Delivery Sizes These changes result in the following changes to the definition of the Markov decision process formulation given in Section 2.2.1 above: Actions Place an order as Q2 + bs Q1 As = {(as , bs )}, as ∈ {0, 1, 2, . . . , a}, bs ∈ {0, 1}. Expected Rewards At the end of the period, after demand has been realized, holding and backorder costs are assessed. r(s, (as , bs )) = G(s + as Q2 + bs Q1 ) + bs Q1 (c1 − c2 ) Transition Probabilities p(j|s, (as , bs )) = P rob{D = s + as Q2 + bs Q1 − j}, j ≤ s + as Q2 + bs Q1 0, j > s + as Q2 + bs Q1 As for the single order quantity described in the previous section, since the rewards (costs) are stationary and bounded and the transition probabilities stationary, a stationary policy over an infinite horizon will minimize the average cost. 2.3.2 Intuitive Starting Point Intuitively, the two delivery size case should be structurally similar to the single delivery size case. Consider that we start with a single delivery size case, Q2 . We know that we are trying to reach the optimal reorder point in each period after ordering, but are often prevented from doing so by delivery size. Now we introduce the smaller delivery size Q1 = 12 Q2 . With this smaller delivery size we should have a better chance at reaching the optimal reorder point, but at added cost. The smaller delivery sizes come at a greater unit cost and so faces the incremental material cost (c1 − c2 )Q1 . With this in mind we hope to find a stationary policy, with a reorder point R, and a finite band of recurrent target states. Since the starting point is only ordering large delivery sizes, there is reasonable expectation that the target states will be a subset of the Q2 target states in the single delivery size case. The decision to order the small delivery size compares the additional material cost to the savings from reaching a lower cost inventory position. The retailer orders enough large delivery sizes to reach the band (R−Q2 , R], and then decides if a small or large delivery size should be ordered to bring the inventory position above R. In the band (R − Q2 , R − Q1 ] there is no decision to be made as only the addition of a large delivery size will cross R. This leaves the band (R − Q1 , R] in which the decision between the two delivery sizes can be made, or 2Q1 possible policies. 15 Chapter 2. Optimal Inventory Replenishment with Two Delivery Sizes Due to the quasiconvexity of the cost function, we conjectured in earlier work, Parkinson (2005), that there would either be a single point in the band (R − Q1 , R] where the optimal decision would switch from the large delivery size to the smaller, or that no small delivery sizes would be ordered. That is, moving up in inventory position from R − Q1 − 1, if it became cost effective at inventory position α to order the smaller delivery size, it would also be cost effective in positions α + 1, . . . , R. If this were true, then there would be only Q1 +1 candidates for the optimal policy once R is established. This policy was given the name the alpha-policy. In the next section we will show that there are 2Q1 candidates for optimal policy, but that the stationary distribution may change with each candidate policy making the choice of when to make the tradeoff to the small delivery size difficult. As a result, the alpha-policy argument does not hold and all 2Q1 candidate policies must be considered. In Section 2.3.5 we show a counterexample showing that the alpha-policy conjecture is false. 2.3.3 Structure of Optimal Policy The analysis of the two delivery size case follows that of the single delivery size case, but differs significantly enough to be demonstrated in detail. Recall that we are seeking a policy that minimizes P v(s) = min(a,b)∈As {r(s, (a, b)) + λ j∈S p(j|s, (a, b))v(j)}. For a given policy πP = {(as , bs )} ∈ As , let G(s + as Q2 + bs Q1 ) = G(s + as Q2 + bs Q1 ) + λ j∈S p(j|s, (a, b))v(j) and let y ∗ be an integer that minimizes G(y). We e + as Q2 + bs Q1 ) = G(s + as Q2 + bs Q1 ) + bs Q1 (c1 − c2 ) want to minimize G(s It should be clear that if s + as Q2 + bs Q1 = y ∗ and bs = 0, then e + as Q2 + bs Q1 ) is minimized at y ∗ . Intuitively this makes sense: if it G(s is possible to reach the minimum cost inventory position with the cheaper, larger delivery sizes, this should be done. Also, as with the single delivery size model, no orders should be placed from inventory positions above y ∗ . e + as Q2 + bs Q1 ) = G(s + as Q2 ) and the If bs is zero for all s, then G(s problem is the single delivery size case. The smaller delivery size allows for the tradeoff of paying the incremental cost of ordering to arrive at an inventory position with a lower one period cost. That is, the tradeoff should be considered if G(s+as Q2 +Q1 ) < G(s+as Q2 ) or G(s+(as −1)Q2 +Q1 ) < G(s + as Q2 ). In the first, the additional delivery size raises the inventory position by Q1 and the second reduces the inventory position by Q1 . This leads us to conclude that if the Q1 integer minimizers of G(y) can be reached by ordering Q2 , no savings can be realized by considering the smaller delivery size. 16 Chapter 2. Optimal Inventory Replenishment with Two Delivery Sizes Recall the definition of the Q2 integer minimizers of G(y) which would be the target states if only the large delivery size is ordered, and divide them in to two groups: the first Q1 indices {yi }, i = 1, . . . , Q1 , and the second Q1 indices {yi }, i = Q1 + 1, . . . , Q2 . We have shown that {yi }, i = 1, Q1 will be in the target inventory positions under the optimal policy. This leaves the possibility that the tradeoff will be worthwhile from the Q1 targets {yi }, i = Q1 + 1, Q2 . Since G(s) contains the future expected costs, and (as we will demonstrate) the stationary probabilities are difficult to compute, it is not simple to evaluate this tradeoff. We can summarize the results above as follows: Proposition 2.1. For the setting with two available delivery sizes Q1 and Q2 = 2Q1 with the unit costs c1 > c2 , the set of recurrent states under the optimal ordering policy will contain the set of Q1 integer minimizers of G(y). Proposition 2.2. For the setting with two available delivery sizes Q1 and Q2 = 2Q1 with the unit costs c1 > c2 , there are 2Q1 candidates for an optimal policy. These policies order into some subset of the set of Q2 integer minimizers of G(s). We are in the position of having the structure of the optimal policy, but it is not clear how to calculate it. For the single delivery size, calculation was simplified by realizing that the uniform stationary distribution implies a long run average cost that is fairly easy to compute. In the next section we will show that this property does not hold for the two delivery size case, making it difficult to evaluate when the small delivery size should be ordered. 2.3.4 Stationary Distribution of a Two Delivery Size Policy Let us now look at the long run average cost of policies of the format described above, to see if this setting provides insight into determining which tradeoffs are necessary to arrive at the optimal policy. As discussed in the previous section, for a given policy π that orders into a finite band of inventory positions, the average expected reward is g π (s) = Pπ∗ r(s). In order to track the incremental cost of ordering the smaller delivery size, we need to differentiate when a state is reached by the small delivery size, 17 Chapter 2. Optimal Inventory Replenishment with Two Delivery Sizes and so we will extend the state space to (s, b). The long run average cost becomes g π (s, b) = Pπ∗ r(s, b) where r(s, b) = G(s) + bQ1 (c1 − c2 ). If the stationary distributions of all the candidate policies happened to be uniform, the evaluation of each of the possible tradeoffs would be the comparison of G(s) to G(s ± Q1 ) + Q1 (c1 − c2 ). Unfortunately this is not the case. Only under the policy that only Q2 is ordered is the stationary distribution guaranteed to be uniform. For the other 2Q1 − 1 policies, the stationary distribution depends on the demand distribution. Starting from the Q2 only policy, there are 2Q1 − 1 policy options to consider, including ordering Q1 for any one of the yQ1 +1 , . . . , yQ2 targets, or a subset of the targets. Each of these tradeoffs affects the stationary distribution of the recurrent target states. The following illustrates how making even one of the possible tradeoffs may result in a non-uniform stationary distribution. Let us consider the candidate policy of ordering only the large delivery sizes. The set of target states will be the Q2 integer minimizers of G(y). Let ym be the largest of the target states. Table 2.2 is the transition probability matrix for this policy. 18 ∞ X ym − (Q2 − 1) P {D = jQ2 } ... ... j=0 ∞ X ym − 1 − Q1 ... P {D = jQ2 + 2} ... j=1 ym − 1 ∞ X P {D = jQ2 + 2} j=0 ym ∞ X P {D = jQ2 + 1} j=0 . . . P = ym − 1 − Q1 ∞ X P {D = jQ2 + Q1 − 2} ∞ X ... j=0 P {D = jQ2 } ... j=0 ∞ X P {D = jQ2 = Q1 } j=0 ∞ X P {D = jQ2 + Q1 − 1} j=0 . . . ym − 1 ym ∞ X j=0 ∞ X j=0 P {D = jQ2 + Q2 − 2} P {D = jQ2 + Q2 − 1} ... ... ∞ X j=0 ∞ X P {D = jQ2 + Q1 } P {D = jQ2 + Q1 + 1} j=0 ... ... ∞ X j=0 ∞ X P {D = jQ2 } P {D = jQ2 + 1} j=0 ∞ X P {D = jQ2 + Q2 − 1} j=0 ∞ X P {D = jQ2 } j=0 Table 2.2: Transition Probability Matrix for Only Ordering Large Delivery Sizes 19 Chapter 2. Optimal Inventory Replenishment with Two Delivery Sizes ym − (Q2 − 1) To State → From State ↓ Chapter 2. Optimal Inventory Replenishment with Two Delivery Sizes Let us now assume that ym − 1 is one of the yQ1 +1 , . . . , yQ2 integer minimizers of G(y) so we want to consider ordering the small delivery size to move to the target inventory position ym − 1 − Q1 which is one of the Q1 integer minimizers of G(y). In order to track the incremental cost of ordering the smaller delivery size, we need to differentiate when a state is reached by the small delivery size and we will use hQ1 i to denote a target reached by the small delivery size. Thus when ym − 1 − Q1 is reached by a small delivery size it will be referred to as state ym − 1 − Q1 hQ1 i. With this one change in ordering policy the probability transition matrix becomes that shown in Table 2.3 20 ym − (Q2 − 1) ∞ X ym − (Q2 − 1) P {D = jQ2 } ym − 1 − Q1 ... ... j=0 ∞ X P {D = jQ2 + Q2 + 2} j=0 ym − 1 − Q1 hQ1 i ∞ X P {D = jQ2 + 2} ... ... j=0 ym ∞ X P {D = jQ2 + 1} j=0 . . . ym − 1 − Q1 ym − 1 − Q1 hQ1 i ∞ X j=0 ∞ X P {D = jQ2 + Q1 − 2} P {D = jQ2 + Q1 − 2} ∞ X ... j=0 ∞ X ... j=0 P {D = jQ2 } P {D = jQ2 } j=0 ∞ X j=0 ∞ X P {D = jQ2 + Q1 } P {D = jQ2 + Q1 } ... ... j=0 ∞ X j=0 ∞ X P {D = jQ2 + Q1 − 1} P {D = jQ2 + Q1 − 1} j=0 . . . ym ∞ X j=0 P {D = jQ2 + Q2 − 1} ... ∞ X j=0 P {D = jQ2 + Q1 + 1} ∞ X j=0 P {D = jQ2 + 1} ... ∞ X P {D = jQ2 } j=0 Table 2.3: Transition Probability Matrix for Ordering a Small Delivery Size at Only One Inventory Position 21 Chapter 2. Optimal Inventory Replenishment with Two Delivery Sizes To State → From State ↓ Chapter 2. Optimal Inventory Replenishment with Two Delivery Sizes To summarize the changes in the transition matrix resulting from the change in the policy • state ym − 1 − Q1 hQ1 i has been added. The entries of this row will be exactly the same as those of the row for state ym − 1 − Q1 since the probability transitions are independent of how the originating state was reached. • the values in the column for state ym − 1 have been moved to the column ym − 1 − Q1 hQ1 i since the latter is now reached instead of the former • the values in the column for state ym − 1 have been set to zero as there is now no way to reach this state under the revised policy. In the new transition matrix, the probabilities in each row still sum to one – with probability one leaving any single state will end up at one of the other states, however, the probabilities in each column might not. More importantly, the probabilities in each column do not sum to the same value. The solution to q T Pπ = q now depends on the demand distribution and we do not see any general closed-form for stationary distributions PQ2expression −1 here. The cost for this policy is i=0,i6=1 q(ym − i)G(ym − i) + q(ym − 1 − Q1 hQ1 i)(G(ym − 1 − Q1 ) + Q1 (c1 − c2 )) and the tradeoff between holding and material costs is not transparent. Clearly, the stationary distribution will only be uniform under very specific demand distributions. Let us now consider the policy when the small delivery size is ordered whenever necessary to reach the Q1 integer minimizers of G(·). In this case, while each column of the probability transition matrix does not sum to one, the sum of each pair of columns yi , yi hQ1 i does sum to one and the result is that for each i = 1, . . . , Q1 q(xi ) + q(xi hQ1 i) = Q11 . The long run P 1 PQ1 average cost for this policy is Q11 Q i=1 G(xi ) + i=1 q(xi hQ1 i)Q1 (c1 − c2 ). This will be used to create upper and lower bounds on the optimal cost in Section 2.3.7. Clearly the stationary distribution is not guaranteed to be uniform and computation (solving q = q T P ) is necessary to determine the long run average cost of each of the 2Q1 candidate policies. 2.3.5 Counterexamples to Optimality of Alpha-Policy The previous section suggests that the optimal policy will not always be an alpha-policy. Here we give two examples demonstrating that alpha-policies 22 Chapter 2. Optimal Inventory Replenishment with Two Delivery Sizes are not always optimal. In the first we show that the selection of R can change the appeal of the alpha-policies, and in the second we show a case where even the selection of R does not make the alpha-policy optimal. For the first counterexample, consider the following instance: • Q1 = 3 • h = 25, b = 50 • c1 − c2 = 5 • demand is uniformly distributed over the integers [47, 56], U (47, 56), 0.1 for x = 47, 48, . . . , 56 i.e., P {D = x} = 0 otherwise To compute the optimal policy we must identify the Q2 minimizers of the cost function then evaluate the 2Q1 candidate policies. In this example we use the following analysis to identify the minimizers. The one period expected holding cost given that the inventory position after ordering is y is a combination of the expected holding and backorder costs to be paid at the end of the period after demand D is realized. Under linear costs, a unit holding cost of h is applied to each item held, and unit backorder cost of b is applied to each item in shortfall. Using (y)+ = max{y, 0} and E[y] as the expected value of y, the cost function is given by G(y) = E[h(y − D)+ + b(D − y)+ ] G(y) = (h + b)(yP {D ≤ y} − y X dP {D = d}) − b(y − E[D]) d=0 In this analysis we will use the following: G(y) − G(y − 1) = (h + b)(yP {D ≤ y} − y X dP {D = d}) + b(E[D] − y) d=0 − [(h + b)((y − 1)P {D ≤ (y − 1)} − y−1 X dP {D = d}) + b(E[D] − (y − 1))] d=0 = (h + b)(yP {D ≤ y} − (y − 1)P {D ≤ (y − 1)} − yP {D = y}) − b = (h + b)(yP {D = y} + P {D ≤ y − 1} − yP {D = y}) − b = (h + b)P {D ≤ y − 1} − b 23 Chapter 2. Optimal Inventory Replenishment with Two Delivery Sizes If y ∗ is a minimizer then G(y ∗ ) < G(y ∗ + 1) and G(y ∗ ) < G(y ∗ − 1). Starting with the first, G(y ∗ ) − G(y ∗ + 1) < 0 (h + b)P {D ≤ y ∗ } − b > 0 P {D ≤ y ∗ } > b h+b And the second G(y ∗ ) − G(y ∗ − 1) < 0 (h + b)P {D ≤ y ∗ − 1} − b < 0 P {D ≤ y ∗ − 1}) < b h+b So if G(y) minimized at y ∗ then P {D ≤ y ∗ − 1} < b < P {D ≤ y ∗ } h+b b = 0.66, and P {D ≤ 53} = .7 and In our example, P {D ≤ 52} = 0.6, h+b we identify 53 as the minimizer. Evaluating G(y) in the vicinity of y = 53 G(50) = 120 G(51) = 100 G(52) = 87.5 G(53 = y ∗ ) = 82.5 and establish our Q2 minimizers are [51, 56], we see G(54) = 85 G(55) = 95 G(56) = 112.5 G(57) = 137.5 G(53) < G(54) < G(52) < G(55) < G(51) < G(56), and R = 50. From the demand distribution, the entries in the probability transition matrices are ∞ ∞ X X P {D = iQ2 } = 0.2 P {D = 1 + iQ2 } = 0.2 i=0 ∞ X i=0 P {D = 2 + iQ2 } = 0.2 i=0 ∞ X i=0 ∞ X P {D = 3 + iQ2 } = 0.1 i=0 P {D = 4 + iQ2 } = 0.1 ∞ X P {D = 5 + iQ2 } = 0.2 i=0 We have 23 = 8 candidates for the optimal policy as listed in Table 2.4 24 Chapter 2. Optimal Inventory Replenishment with Two Delivery Sizes Policy 1 2 3 4 5 6 7 8 51 1 1 1 1 1 1 1 1 52 1 1 1 1 1 1 1 1 Target Inventory Positions 53 54 55 56 53 hQ1 i 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 54hQ1 i 52hQ1 i 1 1 1 1 1 1 1 1 Table 2.4: Candidates for Optimal Policy: Counterexample 1 Across the top are first listed the Q2 inventory positions that minimize the cost function, followed by the Q1 inventory positions that can be reached by ordering the smaller delivery size. In this example inventory positions 51 through 56 are the 6 minimizers, then 53, 54, 52 can also be reached by ordering small delivery sizes. In the first policy only large delivery sizes are ordered and all inventory positions 51 through 56 are reached. In the second policy, the choice has been made to order a small delivery size to reach inventory position 53 instead of reaching 54 with the large delivery size. The remaining policies enumerate the combination of choices for small delivery size. Note that the alpha-policies allow a switch from the large delivery size to the small delivery size at or below inventory position R, in the band (R − Q1 , R]. In this case, only two of the Q1 = 3 minimizers can be reached by ordering a small delivery size from at or below inventory position R = 50. The alpha-policies are the subset of the candidates shown in Table 2.5 Policy 1 2 6 51 1 1 1 52 1 1 1 Target Inventory Positions 53 54 55 56 53 hQ1 i 1 1 1 1 1 1 1 1 1 1 1 54hQ1 i 52hQ1 i 1 Table 2.5: Alpha Policy Candidates: Counterexample 1 25 Chapter 2. Optimal Inventory Replenishment with Two Delivery Sizes The stationary distributions for each policy, and the average cost are given in the Table 2.6 . Policy 1 2 3 4 5 6 7 8 51 0.1667 0.1687 0 0.1852 0 0.1852 0 0 52 0.1667 0.1852 0.1687 0.1831 0.1852 0.2 0.1852 0.2 53 0.1667 0.1831 0.1852 0.1646 0.2 0.1815 0.1852 0.2 Target Inventory Positions 54 55 56 53 hQ1 i 0.1667 0.1667 0.1667 0 0.1646 0.1481 0 0.1502 0.1831 0.1646 0.1481 0 0.1481 0 0.1687 0 0.1815 0.1481 0 0.1333 0.1481 0 0 0.1519 0.1667 0 0.1481 0 0.1667 0 0 0.1333 54 hQ1 i 0 0 0.1502 0 0.1519 0 0.1667 0.1667 52 hQ1 i Avg Cost 0 93.77 0 90.89 0 92.92 0.1502 95.08 0 90.39 0.1333 92.06 0.1481 94.17 0.1333 91.50 Table 2.6: Stationary Distribution and Average Cost of Candidate Policies: Counterexample 1 Policy 5 has the lowest cost at 90.93, and we conclude that it is optimal to order the large delivery size at every inventory position at or below 49, and to order the small delivery size at inventory positions 50 and 51. The example is illustrated in Figure 2.2 below which shows G(y) for this counterexample 1. The Q1 = 3 minimizers are shown as solid circles {52, 53, 54} and the remaining 3 minimizers are shown as open circles {51, 55, 56}. R is indicated with the vertical line at inventory position 50. Along the bottom are indicators of the delivery size ordered from an inventory position, the open triangles and squares, and the delivery size ordered to arrive at a target state, closed triangles and squares. The open triangles show the Q1 = 3 inventory positions from which Q1 can be ordered to reach one of the Q1 minimizers, and the open squares show the Q1 = 3 inventory positions from which Q2 can be ordered to reach one of the Q1 minimizers. In this run, Q2 was ordered from inventory positions at or below 49 to reach target states {52, . . . , 54, 55}, and Q1 is ordered from inventory positions {50, 51} to reach target state {53, 54}. We see that at inventory position 49 the large batch is ordered and an inventory position outside of the set of Q1 minimizers is reached. 26 Chapter 2. Optimal Inventory Replenishment with Two Delivery Sizes Figure 2.2: Graphical Representation of Optimal Policy. The Q1 = 3 minimizers are shown as solid circles {52, 53, 54} and the remaining 3 minimizers are shown as open circles {51, 55, 56}. The open triangles and squares show the inventory positions where the last order is placed. The triangles indicate where Q1 can be ordered to reach a Q1 minimizer, and the squares indicate where Q2 can be ordered to reach a Q1 minimizer. The solid triangles and squares show the inventory positions in the set of recurrent states for this policy. 27 Chapter 2. Optimal Inventory Replenishment with Two Delivery Sizes While one may argue that instead of being a counterexample, the example above instead makes a case for somehow tweaking the definition of the alpha-policy, the next example will show that the optimal policy can have gaps that no tweaking can cover. For the second counterexample, consider the following instance: • Q1 = 3 • h = 25, b = 50 • c1 − c2 = 5 • Demand distribution, P {D = x} = 0 except for the cases listed below. P {D = 23} = 0.09667 P {D = 27} = 0.08467 P {D = 28} = 0.08467 P {D = 29} = 0.064 P {D = 48} = 61 P {D = 49} = 61 P {D = 53} = 0.10267 P {D = 57} = 0.082 P {D = 58} = 0.082 Following the steps from the previous example, we first identify the minimizer of G(y), y ∗ = 49. 50 2 = < P {D ≤ 49} 75 3 Next we calculate the single period expected holding and back order costs of given inventory position after ordering around the minimizer. G(47) = 348 G(48) = 328 G(49) = 320.5 G(50) = 325.5 G(51) = 330.5 G(52) = 335.5 G(53) = 340.5 G(54) = 353.2 From the above list, the cost function is minimized at 49. We can order the minimizers as follows: G(49) < G(50) < G(48) < G(51) < G(52) < G(53) From the demand distribution, the entries in the probability transition matrices are ∞ ∞ X X P {D = iQ2 } = 0.1667 P {D = 1 + iQ2 } = 0.1667 P {D ≤ 48} < i=0 ∞ X i=0 ∞ X i=0 i=0 P {D = 2 + iQ2 } = 0.07 P {D = 4 + iQ2 } = 0.1667 ∞ X i=0 ∞ X P {D = 3 + iQ2 } = 0.1667 P {D = 5 + iQ2 } = 0.2633 i=0 28 Chapter 2. Optimal Inventory Replenishment with Two Delivery Sizes We again have 23 = 8 candidates for the optimal policy as shown in Table 2.7 Policy 1 2 3 4 5 6 7 8 48 1 1 1 1 1 1 1 1 49 1 1 1 1 1 1 1 1 Target Inventory Positions 50 51 52 53 48 hQ1 i 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 49hQ1 i 50hQ1 i 1 1 1 1 1 1 1 1 Table 2.7: Candidate Policies: Counterexample 2 The stationary distributions for each policy, and the average cost are given in Table 2.8. Policy 1 2 3 4 5 6 7 8 48 0.1667 0.1655 0.1604 0.1344 0.1604 0.1344 0.1344 0.1344 49 0.1667 0.1989 0.1655 0.1604 0.1989 0.1989 0.1604 0.1989 50 0.1667 0.1729 0.1989 0.1655 0.1989 0.1729 0.1989 0.1989 Target Inventory Positions 51 52 53 48 hQ1 i 0.1667 0.1667 0.1667 0.1344 0.1604 0.1679 0.1729 0.1344 0.1989 0.1729 0.1344 0.1729 0.1344 0.1989 0.1989 0.1989 49hQ1 i 50hQ1 i Avg Cost 330.083 331.608 0.1679 329.634 0.1679 330.275 0.1344 331.293 0.1604 332.073 0.1729 0.1344 329.774 0.1344 0.1344 331.683 Table 2.8: Stationary Distribution and Average Costs: Counterexample 2 While the difference is small, policy 3 is the lowest cost policy. It qualifies as a non-alpha-policy since one should order the small delivery size at inventory position 52 − 6 = 46, but at inventory position 47 the large delivery size should be ordered. This creates a gap in the set of recurrent states at inventory position 46. This example shows the intuition behind the alpha-policy that once a transition to the small delivery size is made it is not reversed, does always not hold, thus the alpha-policy conjecture is false. 29 Chapter 2. Optimal Inventory Replenishment with Two Delivery Sizes The example is illustrated in Figure 2.3 below which shows G(y) for counterexample 2. The Q1 = 3 minimizers are shown as solid circles {48, 49, 50} and the remaining 3 minimizers are shown as open circles {51, 52, 53}. R is indicated with the vertical line at inventory position 47. Along the bottom are indicators of the delivery size ordered from an inventory position, the open triangles and squares, and the delivery size ordered to arrive at a target state, closed triangles and squares. In this run, Q2 was ordered from inventory positions {42, . . . , 45, 47} to reach target states {48, . . . , 51, 53}, and Q1 is ordered from inventory positions 46 to reach target state 49. Figure 2.3: Graphical representation of optimal policy. The Q1 = 3 minimizers are shown as solid circles and the remaining 3 minimizers are shown as open circles. The open triangles and squares show the inventory positions where the last order is placed. The triangles indicate where Q1 can be ordered to reach a Q1 minimizer, and the squares indicate where Q2 can be ordered to reach a Q1 minimizer. The solid triangles and squares show the inventory positions in the set of recurrent states for this policy. Note that at inventory position 46 a small delivery size is ordered, but at inventory position 47 the large delivery size is ordered. 30 Chapter 2. Optimal Inventory Replenishment with Two Delivery Sizes 2.3.6 An Algorithm for Computing Exact Solutions As demonstrated, the stationary distribution of a policy depends on the demand distribution. We have shown that there are 2Q1 possible candidates for optimal policy, and each candidate can be evaluated by finding the stationary distribution and evaluating the long run average cost. Let us denote the 2Q1 possible candidates for optimal policy as D∗ . For small values of Q1 , evaluating each candidate policy may be feasible, but as Q1 grows, this quickly becomes infeasible. For Q1 = 10 there are 1024 candidate policies. Ideally one would look at each of the Q1 possible inventory positions where the small delivery size can be ordered, and decide if the tradeoff is warranted. Instead, we turn to conventional methods for solving Markov decision processes. Puterman (1994) describes some algorithms used for solving Markov decision processes. The basis for these algorithms is the Bellman or optimality equations. Since we have a finite set of candidate policies, we choose the policy iteration of Howard (1960). Each step of policy iteration uses a modified form of the Bellman equations to search through candidate policies for one that has an improvement in value. We also show that the structure of policy iteration provides a separation of the problem that allows for evaluation of each of the Q1 tradeoffs individually. In a Markov process states are classified by the expected traffic to the states. The traffic to a state s is described with two random variables, the number of visits, vs and the time of the first visit τs . If the chain starts at state s then τs represents the time of the first return visit. A state is recurrent if it is expected that it will be visited repeatedly. That is, E[τs ] < ∞. Otherwise, if E[τs ] = ∞ the state is transient and the state will rarely be visited. Using the number of visits, a state is recurrent if and only if E[vs ] = ∞, that is, there is nothing that stops visits to the state. Recurrent states are then grouped into sets based on accessibility between the states. Two states are accessible if there is a non-zero probability of transitioning between them in any number of periods. Formally, state i is accessible from state j if pn (i|j) > 0 for some n ≥ 0. For a finite state space where all states are recurrent and accessible, the Markov chain is referred to as unichain. Since the probability matrix for our problem consists of a single recurrent class for each of the candidate policies, the model can be classified as unichain. The Unichain Policy Iteration Algorithm provided by Puterman converges to an optimal solution in a finite number of iterations for a unichain model. The algorithm is started with an arbitrary decision rule/policy d, for which the stationary distribution Pd∗ and average reward 31 Chapter 2. Optimal Inventory Replenishment with Two Delivery Sizes g = Pd∗ rd are computed. Next 0 = rd − ge + (Pdn − I)hn and Pd∗n hn = 0 are solved to obtain the vector h, referred to as the bias. The bias is the expected difference between the reward in each state and the average reward under a given policy. In the policy improvement step dn+1 ∈ arg max {rd + Pd hn } d inD is computed. In this step, the bias from the previous iteration is used to estimate the average reward of the other candidate policies, and the best is selected. The algorithm terminates when there is no improvement possible. Since the algorithm only proceeds if there is an improvement in average value, no candidate policy is evaluated more than once and the algorithm is guaranteed to converge to the optimal solution if there is a finite number of policies. To start the Unichain Policy Iteration Algorithm a decision rule must be chosen and the set of target states determined. We have already shown that the optimal policy will only target some subset of the Q2 integer minimizers of G(·). Since Q1 of these may be reached both by ordering Q2 or by ordering the additional Q1 , there are Q2 + Q1 target states if we also include in the state definition the number of small delivery sizes that must be ordered to reach it. These 3Q1 states are the • the Q2 integer minimizers of G(·), {x1 , . . . , xQ2 }, and • the Q1 states reached by ordering a small delivery size {x1 hQ1 i, . . . , xQ1 hQ1 i} By choosing the policy in which no small delivery sizes are ordered, the P 2 calculations in the first iterations can be simplified. First, g = Q G(x i) i=1 is easily calculated, and second, the stationary distribution P ∗ is the square matrix with all entries in the first Q2 columns Q12 and zero in all remaining entries. The reward rd i is the cost for spending a period in state i, under policy d. In this model, the cost for spending the period at one of the Q2 integer minimizers after ordering a large delivery size is policy independent and rd xi = rxi = G(xi ). Similarly, for an inventory position reached by ordering a small delivery size, rd xi hQ1 i = rxi hQ1 i = G(xi ) + Q1 (c1 − c2 ). The calculation of h is further simplified by the observation that since rxj hQ1 i = rxj + Q1 (c1 − c2 ), it follows that hn xj hQ1 i = hn xj + Q1 (c1 − c2 ) for j = 1, . . . , Q1 . 32 Chapter 2. Optimal Inventory Replenishment with Two Delivery Sizes Once h is calculated, the policy improvement step is next, and here is where each tradeoff can be individually evaluated, thus reducing the need to evaluate all 2Q1 policies. Since rdi does not change with policy, policy improvement is reduced to inspecting dn+1 ∈ arg maxd∈D {Pd hn }. For each state i, the improvement of policy d over dn is Q2 X (pd i,xj − pdn i,xj )hn xj + j=1 Q1 X (pd i,xj hQ1 i − pdn i,xj hQ1 i )hn xj hQ1 i . j=1 Nominally, this is only 3Q1 comparisons, instead of the 2Q1 −1 other policies in D∗ . Fortunately a few simplifications can be made to reduce the number of comparisons down to Q1 as a result of the structure of the problem. Recall the description of the construction of the transition matrices for policies where small delivery sizes are ordered, described above, using the transition matrix for the ordering only large delivery size policy as a starting point. For replacement of xj with state xi hQ1 i, a row is added with the same probabilities as xi . Next, a new column taking the entries of the xj column is added and the entries in column xj set to zero. This construction highlights the following: • for any policy d ∈ D∗ , pd i,xj = pdn i,xj for j = 1, . . . , Q1 . The transition probabilities from any state to one of the Q1 minimizers does not change with policy since these columns are unchanged. Now, for each state i, the improvement of policy d over dn is Q2 X (pd i,xj − pdn i,xj )hn xj + j=Q1 +1 Q1 X (pd i,xj hQ1 i − pdn i,xj hQ1 i )hn xj hQ1 i . j=1 • for each pair of states xj hQ1 i = xk ± Q1 , j = 1, . . . , Q1 , k = Q1 + 1, . . . , Q2 , only one can be in the final set of recurrent states. When comparing policy d to dn either xk in d but the switch to xj hQ1 i has been made in policy dn . In this case pdn i,xj hQ1 i = pd i,xk and pd i,xj hQ1 i = pdn i,xk = 0 The contribution of this pair of states to potential improvement of policy dn is (pd i,xk − pdn i,xk )hn xk + (pd i,xj hQ1 i − pdn i,xj = pd i,xk hn xk − pdn i,xj hQ1 i )hn xj hQ1 i hQ1 i hn xj hQ1 i = pd i,xk (hn xk − hn xj hQ1 i ) 33 Chapter 2. Optimal Inventory Replenishment with Two Delivery Sizes or xj hQ1 i in d but the switch to back to xk has been made in policy dn . In this case pdn i,xk = pd i,xj hQ1 i and pd i,xk = pdn i,xj hQ1 i = 0. The contribution of this state pair to policy improvement becomes pd i,xj hQ1 i (−hn xk + hn xj hQ1 i ) or xj in both pdn i,xj hQ1 i = pdn i,xk = 0 and pdn i,xk = pd0 i,xk . With no change in policy for these states there is potential improvement. Combining these, policy improvement at the end of the first iteration, can be reduced to looking at the Q1 differences hn xj <Q1 > − hn xk j ∈ [Q1 + 1, Q2 ]. If hn xj <Q1 > − hn xk < 0 then the tradeoff should be made to take the small delivery size when at state xk as it will result in a reduction in the long run cost of the policy. In further iterations, evaluation of the potential contribution for each pair of points as above, determines if there is a policy with lower costs. This simplification is illustrated using the following example in which Q1 = 2 and a cost function G(·) such that G(20) < G(21) < G(22) < G(19) < G(y) y ∈ Z, y 6= {21, 22, 23, 24}. In this example there will be 4 candidate policies with the following probability transition matrices. Policy: Order all Q2 19 20 21 22 20 hQ1 i 21hQ1 i 19 0.1 0.4 0.3 0.2 0.4 0.3 20 0.2 0.1 0.4 0.3 0.1 0.4 21 0.3 0.2 0.1 0.4 0.2 0.1 22 20hQ1 i 21hQ1 i 0.4 0 0 19 0.3 0 0 20 0.2 0 0 21 0.1 0 0 22 0.3 0 0 20 hQ1 i 0.2 0 0 21hQ1 i Policy: Order Q1 at 21 19 20 21 22 20 hQ1 i 21hQ1 i 19 0.1 0.4 0.3 0.2 0.4 0.3 20 0.2 0.1 0.4 0.3 0.1 0.4 21 0.3 0.2 0.1 0.4 0.2 0.1 Policy: Order Q1 at 19 19 0 0 0 0 0 0 20 0.2 0.1 0.4 0.3 0.1 0.4 21 0.3 0.2 0.1 0.4 0.2 0.1 22 20hQ1 i 21hQ1 i 0.4 0 0.1 0.3 0 0.4 0.2 0 0.3 0.1 0 0.2 0.3 0 0.4 0.2 0 0.3 Policy: Order all Q1 22 20hQ1 i 21hQ1 i 0 0.4 0 19 0 0.3 0 20 0 0.2 0 21 0 0.1 0 22 0 0.3 0 20 hQ1 i 0 0.2 0 21hQ1 i 19 0 0 0 0 0 0 20 0.2 0.1 0.4 0.3 0.1 0.4 21 0.3 0.2 0.1 0.4 0.2 0.1 22 20hQ1 i 21hQ1 i 0 0.4 0.1 0 0.3 0.4 0 0.2 0.3 0 0.1 0.2 0 0.3 0.4 0 0.2 0.3 If we start the algorithm with the policy of ordering all Q2 as d0 we can compare the improvement of the other policies as (Pd − Pd0 )h0 . Let us use h0s to indicate the state s component of the bias vector from the 34 Chapter 2. Optimal Inventory Replenishment with Two Delivery Sizes initialization step. The improvement for ordering the small delivery size at inventory position 19 is −.1 −.4 −.3 −.2 (h019 − h021hQ1 i ), −.4 −.3 for ordering the small delivery size at inventory position 22 −.4 −.3 −.2 −.1 (h022 − h020hQ1 i ), −.3 −.2 and for always ordering small delivery sizes −.4 −.1 −.3 −.4 −.2 −.3 −.2 (h019 − h021hQ1 i ) + −.1 −.3 −.4 −.3 (h0 − h0 ). 22 20hQ1 i −.2 Clearly, if h019 − h021hQ1 i < 0 then the tradeoff to the smaller delivery size at inventory position 19 should be made for the next iteration. Similarly, if h022 − h020hQ1 i < 0 then the tradeoff to ordering the smaller delivery size to reach inventory position 22 should be made for the next iteration. By evaluating the difference in bias between pairs of states for each opportunity to order the smaller delivery size, the policy that satisfies dn+1 ∈ arg maxd inD {rd + Pd hn } is quickly found. The improvement for each tradeoff consideration is additive and the policy improvement step is reduced to the consideration of Q1 tradeoffs. While the evaluation of each step is fast, the consideration of Q1 tradeoffs, it may be necessary for the algorithm to evaluate all of the 2Q1 possible candidates for optimal policy before completion. In practice, however, this algorithm converges quickly. In the eighteen examples worked to date with 35 Chapter 2. Optimal Inventory Replenishment with Two Delivery Sizes linear cost functions, 4 converge in one step and the others converge in at most 4 steps. These examples are discussed in Section 2.4 In this section we have shown that one pass of the Unichain Policy Iteration Algorithm can easily provide a method of evaluating the tradeoff of ordering a small delivery size at each of the inventory positions where it may be considered. Unfortunately we can’t prove that it will converge in one pass of the algorithm. If convergence is not reached in a small number of iterations, then it may be sufficient to bound the results. The following section provides easily calculated upper and lower bounds on the value of the optimal policy. 2.3.7 Upper and Lower Bounds Faced with the potential evaluation of many iterations of the Unichain Policy Iteration Algorithm or finding the stationary distribution of 2Q1 policies to find the optimal policy, we have developed upper and lower bounds on the cost of optimal policies. This section describes the bounds and demonstrates that the gap between the bounds is small for reasonable examples. Proposition 2.3. For the setting with twoPavailable delivery sizes Q1 and 1 Q2 = 2Q1 with the unit costs c1 > c2 , Q11 Q i=1 G(xi ) is a lower bound on the long run average cost of the optimal policy. Proof. When only Q1 is available to be ordered, the optimal policy is to order to the set of PQ11 integer minimizers of G(·) resulting in a long run average cost of Q11 Q i=1 G(xi )+ ordering. Adding the possibility of ordering Q2 allows for the tradeoff of higher holding cost, P 1 for a savings in ordering cost. While there is not a policy with cost Q11 Q i=1 G(xi ), it is a lower bound on the policy. Proposition 2.4. For the setting with two availablePdelivery sizes Q1 and 1 Q2 = 2Q1 with the unit costs c1 > c2 , min{ Q11 Q i=1 G(xi ) + Q1 (c1 − P 2 c2 ), Q12 Q i=1 G(xi )} is an upper bound on the long run average cost of the optimal policy. Proof. The long run averagePcost of the optimal policy when only the large 2 delivery size is ordered is Q12 Q i=1 G(xi ), thus this is a feasible solution and is an upper bound. When Q1 is ordered so that the set of Q1 integer minimizers 1 PQ1 of G(·) are always reached, the long run average cost is Q1 i=1 G(xi ) + PQ1 1 i=1 q(xi hQ1 i)Q1 (c1 − c2 ) for 0 ≤ q(xi hQ1 i) ≤ Q1 as shown at the P P Q1 1 end of Section 2.3.4. Since Q11 Q i=1 G(xi ) + i=1 q(xi hQ1 i)Q1 (c1 − c2 ) ≤ 36 Chapter 2. Optimal Inventory Replenishment with Two Delivery Sizes 1 Q1 PQ1 + Q1 (c1 − c2 ) this is an upper bound to this policy. Thus, the minimum of the two upper bounds is an upper bound on the optimal policy. i=1 G(xi ) The gap between the upper and lower bounds is min{ Q1 Q1 Q2 Q1 1 X 1 X 1 X 1 X G(xi )+Q1 (c1 −c2 )− G(xi ), G(xi )− G(xi )} Q1 Q1 Q2 Q1 i=1 i=1 i=1 i=1 Q2 Q1 X X 1 = min{Q1 (c1 − c2 ), ( G(xi ) − G(xi ))}. Q2 i=Q1 +1 i=1 The gap is the minimum of the incremental material cost of the small delivery size, and the average incremental one period costs under the large delivery size. This is the fundamental tradeoff and the extremes are the basis P for the bounds. Unfortunately, this gap could be arbitrarily bad 1 if Q11 Q i=1 G(xi ) is close to zero. However, it seems unlikely that P Q2 1 1 PQ1 i=1 G(xi ) is close to i=Q1 +1 G(xi ) would be excessively large if Q1 Q2 zero. The following numerical study provides some evidence that the gap is can be small in practice, although some instances had a gap 32% of the optimal cost. 2.4 Numerical Study An numerical study was conducted to demonstrate the magnitude of the gap between the upper and lower bounds and compare the solution techniques. A linear holding and backorder model was considered where the one period expected cost after ordering is given by G(y) = hE[y −d]+ +bE[d−y]+ . The demand was assumed to follow either a normal or uniform distribution. The parameters for the runs were chosen to create instances that varied not only with delivery size, but also with relative costs and demand. Small delivery size ranged from 3 to a maximum of 20 and the incremental cost of ordering the small delivery size, c2 − c1 was either 1,3 or 5. Two levels of holding and back order costs were used; low, approximately 5, and high, 20 to 30. The spread of the demand distribution relative to the delivery size was also varied. Ten instances were selected and run with uniformly distributed demand for Runs 1 through 10, then with normally distributed demand for Runs 11 through 20. 37 Chapter 2. Optimal Inventory Replenishment with Two Delivery Sizes In each run, we computed the upper and lower bounds and solved the instance using the following three algorithms 1. Computation of stationary distributions and total cost of each of the 2Q1 candidates for optimal policy 2. Unichain Policy Iteration Algorithm (UPIA) as described in Section 2.3.6 3. Computation of stationary distributions and total cost of the up to Q1 + 1 alpha-policy candidates Table 2.9 shows the results for runs with uniform demand distribution and Table 2.10 the results for the runs with normal demand distribution. Matlab was used for the computations. We use the following notation in the presentation of the results. ∗ g is the average cost of the optimal policy (obtained through enumeration) P 1 UB Q1 Q11 Q i=1 G(xi ) + Q1 (c1 − c2 ) 1 PQ2 UB Q2 Q2 i=1 G(xi ) UB min{ UB Q1 , UB Q2 } 1 PQ1 LB i=1 G(xi ) Q1 P ∗ g the average cost of the policy found through the UPIA g P1 is the average cost of the policy from one iteration of the UPIA Note: The 0th iteration will have value UB Q2 # It is the number of iterations required for the UPIA to converge Note: Starting with ordering all Q2 is considered the 0th iteration. One additional iteration is required to confirm convergence. gα is the average cost of the policy found from evaluating only at most Q1 + 1 alpha-policies In the cases run, the gap between the upper and lower bound ranges between 0.12% and 32% of the optimal solution and did not confirm the proposition that the gap between bounds is small. The success of policy iteration is also demonstrated with very few iterations required for the algorithm to converge. Of the instances run, some converged after one iteration and the maximum number of iterations required was 4. In the instance with Q1 of 20, UPIA converged in 3 iterations and we were unable to complete the 220 = 1 048 476 evaluations necessary for full enumeration. 38 Chapter 2. Optimal Inventory Replenishment with Two Delivery Sizes Run 1 2 3 4 5 6 7 8 9 Q1 c1 c2 3 20 19 5 20 19 7 15 12 3 20 15 15 20 19 5 20 19 10 20 19 15 20 19 20 20 19 h b 30 30 30 30 22 33 20 20 30 30 5 3 5 3 5 3 5 3 Demand Uniform (65,135) (90,110) (77,123) (85,115) (90,110) (40,60) (40,60) (40,60) (40,60) g∗ UB Q1 UB Q2 UB LB 533.60 535.68 533.73 533.73 532.68 162.38 165.00 169.29 165.00 160.00 319.22 333.45 319.59 319.59 312.45 156.88 170.27 156.88 156.88 155.27 192.38 198.81 261.67 198.81 183.81 21.21 25.05 21.21 21.21 20.05 UB-LB U B−LB g∗ U B−LB LB 1.06 0.2% 0.2% 5.00 3.08% 3.13% 7.14 2.24% 2.28% 1.61 1.03% 1.04% 15.00 7.8% 8.16% 1.17 4.74 9.77 15.21 5.5% 19.1% 32.2% 5.82% 22.4% 42.1% 58.6% gP ∗ g P1 P1 g − g∗ # It 533.60 533.60 0 2 162.38 162.38 0 2 319.22 319.22 0 2 156.88 156.88 0 1 192.38 193.29 0.914 4 21.21 21.21 0 1 24.81 24.81 0.000 2 30.37 30.47 0.093 3 gα g − g∗ 533.60 0 167.38 5.0 319.50 0.275 156.88 0 242.12 49.74 21.21 0 25.31 0.505 31.62 1.24 Solve Time (s) UPIA 0 all 2Q1 0.016 α 0 0 0.015 0.016 0.016 0.141 0.016 0 0 0 0.031 531.20 0.047 0 0.016 0 0.015 0.031 0.063 2.032 529.05 0.015 0.062 α 24.81 31.21 25.96 25.96 21.21 30.37 38.22 32.99 32.99 23.22 45.96 41.17 41.17 25.96 35.39 35.70 3 Table 2.9: Runs with Uniform Distribution The UPIA proved most effective for the large values of Q1 where complete enumeration was time consuming or not possible (Runs 9 and 18). These runs also demonstrate that the minimum value alpha-policy is not always optimal, with the largest gap of 25% from optimality. Note that these runs provide further counterexamples to the alpha-policy conjecture beyond those in Section 2.3.5 Table 2.11 illustrates the first step of the UPIA for three instances that converged in one iteration. Recall that the 0th iteration assumes only the 39 Chapter 2. Optimal Inventory Replenishment with Two Delivery Sizes Run 10 11 12 13 14 15 16 17 18 Q1 c1 c2 3 20 19 5 20 19 7 15 12 3 20 15 15 20 19 5 20 19 10 20 19 15 20 19 20 20 19 h b 30 30 30 30 22 33 20 20 30 30 5 3 5 3 5 3 5 3 Demand Normal (100,35) (100,10) (100,23) (100,15) (100,10) (50,10) (50,10) (50,10) (50,10) g∗ UB Q1 UB Q2 UB LB 739.95 742.09 739.96 739.96 739.09 243.95 246.95 249.02 246.95 241.95 477.50 492.78 477.57 477.57 471.78 240.51 254.44 240.51 240.51 239.44 268.69 276.29 320.29 276.29 261.29 31.57 35.64 31.57 31.57 30.64 UB-LB U B−LB g∗ U B−LB LB 0.87 0.12% 0.12% 5.00 2.05% 2.07% 5.79 1.21% 1.23% 1.07 0.44% 0.45% 15.00 5.58% 5.74% 0.94 3.56 7.47 12.12 2.96% 10.2% 19.2% 3.05% 11.3% 22.6% 34.5% gP ∗ g P1 P1 g − g∗ # It 739.95 739.95 0 2 243.95 243.95 0 2 477.50 477.50 0 2 240.51 240.51 0 1 268.69 268.69 0 2 31.57 31.57 0 1 34.83 34.83 0 2 38.95 38.96 0.013 3 gα gα − g∗ 739.95 0 246.49 2.538 477.57 0.0747 240.51 0 296.88 28.19 31.57 0 35.06 0.222 39.67 0.71 Solve Time (s) UPIA 0.015 0 all 2Q1 α 0 0 0.032 0 0.016 0.141 0.015 0 0.016 0 0.032 546.3 0.109 0 0.016 0 0.015 2.313 0.015 0.031 530.4 0.078 34.83 41.57 35.13 35.13 31.57 38.95 48.07 40.54 40.54 33.07 55.13 47.24 47.24 35.13 43.56 43.63 3 0.062 Table 2.10: Runs with Normal Distribution large delivery size is ordered and has uniform stationary distribution of 1 Q2 for each state. The table shows the h values for each of the integer minimizers where the decision to order the small delivery size can be made, and compares it to the h value at its tradeoff pair. Note yj denotes the states that are candidates for tradeoffs, i.e., {yj } = {xi } i ∈ [Q1 + 1, Q2 ] Where the difference in h value, i.e., ∆ = hxi hQ1 i − hy1 is negative, a switch to ordering the small delivery size is made. The lower section of the table enumerates the 2Q1 candidate policies and value, identifying them by the 40 Chapter 2. Optimal Inventory Replenishment with Two Delivery Sizes which tradeoffs, if any, to the small delivery size have been made. In the three instances shown, the first iteration policies match the optimal policies found through enumeration. Run 1 Normal (100,35) Q1 = 3, c1 = 20, c2 = 19 h = 30, b = 30 4 Normal (100,15) Q1 = 3, c1 = 20, c2 = 15 h = 20, b = 20 Uniform (91, 105) Q1 = 3, c1 = 20, c2 = 19 h = 30, b = 30 hy1 hxi hQ1 i ∆ 0.552 2.151 1.6 hy1 hxi hQ1 i ∆ 1.75 13.59 11.84 hy1 hxi hQ1 i ∆ 11.02 -2.68 -13.70 order Q1 hy2 hxi hQ1 i ∆ 0.163 2.344 2.18 hy2 hxi hQ1 i ∆ -0.34 14.63 14.97 hy2 hxi hQ1 i ∆ 2.13 -1.8 -3.93 order Q1 1.913 hy3 1.877 hxi hQ1 i -0.036 order Q1 ∆ 1.913 13.58 11.67 hy3 hxi hQ1 i ∆ 0.74 -0.41 -1.15 order Q1 Policy Value 240.51 242.48 243.01 242.47 244.98 244.45 244.97 246.94 Tradeoff to Q1 none y1 y2 y3 y1 &y2 y1 &y3 y2 &y3 y1 &y2 &y3 hy3 hxi hQ1 i ∆ Tradeoff to Q1 none y1 y2 y3 y1 &y2 y1 &y3 y2 &y3 y1 &y2 &y3 Policy Value 739.96 740.22 740.32 739.95 740.58 740.22 740.31 740.58 Tradeoff to Q1 none y1 y2 y3 y1 &y2 y1 &y3 y2 &y3 y1 &y2 &y3 Policy Value 118.33 115.88 117.63 118.13 115.28 115.53 117.4 114.93 Table 2.11: Details of Policy Iteration For the example in the left of Table 2.11, Run 1, policy iteration chooses the policy that only makes the tradeoff for the small delivery size at y3 , for the example in the middle, Run 4, policy iteration chooses never to make the tradeoff. In the third example, policy iteration indicates the best policy makes the tradeoff whenever possible — again supported by calculating the value of all policies. These policies are all demonstrated to be optimal in the lower rows. 41 Chapter 2. Optimal Inventory Replenishment with Two Delivery Sizes 2.5 Further Extensions One natural extension is to the case where Q2 is not restricted to 2Q1 . In this case, if only Q2 were ordered, the set of Q2 integer minimizers of G(·) would be the targets and form the band (R, R + Q2 ]. The long run average cost for this policy is easily computed as the stationary distribution is Q12 . With no differences in unit costs between delivery sizes, the lowest cost policy will order into the set of P Q1 integer minimizers of G(·). If Q2 is an 1 1 integer multiple of Q1 then Q1 Q i=1 G(xi ) still is a lower bound. If Q2 is not an integer multiple of Q1 then it is not guaranteed that the stationary distribution of the policy ordering to the set of Q1 integer minimizers of G(·) will be constant over all states. For any two delivery size sets Q1 < Q2 and unit costs ci such that it is reasonable to order both, there are 2Q2 −Q1 candidates for optimal policy. For each of the Q2 − Q1 integer minimizers of G(·) {xQ1 +1 , . . . , xQ2 } a decision is made to order Q1 to reach an inventory position with lower value for G(·) at the cost of an incremental delivery size. Starting with a policy of ordering only the large delivery size, and a set of target states including • the Q2 integer minimizers of G(·) • the Q2 − Q1 states reached by ordering a multiple of Q1 . These states will be a subset of the above group, but are differentiated with the notations < Q1 > denoting that a small delivery size must be ordered to reach them. The same policy iteration technique can be used to identify the optimal policy. Table 2.12 shows the results of some runs with two delivery sizes, Q2 6= 2Q1 . These runs illustrate the extension to arbitrary delivery sizes and in these examples the Unichain Policy Iteration converged to the optimal solution within two iterations. 2.6 Conclusions In this work we studied the problem of optimal inventory policy when orders are made up of a combination of two delivery sizes. We found that the optimal solution was highly dependent on the demand distribution, and so optimal solutions do not in general have a simple form. Our analysis 42 Chapter 2. Optimal Inventory Replenishment with Two Delivery Sizes shows how to identify the candidate policies and provides easily calculated upper and lower bounds on the optimal solution. Finally, we demonstrate that Unichain Policy Iteration Algorithm can leverage the structure of the problem to reduce the policy improvement step to a simple subtraction for each of the states where the option to order a smaller delivery size exists. In the runs in the numerical study, Unichain Policy Iteration Algorithm converged with at most 4 iterations. 2.7 Acknowledgments Research supported in part by research grants from the Natural Sciences and Engineering Research Council of Canada (NSERC). We thank Jeannette Song, Mahesh Nagarajan, Tim Huh and Eric Cope for their valuable help. 43 Chapter 2. Optimal Inventory Replenishment with Two Delivery Sizes Q1 Q2 c1 c2 2 6 20 19 2 5 20 19 2 7 15 12 2 8 20 15 7 21 20 19 h b 30 30 30 30 22 33 20 20 30 30 Demand Uniform (65,135) (90,110) (77,123) (85,115) (90,110) g∗ UB Q1 UB Q2 UB LB 533.50 534.61 533.73 533.73 532.61 158.48 159.86 160.00 159.86 157.86 312.45 316.22 312.45 312.45 310.22 158.39 165.16 158.39 158.39 155.16 167.52 169.86 209.52 169.86 162.86 UB-LB U B−LB g∗ U B−LB LB 1.13 0.21% 0.21% 2.00 1.26% 1.27% 2.22 0.71% 0.72% 3.23 2.04% 2.08% 7.00 4.18% 4.3% gP ∗ g P1 P g 1 − g∗ # It 533.50 533.50 0 2 158.48 158.48 0 2 312.45 312.45 0 1 158.39 158.39 0 1 167.52 167.52 0 2 0.016 0 0.016 0.016 0.016 0.031 0.015 125.56 Solve Time (s) UPIA 0 all 2Q2 −Q1 0 g∗ UB Q2 gP ∗ g P1 # It UPIA average cost of the optimal policy obtained through enumeration PQ2 1 i=1 G(xi ) Q2 average cost of the policy found through the Unichain Policy Algorithm average cost of the policy found from one iteration of the Unichain Policy Iteration Algorithm. Note: The 0th iteration will have value UB Q2 number of iterations of the Unichain Polity Iteration Algorithm to convergence. Note: Ordering all Q2 is considered the 0th iteration. One additional iteration is required to confirm convergence. Unichain Policy Iteration Algorithm Table 2.12: Runs with Two delivery sizes, Q2 6= 2Q1 44 Chapter 3 Sequence Optimization in Block Cave Mining 3.1 Introduction When ore bodies such as diamonds, gold and copper lie far below the surface, traditional open pit mining techniques are neither cost effective nor environmentally friendly. Instead, deep ore bodies can be excavated using a technique called block cave mining, Gertsch (1998). Vertical shafts are sunk below the ore body and a network of horizontal tunnels are dug as a platform for the mining operation. Blasting is done above the tunnels causing the rock to fall in a controlled manner into the tunnels. The cave formed by the initial blasting weakens the rock above, and the rock continues to fall. Figure 3.1 shows a cross section of a block cave mine. It is well known that the traditional open pit mine design problem can be represented in a network flow framework, Ahuja et al (1993), and solved using the Lerchs Grossman algorithm, Lerchs & Grossman (1965) or mincut/ max flow algorithms, i.e., Yegulalp and Arias (1992), Giannini et al (1991). The ore body is divided into blocks via a 3-dimensional grid, each with an estimated value based on the ore composition. The objective is to determine the blocks to extract to maximize the total value, subject to physical constraints. These constraints include having to remove the blocks above a desired block and maintaining a pit shape that contains stable slopes. At a first approximation, a block cave mine can be thought of as an inverted open pit mine. The ore body is divided into blocks, with each vertical column of blocks accessible from a draw point along the horizontal tunnels. Ore at the bottom of a column must be extracted before ore at the top can be reached, and in order to keep the caves stable, neighbouring draw points must extract at a similar rate. Block cave mining differs significantly from open-pit mining in two ways. First, the amount extracted from each draw point in each time period must 45 Chapter 3. Sequence Optimization in Block Cave Mining Figure 3.1: A cross section of the underground workings at the Dutoitspan mine, Kimberley, South Africa. Adapted from De Beers. (American Museum of Natural History (2008)) exceed a minimum amount in order to keep the rock flowing, and must not exceed a maximum amount for safety considerations. Secondly, having the draw point locations tied to both the tunnel structure and the blasting procedures used imposes additional constraints on the order in which the draw points are brought into production. Typically a block cave mine will consist of 200 to 2000 draw points and will be mined for 10 to 20 years. Gemcom (2006), our industrial partner in this work, provides specialised mining productivity solutions that help their clients plan, monitor and improve their mining operations. They have worked with every major mining company in the world and have extensive experience planning block cave mines. To support this work they have developed PC-BC, a program for designing and evaluating block cave mining 46 Chapter 3. Sequence Optimization in Block Cave Mining operations, Diering (2000). Given the layout of draw points, the order in which they should be opened, and detailed information on the ore body such as grade and grade distribution which would be important to define the NPV values mentioned later, the program produces a production plan that tries to maximize the net present value of the material extracted from the mine. The model takes into consideration both vertical mixing, within the rubble in a draw point, and horizontal mixing, across draw points. This work is designed to assist in providing a required input to the PCBC program; the order in which the draw points are opened. This step, known as sequence optimization, is currently a manual process that relies on the expertise of the mine planner. The sequence lists which draw points will be opened in each year of operation and is used by PC-BC to plan each year of mine operation. The current process first manually develops a sequence, then evaluates the sequence using PC-BC, then tries to improve the sequence, etc. As it is practical to evaluate only a small subset of sequences, planners do not know how much value they are missing out on, and how much effort should be invested in the iterative process. The objective of this work is to find an optimal opening sequence in an automated manner. 3.1.1 Current Practice Currently sequence optimization is done by an experienced person on a trial and error basis. It is assumed once a draw point is opened ore will be withdrawn at a constant rate and a lifetime can be calculated for each draw point. A value can also be assigned to each draw point which takes into consideration the estimated lifetime and the ore profile. The feasible opening sequences are created by incorporating draw point value, knowledge of physical constraints and experience, then evaluated. A series of opening sequences are evaluated and the results used to try to create an opening sequence that improves on the highest achieved net present value. Unless every possible feasible opening sequence is evaluated, this trial and error approach can not be guaranteed to deliver the plan with the highest net present value. It is also not possible to estimate or put bounds on the highest net present value for a given ore body using the trial and error approach. 3.1.2 Tunnel Details The procedures used in the mine to develop and open the draw points constrains the order in which they can be opened. We now describe the process 47 Chapter 3. Sequence Optimization in Block Cave Mining to aid in the development of the algebraic constraints. A block cave mine consists of the tunnels and access routes in the mining footprint, the “draw bells” where the blasting is done and the “draw points” where the ore is removed by the trucks. As the broken rock is removed from draw points, the rock above the draw bell falls into the open space creating more broken rock to be removed. All the blasting, rock removal and other mining activities take place at the mining footprint below the ore body that is being mined. See Figure 3.2 Figure 3.2: Cross Section Picture of Block Cave Mine The draw points are arranged on a network of parallel, horizontal tunnels at the bottom of the mining area. This is known as the Extraction Level. Between the tunnels a series of draw bells are constructed, each accessed by two draw points. Figure 3.3 shows draw point and draw bell orientation to tunnels. The rock will fall from above into the draw bells, then be removed through the draw points and transported along the tunnels to the conveyor system which takes the rock out of the mine. Approximately 18 meters (m) above the Extraction Level is the Undercut Level where the initial blasting is done. The steps of Undercut Level development are summarized in Figure 3.4. Horizontal Undercut Tunnels are dug above and parallel to the tunnels in the Extraction Level. Depending on the mine layout, the Undercut Tunnels are either directly aligned with the Extraction Tunnels, or at some consistent offset. When a draw point is 48 Chapter 3. Sequence Optimization in Block Cave Mining Figure 3.3: Plan of Extraction Level of Block Cave Mining Operation scheduled to start producing, blasting is done to free the rock above so it can fall into the draw bell. Blast holes are drilled in the walls and ceiling of the Undercut Tunnel above the draw bell, making a ring or fan shape around the tunnel. The holes are packed with explosives and the blasting causes the rock to fall down to the draw bell below. See Figure 3.5. For the purposes of this work, the constraints of the blasting process can be summarized as follows: blasting starts at one point in the undercut tunnel and works along the tunnel, in both directions, from that point. Since the draw point opening is synonymous with blasting, the draw point opening sequence must also progress along the tunnel from the initial start point. At any time, the opened draw points must form a contiguous group with no gaps between opened draw points. This is illustrated in Figure 3.6 which depicts three single tunnels and the draw points along them. For each example, the number in each square refers to the period in which the draw 49 Chapter 3. Sequence Optimization in Block Cave Mining Figure 3.4: Steps in Undercut Development. Adapted from Resolution Copper Mining (2009) Figure 3.5: Example of Conventional Undercut Layout Showing Blasting Pattern. Barber (2000) point is opened. In the first example, a draw point in the middle region is the initial start point and the draw points to the right are opened until the end of the tunnel is reached, then those to the left of the start point are opened. In the second example, the opened draw points progress in both directions. Both these examples form contiguous groups and meet the blasting constraints. The third example, at the bottom of the figure, breaks contiguity in the third period by opening a lone draw point at the left hand 50 Chapter 3. Sequence Optimization in Block Cave Mining end of the tunnel. Figure 3.6: Three Examples of Single Tunnel Opening Sequences. Each square represents a draw point on a tunnel and the number in it the period in which the draw point is opened. The top two maintain a contiguous group of opened draw points, the bottom one does not. In addition to the constraints on the order in which draw points within a tunnel can be opened, there are also some across-tunnel constraints. In general, the mining operation is thought of as a cave and the structural conditions of the cave roof govern the falling rock and the propagation of the cave upwards. The objective is a contiguous cave that grows out from a starting point. The shape should be convex so as not to leave isolated pockets of rock. A diamond shape cave is commonly created as this has resulted in good rock behavior. As a result, the opening sequence of draw points must consider what is happening in adjacent tunnels. In the PC-BC program this is covered under neighbour constraints which prevent one draw point from producing much more than its within-tunnel and adjacent-tunnel neighbours. To create/preserve the cave structure, some coordination is needed between the tunnels. When work begins in an adjacent tunnel, it will be next to rings that have already been developed, see Figure 3.7. The difference in the number of blasting rings between adjacent tunnels is controlled so that the lead or lag does not exceed an upper limit. Generally it is acceptable to have a lead of one draw bell (approximately 7 rings), but a lead of two draw bells (15 rings) may not be acceptable. While the cave must grow in some sort of contiguous fashion, it is not always a constraint that there be only one cave. In rare cases such as very 51 Chapter 3. Sequence Optimization in Block Cave Mining Figure 3.7: Plan of Undercut Level of Block Cave Mining Operation large mines, mining can begin in two separate tunnels, and the caves merge together over time. In practice, a cave is initiated by opening multiple adjacent draw points in order to create the conditions that cause the rock to drop and it is not common to have multiple caves forming. 3.1.3 Cave Shape In current practice, the draw point opening pattern grows in a diamond shape, aligned with the tunnels as illustrated in Figure 3.8. A chevron pattern is also common. This diamond/chevron pattern conserves the within-tunnel constraint of one starting point, and also links the adjacent tunnels together. Mining engineers prefer the diamond shape opening patterns over other shapes for operational reasons. One reason is the known as the “hydraulic radius”. The hydraulic radius is the ratio of the opening area to perimeter, 52 Chapter 3. Sequence Optimization in Block Cave Mining Figure 3.8: Typical Diamond Pattern of Opened Draw Points. The solid dark draw points are opened in the first period, the white ones in the second period and the hatched ones in the third period. Diamond shapes have been drawn around the draw points opened in each period. and gives an indication of how well the rock will break apart and fall into the draw points. Diamond shaped openings result a hydraulic radius in the desired range for good caving. The diamond shapes also ensure that mining progresses into new rock in a straight-edged wedge shape. While the interior angle of this wedge can vary, in practice values below 40 degrees have not been used. 53 Chapter 3. Sequence Optimization in Block Cave Mining 3.1.4 Other Considerations There are other considerations that can arise. At this time they are outside of the scope of this thesis, but will be listed here. Cave Limitations Some mines may have additional constraints imposed that are unique to that mine. For example, there may be surface activity above part of the ore body that must be completed before the block cave mining can begin. As a result, the options for the starting point would be restricted to a smaller subset of draw points not below the surface activity. There may also be additional data available from sampling that indicates additional information or constraints. This could also take the form of setting the starting draw points for the mine. Grade Limitations Some ore bodies are contaminated with other minerals. While this is usually not a concern with mine planning, it may be a concern in the final product leaving the downstream processing plant. A planning schedule may have to take into account that on average the final product may have limits on the contaminant content. 3.1.5 Assumptions A starting assumption in this work is that we can estimate a lifetime for each draw point by assuming a constant draw rate. The draw rate is the amount of broken rock that is removed in each period from each draw point. In reality this is not the case as the PC-BC program helps to optimize the draw rates of each draw point. In addition, production rate increases over time since the rock removed generally becomes finer over time. That is, when a draw point is first opened, the rock removed is large, coarse pieces, but over time the rock pieces become smaller. 3.1.6 Previous work Models have been created of many stages of the mining process. Newman (2000) reviews mine planning applications of operations research and notes the literature dates back to the 1960s. The work had been conducted on surface and underground mine planning problems that include equipment selection and long- and short-term production scheduling. Work in block cave 54 Chapter 3. Sequence Optimization in Block Cave Mining mining includes Diering (2006), who models the rock flows from individual draw point, taking into account lateral mixing. Rubio (2002) maximizes NPV in a block cave mine using draw rate as one constraint. The cave area is set by specifying the area of a contiguous cave. Weintraub (2008) creates a global model of all the copper mines owned by CODELCO, the Chilean state copper mine. All these models assume that the opening sequence is known and use it as a constraint that is determined prior to running their models. 55 Chapter 3. Sequence Optimization in Block Cave Mining 3.2 Model Framework The problem that we consider is the following. Given the total net present value (NPV) and expected lifetime of each draw point, determine which draw points should open in each period to maximize the total NPV of the mine subject to the physical and operational constraints. Due to discounting, the total NPV of the mine is a function of the opening sequence. Intuitively, opening draw points with higher total NPV first increases the total NPV of the mine over starting with draw points of lower value. In the absence of any constraints, opening all positively valued draw points in the first period, and never opening the negatively valued ones, maximises the total NPV of the mine. The total NPV of draw point i is the net present value if it was started at time 1 and ore was extracted at a constant rate until it was exhausted, and is represented by wi . “Total” emphasises that this value represents the cumulative value over all the periods from start to completion. The time to exhaustion/completion is its duration, pi . A discount rate, λ, is used to calculate the net present value. Typically, a discount rate of between 7 to 12 % per year has been used, but values as high as 15% have been used, Diering (2006). To clarify, the period when the first material is removed from a draw point is the period in which a draw point is started. In this model, the draw point will be active for the next pi periods including the period in which it was started. A draw point will be called open or opened in the period it is started and each following period until the end of the time horizon. In this model we are interested in how many draw points are active in each period, and the contiguity of opened draw points. 3.2.1 Data Data : pi duration of draw point i wi total NPV value of draw point i if started immediately λ discount rate, 0 ≤ λ ≤ 1. m capacity constraint on number of draw points active at any time 3.2.2 Decision Variables Decision Variables Si,t Binary decision variable modelling the start date for draw point i, over the lifetime of the mine, 1 ≤ t ≤ T . First 56 Chapter 3. Sequence Optimization in Block Cave Mining period is 1. Si,t = 1 if draw point i is started in period t, 0 otherwise Typically the plans are done with a period of length one year or six months. We will assume one period is a year long, but essentially the same formulation works for any period length. 3.2.3 Objective The objective function is the sum of draw point values, discounted to the start date. X X max V alue = wi Si,t e−λ(t−1) i t for some discount rate λ. 3.2.4 Constraints Constraints : Start Once To control the openings, each draw point can be started at most once X Si,t ≤ 1 ∀ draw points i (3.1) t Note: If t X Si,u = 1 then the draw point was started at some u=1 period {1, . . . , t} and is considered open. Otherwise, this sum equals zero. Global Capacity We assume a global capacity constraint, whereby at most m draw points can be active at any time. In general, this capacity is linked to the capacity of the downstream ore processing plant. X t X i u=t−pi +1 Si,u ≤ m ∀ periods t (3.2) Due to the high cost of developing each draw point, approximately $300,000, production will typically ramp up over the first few years. This can be handled by imposing a maximum number, mt to open each year t, i.e., m1 = 50, m2 = 120, mt = 200 for t ≥ 3. 57 Chapter 3. Sequence Optimization in Block Cave Mining Ideally, the ramping rate would be dictated by the model and be a function of the opening cost and the initial available money. A budget constraint for each period would ensure that the draw point opening expenditure did not exceed the available money plus any revenue generated. In reality, the opening cost is incurred at time zero, and the revenues start to flow back at time 1. In the current structure of the model, the revenue has been lumped to a NPV to the end of the first period after opening – so all the revenues would occur at time 1 instead of distributed over the assumed life of the draw point. To approximate, we could take the total NPV and distribute it evenly over the life of the draw point. If more detailed information of the revenue variation over time was available it could also be included. Tunnel Development The mechanics of Block Cave Mining impose two main tunnel development constraints, one within-tunnel and the other across-tunnels. Additional Constraints : Within-Tunnel Contiguity Within a tunnel, the draw point opening starts at one point in a tunnel, then progresses successively in both directions along the tunnel. There are no gaps left between open draw points. Across-Tunnel Contiguity Between tunnels, the opened draw points must form some sort of contiguous cave. Thus when the first draw point in a tunnel is opened, it must be next to opened draw points in an adjacent tunnel. These constraints are difficult to write in algebraic form. Several options will be presented in this work. 3.2.5 Unconstrained Sequence Optimization as a Draw Point Scheduling Model The unconstrained version has the Start Once (3.1) and Global Capacity (3.2) constraints, but none of the cave defining constraints. A parallel can be drawn between the unconstrained version of this problem and machine scheduling problems. The unconstrained version of sequence optimization can be viewed as scheduling the n draw points as “jobs” on m “parallel, identical machines”. 58 Chapter 3. Sequence Optimization in Block Cave Mining P This problem is well studied in a non-discounted setting as pk wj Cj . See Pinedo (1995) for notation explanation and introduction to scheduling problems. The objective is to schedule jobs, all available from the first time period, to minimize the weighted sum of completion times. For the case of a single machine, the optimal solution can be found using Smith’s Rule or the Weighted Shortest Processing Time (WSPT) rule: sequence the jobs in nondecreasing order of their ratio of processing time to weight, wpii . When there is more than one machine, the problem is NP Hard. Results from this problem show that once jobs have been allocated to a machine, they are scheduled using WSPT on that machine. We look at the discounted problem under two conditions, one of uniform draw point duration, and the other when there is no assumption of uniformity. Uniform Draw Point Duration pi = p In the case of uniform draw point duration of p periods, the problem reduces to deciding which m draw points will open in periods 1, p + 1, 2p + 1 . . . In the absence of discounting, it does not matter in what order they are opened. With discounting, it is intuitive that to maximize NPV, the most valuable draw points should be opened first. We now show that a greedy algorithm on wi , the total value of the draw point, is optimal. Proposition 3.1. When scheduling jobs of uniform duration p on m machines, greedily choosing jobs in non-increasing order of wi , the total value of the job, maximizes NPV, with discount factor λ. Proof. It is clear that for non-negative wi it is optimal to start each job as soon as possible, so it would never be optimal not to start m jobs in the first period, or to ever have less than m jobs active in any given period (except at the end if there are no remaining jobs). We order the jobs by value wi , in non-increasing order, w1 ≥ w2 ≥ . . . ≥ wJP . In period 1 start the m most valuable jobs, i.e., w1 , . . . , wm for a value of m p+1, the next m most valuable jobs are started k=1 wk , then, P2min period−λp for a value of k=m+1 wk e , etc. Consider any alternate sequence of the same job values (s1 , s2 , . . . , sn ), and the resulting schedule defined similarly, Pby starting the first m jobs, i.e., those with values s1 , . . . , sm for a valueP of m k=1 sk , then, in period p+1, the −λp , etc. next m jobs are started for a value of 2m k=m+1 sk e We will switch the first job in the alternate sequence, that does not match with greedy sequence. Denote this as job r. We have that si = 59 Chapter 3. Sequence Optimization in Block Cave Mining wi for i = 1 . . . r − 1 and sr < wr . Let sα = wr and we shall switch the positions of sr and sα . If sr and sα are originally in the same period, then there is no change to the value of the sequence by making the switch. If sr was originally in period q + 1 and sα was originally in a later period q + β + 1, β > 0, then the value of the sequence will change by −sr e−λq + sr e−λ(q+β) − wr e−λ(q+β) + wr e−λq . Reorganizing, the change in objective value is (wr − sr )e−λq + (sr − wr )e−λ(q+β) = (wr − sr )(e−λq − e−λ(q+β) ). For λ ≥ 0 we conclude that this switch is an improvement to the value of the sequence and the thus the algorithm gives us an optimal solution. Non-uniform Draw Point Duration This problem is more difficult. In this form P the problem is similar to the NP-hard, Garey (1979), scheduling problem pk wj Cj , which minimizes the sum of delay penalties calculated as the time to completion (Cj ) times a per period penalty (wj ). To see the similarity, assign each draw point i value wdelayi , the total value of the draw point at the end of its lifetime pi , given that P the draw point was started at period 1. Thus our problem becomes maxi wdelayi e−λCi we P . If instead T −Ci i had a linear discount T −C the problem would become max = wdelay i T T P P P i Ci wdelayi − wdelayi T , and the similarity to minj wj Cj becomes clear. However, the discounting in our problem is not linear so a different solution is needed. Let’s start with a single machine and a modified greedy algorithm. Proposition 3.2. When scheduling jobs of non-uniform duration pi , each with total value wi , on a single machine, greedily choosing jobs in nonwi increasing order of 1−e λpi , maximizes NPV, with discount factor λ. wi Proof. Order the jobs by ratio of 1−e λpi in non-increasing order and run the jobs in that order. Again it is clear that the machine should always be busy for non-negative wi . The total objective value will be w1 + w2 e−λp1 + w3 e−λ(p1 +p2 ) + . . . + wJ e−λ(p1 +...+pJ−1 ) . Consider any alternate sequence of the same ratios ( 1−es1λd1 , 1−es2λd2 , . . . , 1−esnλdn ). Running the jobs in this or- der will result in an objective value of sw1 + s2 e−λd1 + s3 e−λ(d1 +d2 ) + . . . + sJ e−λ(d1 +...+dJ−1 ) . Let r be the position of the first ratio in the alternate sewi quence that differs from the greedy sequence. We have that 1−esiλdi = 1−e λpi w sr sα wr i for i = 1 . . . r−1 and 1−eλdr < 1−eλpi . Let 1−eλdα = 1−eλpr . Unless both jobs have the same duration, or are immediately adjacent to each other, i.e., α = r+1, more than these two terms in the objective function will change. Let us 60 Chapter 3. Sequence Optimization in Block Cave Mining first consider the case of adjacent jobs, α = r+1. The objective function will change by −sr e−λ(p1 +...+pr−1 ) −wr e−λ(p1 +...+pr−1 +pα ) +sr e−λ(p1 +...+pr−1 +pr ) + wr e−λ(p1 +...+pr−1 ) = e−λ(p1 +...+pr−1 ) (wr + sr e−λpr − sr − wr e−λpα ) = e−λ(p1 +...+pr−1 ) (wr (1 − e−λpα ) − sr (1 − e−λpr )) > 0. For λ > 0 this is an improvement to the value of the sequence. Next look at the case where α = r + 2. This switch can be broken down into two adjacent job switches. First change the order of job α = r + 2 and job r + 1. We have shown above that this switch improves the objective value. Now switch jobs α and r. wr+1 wr wα Since 1−e ≥ 1−e λpr ≥ λpα we can again use the argument above to 1−eλpr+1 show that the switch improves the value of the alternate sequence. Since there is improvement from both of the two adjacent switches we conclude there is improvement from the switching of jobs r and α. A similar argument can be made for any switch of jobs and we conclude that the ordering wi maximises the objective value. 1−eλpi Next, move on to scheduling on multiple machines. We have shown above that once the jobs are allocated to each machine, there is a greedy algorithm that is optimal. It is the allocation of jobs to machines that is the problem. It can be shown by the following counter-example that using WSPT ordering to assign to multiple machines is not optimal when minimizing weighted completion times. Example: Schedule the 8 jobs shown below on three parallel machines; the jobs have been sorted in nonincreasing wpii wi Job Number wi pi pi 1 10 2 5 2 4 1 4 3 9 3 3 4 3 1 3 5 8 4 2 6 9 5 1.8 7 5 3 1.667 8 7 5 1.4 A scheduling of the jobs as ordered above is as follows. In period 1 start job 1 on machine 1, job 2 on machine 2 and job 3 on machine 3. In period 2 job 2 will be completed with a value of 2, and job 4 is started on machine 2. In period 3 jobs 1, value 20, and 4, value 6, are completed and job 5 is started on machine 1 and job 6 started on machine 2. In period 4 job 3 is completed, value 27, and job 7 started on machine 3. In period 7 job 5 is completed, value 48, and job 8 started on machine 1. Finally, job 7 is 61 Chapter 3. Sequence Optimization in Block Cave Mining completed in period 6 for a value of 30, job 6 is completed in period 7 for a value of 63 and job 8 is completed in period 11 for a value of 77. The total weighted completion time is 275. The job schedule and weighted completion times are illustrated below. Period 1 2 3 4 5 6 7 8 9 10 11 Machine 1 job 1 1 5 5 5 5 8 8 8 8 8 total wi C i 20 48 77 145 Machine 2 2 4 6 6 6 6 6 wi C i 4 6 63 73 Machine 3 3 3 3 7 7 7 wi C i 27 30 57 275 Now let’s switch the scheduled order of jobs 6 and 7. There is no change in the schedule until period 3. Instead of starting job 6 on machine 2, job 7 is started in its place. In period 4 job 3 is completed and job 6 is started on machine 3. In period 5 job 7 is completed with a new, lower, value of 25, and job 8 is started on machine 2. Next job 6 completes in period 8 for an increased value of 72 and job 8 is completed in period 10 for a reduced value of 70. Overall the weighted completion time is reduced to 272, an improvement over the original schedule as illustrated below. Period 1 2 3 4 5 6 7 8 9 10 11 Machine 1 job 1 1 5 5 5 5 total wi C i 20 48 68 Machine 2 2 4 7 7 7 8 8 8 8 8 wi C i 4 6 25 70 105 Machine 3 3 3 3 6 6 6 6 6 wi C i 27 72 99 272 This example demonstrates that using WSPT is not optimal for scheduling jobs, with non-uniform duration, on multiple machines when minimizing weighted completion time. The same example can be used to show that a greedy algorithm on the i ratio 1−ew−λp does not maximize total NPV when λ = 0.1. Switching the i order of jobs 6 and 7 increases the NPV from 47.180 to 47.272 as shown below. 62 Chapter 3. Sequence Optimization in Block Cave Mining Job Number wi pi wi pi wi 1−e−0.1∗pi 1 2 3 4 5 6 7 8 10 4 9 3 8 9 5 7 2 1 3 1 4 5 3 5 5 4 3 3 2 1.8 1.667 1.4 55.167 42.033 34.725 31.525 24.266 22.873 19.291 17.790 total Value with greedy 10 4 9 2.175 6.550 7.369 3.704 3.842 47.180 Value with jobs 6 and 7 switched 10 4 9 2.175 6.550 6.667 4.094 4.246 47.272 The forgoing shows that even simple variations of this problem are challenging. This motivates investigating the use of general integer programming methods. 63 Chapter 3. Sequence Optimization in Block Cave Mining 3.3 Single Tunnel We begin by capturing the physical neighbour constraints on opening the draw points begun within a single tunnel, as the Within-Tunnel Contiguity constraint is a more manageable starting point than the combined WithinTunnel Contiguity and Across-Tunnel Contiguity constraints of the twodimensional problem. We assume that a tunnel consist of n draw points and they are numbered in increasing order along the tunnel, so the neighbours of draw point i are i − 1 and i + 1. The draw points at the end of the tunnel, 1 and n, each only have a single neighbour, 2 and n − 1 respectively. The problem with formulating the single tunnel problem is to write algebraic constraints that result in a single, contiguous cave in each period. The cave includes previously opened draw points as well as any that may open in that period. A simple approach is to say that a given draw point can be started in a period if either of its immediate neighbours was started in the current or any previous period, i.e., if considered open in the current period. For the draw points at the ends of the tunnel there is only one immediate neighbour. This constraint works if the global capacity constraint is m = 1, but if multiple draw points can be started in a period then they can be next to each other and act as open neighbours for each other, since the sum is up to and including the current period. This can be avoided by only counting draw points opened in previous periods, i.e., sum to previous period. Unfortunately, this would only allow at most two draw points to start each period, one on each end as no others would have an neighbour open in an earlier period. Two formulations were developed for the single tunnel problem. The alternating constraints formulation, see Section 3.3.1, sums up each odd subset of draw points in a tunnel. If there is only one cave along the tunnel, each sum will be at most one. If there are multiple tunnels, at least one sum will exceed one. There are a large number of these constraints, but we show a subset that is sufficient if binary variables are used. In the extended formulation, see Section 3.3.2, an additional binary decision variable is added that compares each draw point to its neighbour and indicates if that draw point is the first in a new cave. By allowing only one transition to a new cave, a single, contiguous cave is achieved. 64 Chapter 3. Sequence Optimization in Block Cave Mining 3.3.1 Alternating Constraints Formulation Single Tunnel: Alternating Constraints This section describes a model that uses a set of constraints we have called alternating constraints to achieve the Within-Tunnel Contiguity constraint. Mathieu Van Vyve (2005) suggested the following set of constraints to describe the convex hull of the Within-Tunnel Contiguity constraint. Decision Variables xi,t ≥ 0 xi,t = Pt u=1Si,u point is open in period t. xi,t = and indicates if the draw 1 if draw point i is open 0 otherwise and the following alternating inequalities: for any increasing sequence i(1) < i(2) < . . . < i(k) ∈ {1, . . . , n} with k positive, odd integer k X (−1)j+1 xi(j),t ≤ 1. j=1 The k = 1 constraints, xi,t ≤ 1 are satisfied by the definition that 0 ≤ Si,t ≤ 1 and the Start Once (3.1) constraint. Van Vyve further suggested that if Si,t is forced to binary, that k = 3 is sufficient, see Theorem 3.1 below, but the higher values of k are necessary for the linear relaxation. This is further discussed below. For a single tunnel of n draw points there are n3 = (n)(n−1)(n−2) ,k=3 6 n constraints plus j , k = j constraints for j = 5, 7, . . . , n if integrality of Si,t is not specified, in each period. Theorem 3.1. The k = 3 alternating constraints ensure the Within-Tunnel Contiguity constraint if Si,t are binary. Let C = {x ∈ {0, 1}n | x is contiguous } and let A = {x ∈ {0, 1}n | xi,t − xj,t + xl,t ≤ 1 ∀ i < j < l ∀ t}. We claim that A = C. Proof. It suffices to show C ⊆ A and A ⊆ C. Start with the first, that every consecutive vector satisfies all the (k = 3) alternating constraints. Given x ∈ C, for any pair xi,t , xj,t i < j there are 4 cases (xi,t = 0, xj,t = 0) Then for any xl,t , l > j the alternating sum will be ≤ 1 (xi,t = 0, xj,t = 1) Then for any xl,t , l > j the alternating sum will be ≤ 0 65 Chapter 3. Sequence Optimization in Block Cave Mining (xi,t = 1, xj,t = 1) Then for any xl,t , l > j the alternating sum will be ≤ 1 (xi,t = 1, xj,t = 0) Then since x ∈ C, for any l > j xl,t = 0 and the alternating sum will be zero. Thus x ∈ A and C ⊆ A. For the reverse, take x ∈ A. By definition, x satisfies all of the k = 3 alternating constraints. For x ∈ {0, 1}n the only way to fail the alternating constraint is to have the triple {xi,t , xj,t , xl,t } = {1, 0, 1}, ∀ i < j < l. Thus x ∈ A means that for every pair {xi,t , xl,t }, i < l, if xi,t = 1 and xl,t = 1 then each xj,t = 1 ∀ i < j < l. This means that there are no “gaps” in the 1’s, and thus x ∈ C and A ⊆ C. Infact, some of the k = 3 alternating constraints are redundant and not all n3 are required. Lemma 3.2. One sufficient set of k = 3 alternating constraints is xi,t − xi+1,t + xl,t ≤ 1 ∀ i + 1 < l. There are (n−1)(n−2) constraints in each period 2 in this set. Proof. Let B = {x ∈ {0, 1}n | xi,t − xi+1,t + xl,t ≤ 1 ∀ i + 1 < l}. Following the proof for Theorem 3.1, we first recognize that since B is a subset of A, every consecutive vector satisfies the subset of the (k = 3) alternating constraints. For the reverse, take x ∈ B. By definition, x satisfies the subset of the k = 3 alternating constraints. For x ∈ {0, 1}n the only way to fail the alternating constraint is to have the triple {xi,t , xj,t , xl,t } = {1, 0, 1}, ∀ i < j < l. Since x ∈ B, for every pair {xi,t , xl,t }, i < l, if xi,t = 1 and xl,t = 1 then xi+1,t = 1 ∀ i + 1 < l. Thus if xi,t = 1 and xl,t = 1, then xi+1,t = 1. Now we have xi+1,t = 1 and xl,t = 1, so xi+2,t = 1. Thus each pair {xi+u,t , xi+u+1,t } will both be 1 ∀ u < l − 1 − i. This means that there are no “gaps” in the 1’s, and thus x ∈ C and B ⊆ C. Separation of Alternating Constraints There are potentially a huge number of alternating constraints, roughly 2n−1 , the number of odd subsets of {1, . . . , n}.. If these constraints do describe the convex hull of the contiguity constraint, it may be more useful to employ a “Branch and Bound” strategy (Wolsey 1998). This would require 66 Chapter 3. Sequence Optimization in Block Cave Mining a separation algorithm to determine if a produced x vector were in the hull, and if not, which constraints it violated. The following algorithm is a first attempt at a separation algorithm for the k = 3 alternating constraints. The idea is to start with the x(i) entries with the highest value. If adding them together does not get us over the 1 hump, then there is no chance of subtracting something and staying over 1. Algorithm for k = 3 alternating constraints Start with vector x(i), 1 ≤ i ≤ n with i designating the position in the tunnel, i = 1 at one end and i = n at the other. Create d(j), a sorting of x, such that d(j) gives the position in the tunnel of the j th largest value of all x(i), i.e., d(1) = a if x(a) ≥ x(k) ∀ k = 1, . . . , n and d(n) = a if x(a) ≤ x(k) ∀ k = 1, . . . , n Start with the highest two x values, if their sum is not greater than one, stop. If they do, look for an x located between them in the tunnel that, when subtracted from the sum, does not bring the total below 1. If one cannot be found, go on to next highest value and try again. procedure k = 3(x, n) 1: if x(d(1)) + x(d(2)) ≤ 1 then 2: STOP {will never find a triple to violate the constraints.} 3: end if 4: if x(d(1)) + x(d(2)) − x(d(n)) ≤ 1 then 5: STOP {will never find a triple to violate the constraints.} 6: end if 7: for j = 1 : n − 1 do 8: if x(d(j)) + x(d(j + 1)) ≤ 1 then 9: STOP {no more options of finding a triple} 10: end if 11: for l = j + 1 to N do 12: if x(d(j)) + x(d(l)) > 1 then 13: {continue and look for a middle value to complete the triple} 14: for h = min{d(j), d(l)} + 1 to max{d(j), d(l)} − 1 do 15: if x(d(j)) − x(h) + x(d(l)) > 1 then 16: {have violation, record and proceed to next h} 17: end if 18: end for{next h} 19: end if 20: end for{next l} 21: end for{next j} Running Time: The actual running time will be a function of the num67 Chapter 3. Sequence Optimization in Block Cave Mining ber of constraints that are violated. To estimate the worst-case run time, consider that all triples are checked, so the run time is O(n3 ) as there are n3 triples. Related Work In the unit commitment problem Rajan (2005), machines are scheduled to meet known future demand at minimum cost. For some machines, for example electric generators, there are requirements on the minimum time that they must be run, and the minimum time that they must be down before restarting. If the operating status of a generator is modelled using a binary variable for each period, then the problem is to create sequences of ones and zeros that meet the demand and are long enough to meet the minimum up and down times. Lee (2003) developed a set of constraints for the relaxation of this problem which they called the “alternating up inequalities” (for the time periods the machine is running) and “alternating down inequalities” (for the periods the machine is down) and presented a separation algorithm that runs in linear time. Rajan (2005) studied an extension of this problem in which there are start-up and shut-down costs. In order to track these costs, an additional binary variable that indicates if a machine is started up in a given period is needed. This new variable is then used in a set of “turn on/off inequalities” which limit the number of start-ups/shut-downs in each interval of length minimum up time /minimum down time. This new set of inequalities dominates the “alternating up/down” inequalities, and it is also “much smaller in size”. A linear time separation algorithm is presented for use in a branch-and-cut algorithm. This additional variable and set of inequalities is a more general form of the second formulation we developed for the single tunnel problem, which is presented below. 3.3.2 Extended Formulation The idea behind the extended formulation is based on the spatial setting of the draw points in the tunnel and the requirement that a group of opened draw points must have a starting point, or first opened draw point, as one moves along the tunnel in order of increasing draw point number. The beginning of a set of open draw points is characterized by moving from an unopened draw point to an opened one. Setting a limit of only one such transition results in one contiguous set of open draw points. We have called it an extended formulation as an additional decision variable that marks the draw point where the transition occurs, has been added. 68 Chapter 3. Sequence Optimization in Block Cave Mining Single Period In the single tunnel, single period setting, the problem is to determine which contiguous set of draw points to open. Since there is only one period, either a draw point is open/active in the period or remains unopened; we need not be concerned about draw points that are open/inactive. Data : pi , wi , λ, m n number of draw points in tunnel Decision Variables : Si,t Fi,t Binary decision variable modelling the “first” open draw point in the tunnel in period t, 1 ≤ i ≤ n. 1 if draw point i is the first open draw point along the tunnel in period t, Fi,t = 0 otherwise Then since there is only one period and hence no discounting, the objective is max Pn i=1 wi Si,1 Constraints : Start Once (3.1) Global Capacity (3.2) One First Open at most one first open draw point. n X Fi,1 ≤ 1 (3.3) i=1 Definition of First Open i can only be first open draw point if i is open and i − 1 is unopened, 2 ≤ i ≤ n Fi,1 ≥ Si,1 − Si−1,1 (3.4) Definition of First Open, First Draw point in Tunnel draw point 1 can only be first open draw point if it is open. F1,1 ≥ S1,1 (3.5) 69 Chapter 3. Sequence Optimization in Block Cave Mining First Open Must have Started i can not be first open draw point if unopened, 1 ≤ i ≤ n. Fi,1 ≤ Si,1 (3.6) The new binary variable Fi,1 calculates any difference between two draw points in their cumulative operating history, and indicates the beginning of a contiguous group of open draw points. Beginning is defined in the sense of moving down the tunnel of consecutively numbered draw points starting from the lowest numbered one. Thus, Fi,1 = 0 if both, or neither, draw point i and i − 1 are open, or if i is open, but i − 1 is not. This definition is enforced by the Definition of First Open (3.4) constraint. Definition of First Open, First Draw Point in Tunnel (3.5) is included to indicate when the first open draw point is at the beginning of the tunnel. In order to form a contiguous group, the sum of Fi,1 i = 1 . . . n must be 1 (unless nothing is open and the sum would be 0). This single tunnel, single period problem can be framed as a most profitable path problem with the network shown in Figure 3.9 with one unit available to flow from source to sink. The nodes {1, n} represent the draw points in the tunnel and the nodes {10 , n0 } represent the open draw points. Leaving the source, the unit chooses the path Fs,1 , s ∈ {1, n}, the variable Figure 3.9: Network Representation of Single Tunnel in One Period indicating which draw point is the first open in the tunnel in the first period. Next the unit continues along the horizontal path until the node denoting the last open draw point, t, t ∈ {s, n} is open, before travelling to the sink. As a result, arcs Ss,1 , Ss+1,1 , . . . , yt are chosen and set to 1. This problem has an integral optimal solution when the binary constraints on Fi,1 and Si,1 are relaxed. We can show this by first showing that a simple network flow with integral capacities has an integral optimal 70 Chapter 3. Sequence Optimization in Block Cave Mining solution. Then we show a transformation from the simple network to that in Figure 3.9. Figure 3.10 shows a simple network with a source and sink and n nodes. The flow on each of the arcs are given by the variables below. Figure 3.10: Simple Network let xs,i xi,z xi,j Then the max subject to be the flow from the source to draw point i be the flow from draw point i to the sink be the flow from draw point i to draw point j = i + 1 objective is Pn w i=1 i (xs,i + xi−1,i ) Pn i=1 xs,i ≤1 xs,1 = x1,2 + x1,z xs,i + xi−1,i = xi,i+1 + xi,z 2 ≤ i ≤ n send at most one unit of flow from source flow balance on node 1 flow balance on node i The Integral Flow Theorem, Dantzig (1956), states that if the capacities of a network are integers, then there exists an integral maximum flow. The integrality of the optimal solution can also be shown through unimodularity. A linear programming problem of the form max{cx : Ax ≤ b, x ∈ Rn+ } will have an optimal solution that is integral if the matrix A is totally unimodular Wolsey (1998). There are three conditions that are sufficient to show that matrix A is totally unimodular, and we can show they are met with the problem above. First, all elements aij ∈ {−1, 0, +1}, which is true for this problem. Secondly, each column contains a most two nonzero coefficients. If A is constructed in the order of the constraints listed above, we can see that each xs,i will appear in the first row, then again in row for the balance on node i, meeting the condition. Similarly, xi,z only 71 Chapter 3. Sequence Optimization in Block Cave Mining appears in the row for the balance on node i, meeting the condition. Finally, xi,i+1 appears once in the balance on node i − 1 and once in the balance on node i again meeting the condition. The third condition is that there is a partition (M1 , M2 ) of the set M of that each column j containing P rows suchP two nonzero coefficients satisfies i∈M1 aij − i∈M2 aij = 0 (Wolsey (1998)). If we partition A such that M1 is the constraint on the flow from source, and M2 is the remaining constraints, we can satisfy this third condition and show that A is totally unimodular and the problem will have an integral optimal solution. Next we show the translation from the simple network problem Figure 3.10 to the single tunnel, single period problem Figure 3.9. To show the translation, split each node into two and let Fi,1 = xs,i S1,1 = xs,1 Si,1 = xi−1,i + xs,i the flow between two parts of each node And we get the formulation for the single tunnel above, which will also have an integral Pnoptimal solution. Pn max = i=1 wi (xs,i + xi−1,i ) i=1 wi Si,1 subject to P Pn n becomes i=1 Fi,1 ≤ 1 i=1 xs,i ≤ 1 xs,1 = x1,2 + x1,z S1,1 = S2,1 − F1,1 + x1,z S1,1 ≥ S2,1 − F1,1 F1,1 ≥ S2,1 − S1,1 xs,i + xi−1,i = xi,i+1 + xi,z Fi,1 + Si,1 − Fi,1 = Si+1,1 −Fi+1,1 + xi,z 2 ≤ i ≤ n Fi+1,1 ≥ Si+1,1 − Si,1 2 ≤ i ≤ n While this result is promising for the extended formulation, the addition of the Global Capacity Constraint (3.2) can result in fractional solutions. 72 Chapter 3. Sequence Optimization in Block Cave Mining Multiple Periods In a single tunnel, multiple period setting a contiguous set of draw points is opened in the first period. As each active draw point reaches the end of its duration pi , it becomes open but inactive and an additional draw point may be opened. The newly opened draw points in each period together with the previously opened draw points (both active and inactive) must be part of contiguous set. In each period, the capacity constraint parameter m provides an upper limit on the number of active draw points. The formulation becomes Data : n, pi , wi , λ, m T number of periods in the lifetime of the mine Decision Variables : Si,t , Fi,t Now with multiple periods, the discounting appears in the objective function: Pn PT −λt max i=1 wi t=1 Si,t e And the constraints are extended to include multiple periods Constraints : Start Once (3.1) Global Capacity (3.2) One First Open at most one first open draw point, in each period n X Fi,t ≤ 1 ∀ 1 ≤ t ≤ T (3.7) i=1 Definition of First Open i can only be first open draw point if i is open and i − 1 is unopened, 2 ≤ i ≤ n Fi,t ≥ t X Si,u − u=1 t X Si−1,u ∀ 1 ≤ t ≤ T (3.8) u=1 Definition of First Open, First Draw point in Tunnel draw point 1 can only be first open draw point if it is open. F1,t ≥ t X S1,u ∀ 1 ≤ t ≤ T (3.9) u=1 73 Chapter 3. Sequence Optimization in Block Cave Mining First Open Must have Started i can not be first open draw point if unopened, 1 ≤ i ≤ n. Fi,t ≤ t X Si,u ∀ 1 ≤ t ≤ T (3.10) u=1 74 Chapter 3. Sequence Optimization in Block Cave Mining 3.3.3 Computational Results A series of runs were done to test the practicality of computing with the two formulations. CPLEX 9.1 was used on a computer with an Intel Pentium(R) D CPU 3.20 GHz processor with 1.99 GB of RAM for all computational runs. Number of Variables Number of Constraints Example 1: single tunnel 10 draw points time to solve Example 2 single tunnel 23 draw points time to solve k=3 Alternating Constraints D×T (all binary) approx D2 × T Extended Formulation 300 binary 1,120 constraints 2.23 seconds 300 binary, 300 continuous 650 constraints 1.59 seconds 1,150 binary 11,623 constraints 15,701 seconds 1,150 of each 2,423 constraints 1,513 seconds 2D × T (half are binary) approx 2D × T Table 3.1: Table of Single Tunnel Comparison of Alternating Constraints Formulation and Extended Formulation These runs showed that for small tunnels there was little difference in the computation times in CPLEX between the two formulations. However, when the tunnel got longer, the extended formulation solved much quicker than the alternating constraints. Hence we will use the extended formulation as we move towards multiple tunnel formulations in the next section. 75 Chapter 3. Sequence Optimization in Block Cave Mining 3.4 3.4.1 Multiple Tunnels Introduction/Overview With this understanding of the single tunnel problem, we move to the real mine setting of multiple tunnels. In this setting the idea of one contiguous cave still holds within a tunnel, but is expanded to link together activity within adjacent tunnels. The biggest challenge in this work is the definition of an acceptable cave, and writing the constraints to achieve it. It is not fully clear from discussions with Gemcom exactly just what ‘acceptable’ means here, and thus there is some trade-off between mathematical tractability and goodness of the model. The models fall into three main approaches. The first (see Section 3.4.3) is an expansion of the single tunnel formulations discussed in the previous section. Keeping the Within-Tunnel Contiguity conditions, the indices of the decision variables are expanded to reflect multiple tunnels. In particular, a second dimension approximately perpendicular to the tunnels and called “across-tunnel” is added. The same ideas for contiguity (either alternating constraints or extended formulation) are applied in the across-tunnel direction to satisfy the Across-Tunnel Contiguity constraint. The second model, (see Section 3.4.4) comes from a suggestion by Malkin & Wolsey (2006). They propose a series of models, but realize that in the 2dimensional setting (multiple tunnels) it is hard to guarantee the formation of a single contiguous cave. Instead, multiple contiguous caves could be produced. Their suggestion is a model based on a centre point, where only draw points connected to the centre points by open (active or inactive) draw points can be opened. The third set of models, (see Sections 3.4.5— 3.4.7) is motivated by the desire of practitioners (see e.g., Deiring (2006)) to have a “diamond-shaped” cave. The diamond is defined by its vertices and all draw points within the diamond should be open and will be explained more fully later. Finding the relationship between vertices and the draw points is the challenge of this modelling approach. 3.4.2 Computation and Data Sets In order to develop the models proposed in this work, Gemcom made available four data sets to us. See Figure 3.11 for the mining footprints of the four data sets. In the figures, the circles represent the draw points and the zig-zag lines represent the tunnels. These four data sets include one footprint (data set P2) that is (almost) convex, but the others display the 76 Chapter 3. Sequence Optimization in Block Cave Mining Data Set P1 Data Set P2 Data Set P3 Data Set P4 Figure 3.11: Footprints of Data Sets Provided for Model Development variety of footprint shapes that can be encountered in actual mines. Data set P3 has the additional property of two distinct ore bodies. The data sets also demonstrate the size variations encountered in mine planning. Table 3.2 show that these sets range from 236 draw points up to 3181 draw points. Unfortunately this data is confidential and so is not publicly available. Some of these data sets will be referred to in the model description sections later in this document. Multiple data sets were provided to illustrate some of the possible footprint variations, and to aid in model development. Not all of the data sets were used for all models. The performance of all the models will be compared on two additional data sets in Section 3.5. These data sets were not involved in the model development. 77 Chapter 3. Sequence Optimization in Block Cave Mining Data Data Data Data Set Set Set Set P1 P2 P3 P4 Number of Tunnels 31 14 33 39 Total Number of Draw Points 780 236 3181 2387 Table 3.2: Size of Data Sets Provided for Model Development 3.4.3 Adapting Single Tunnel Formulations This section shows the extension of the single tunnel formulations to cover the mine footprint. Unfortunately, the computational results at the end of this section will show that these are not sufficient to guarantee a single contiguous cave. Basic Model Formulation In moving from the one tunnel setting to multiple tunnels, the challenge is to incorporate the Across-Tunnel Contiguity Constraint. This constraint appeared to have direct parallels to the Within-Tunnel Contiguity Constraint in that the cave had to start in one tunnel or a number of parallel tunnels, and gradually grow outwards leaving no gaps. From this idea the single tunnel models were adapted by writing contiguity constraints similar to the Within-Tunnel Constraints, for the across-tunnel direction. We proposed two formulations for single tunnel, alternating constraints and extended formulation. In theory we could apply either of these in the within and across-tunnel directions, for a total of four combinations. The combination of extended formulation for both directions was used since the extended formulation performed so well in the single tunnel setting. The Basic model formulation is shown below. Note that the draw point index is now two dimensional such that ·, dpt is within tunnel and tn, · is across-tunnel. In the draw point numbering example shown in Figure 3.12, the six horizontal tunnels are numbered from bottom to top, and the draw points from left to right. The across-tunnel dimension runs from bottom to top and includes the (tn, dpt), (tn, dpt + 1) pair from each tunnel, for each odd dpt. Data : T, m, λ nTNLS the number of tunnels in the mine 78 Chapter 3. Sequence Optimization in Block Cave Mining Figure 3.12: Example of Tunnel Numbering Fdpttn the number of the first draw point in tunnel tn, 1 ≤ tn ≤ nT N LS Ldpttn the number of the last draw point in tunnel tn, 1 ≤ tn ≤ nT N LS ptn,dpt the expected duration or periods draw point dpt in tunnel tn will remain open, 1 ≤ tn ≤ nT N LS, dpt = F dpttn . . . Ldpttn wtn,dpt the expected total NPV, if started immediately, of draw point dpt in tunnel tn, 1 ≤ tn ≤ nT N LS, dpt = F dpttn . . . Ldpttn Decision Variables : Stn,dpt,t binary variable that is 1 in the period the draw point is opened 79 Chapter 3. Sequence Optimization in Block Cave Mining Ftn,dpt,t binary variable that is 1 if the draw point is the first along the tunnel to be open in the period Gtn,dpt,t binary variable that is 1 if the draw point is the first across the tunnel to be open in the period Total Number of (binary) Decision Variables: nDP T S × T × 3 Objective : max V alue = Ldpt Xtn nTX N LS wtn,dpt X Stn,dpt,t e−λ(t−1) t tn=1 dpt=F dpttn Constraints : Start Once T X Stn,dpt,u ≤ 1 ∀ tn, ∀ dpt (3.11) u=1 Global Capacity nTX N LS Ldpt Xtn t X Stn,dpt,u ≤ m ∀ t (3.12) tn=1 dpt=F dpttn u=t−ptn,dpt Define First Along Tunnel A draw point is the first along the tunnel to be open if it is open but the previous draw point is not Ftn,dpt+1,t ≥ t X Stn,dpt+1,u − u=1 t X Stn,dpt,u ∀ t, ∀ tn, ∀ dpt (3.13) u=1 First Along Tunnel Must be Open A draw point must be open to be the first along the tunnel t X Stn,dpt,u − Ftn,dpt,t ≥ 0 ∀ t, ∀ tn, ∀ dpt (3.14) u=1 Only One First Along Tunnel Each tunnel can only have one first to be open in each period Ldpt Xtn Ftn,dpt,t ≤ 1 ∀ t, ∀ tn (3.15) dpt=F dpttn 80 Chapter 3. Sequence Optimization in Block Cave Mining Define First Across Tunnel A draw point is the first across the tunnel to be open if it is open but the previous draw point is not. The next draw point is either in the same tunnel or the adjacent one Gtn,dpt,t ≥ t X Stn,dpt,u − u=1 t X Stn,dpt+1,u ∀ t, ∀ tn, ∀ dpt odd u=1 (3.16) Gtn,dpt,t ≥ t X u=1 Stn,dpt,u − t X Stn−1,dpt−1,u ∀ t, ∀ tn, ∀ dpt even u=1 (3.17) First Across Tunnel Must be Open A draw point must be open to be the first across the tunnel t X Stn,dpt,u − Gtn,dpt,t ≥ 0 ∀ t, ∀ tn, ∀ dpt (3.18) u=1 Only One First Across Tunnel nTX N LS (Gtn,dpt,t + Gtn,dpt+1,t ) ≤ 1 ∀ t, ∀ dpt odd (3.19) tn=1 Total Number of Constraints : nDP T S+T +nDP T S×T ×4+nT N LS×T ×2 Results We had success with data set P1 (see Section 3.4.2); the formulation produced a single, contiguous cave that grew over time. The progression is shown in Figure 3.13 Unfortunately, the results for data set P2 were not so nice, see Figure 3.14. In the first period the opening pattern has two distinct caves, violating our goal of a single, contiguous cave. As this result does not violate the constraints as written in the formulation, it points out a failing of the formulation. Conceptually, in order to guarantee a single, contiguous cave these constraints would have to be applied across every possible axis through the mining footprint. Our formulation only addresses two axes, along and across the tunnels and hence does 81 Chapter 3. Sequence Optimization in Block Cave Mining First Period Opening Pattern First 6 Periods Opening Pattern First 12 Periods Opening Pattern First 17 Periods Opening Pattern Figure 3.13: Opening Patterns from Data Set P1 using Extended Formulation Within-Tunnels and Across-Tunnels not guarantee a single cave. It apparently would take a much bigger, more complicated formulation to try to truly enforce contiguity. Even though this model does not guarantee a single cave, it could prove useful to the mine planner. In Section 3.5 the performance of this model is compared to others we developed for this problem. If a single cave is produced, it is a feasible solution. If a single cave is not produced, the solution highlights high value draw point groups in the mining footprint. The result can also form an upper bound against which solutions from other models can be compared. The appearance of multiple caves in data set P2, but not in data set P1 led us to rely more heavily on data set P2 than data set P1 for the development of the other models. 82 Chapter 3. Sequence Optimization in Block Cave Mining First Period Opening Pattern First 20 Periods Opening Pattern Figure 3.14: Opening Patterns from Data Set P2 using Extended Formulation Within-Tunnels and Across-Tunnels Contribution The Basic formulation presented in this section is a new idea I developed for this thesis work. The formulation and implementation is new work on the problem of Sequence Optimization and is part of the contribution of this thesis to the research literature. 83 Chapter 3. Sequence Optimization in Block Cave Mining 3.4.4 Malkin and Wolsey’s 2D Integral Formulation Motivated by the possibility of obtaining a disconnected solution from the two-dimensional solution proposed in Malkin (2006), Malkin and Wolsey proposed an integral two-dimensional formulation (Section 3.1). This formulation chooses a starting or centre point with coordinates (r, c), and ensures that the row and column of this centre point form a dominant axis of the opened draw points in each period. They write “The row with the largest interval must be r and the column with the largest interval must be c” and “The region somewhat resembles a diamond” which is good in practice. Formulation The formulation suggested by Malkin and Wolsey, translated to our standard variables, is as follows Data : r, c the coordinates of the centre point, where the diamond axes intersect Decision Variables : Stn,dpt,1 Constraints : The column-wise neighbour between the draw point and column c must be open to open draw point tn, dpt: Stn,dpt,1 ≤ Stn,dpt+1,1 ∀ tn, ∀ dpt < c (3.20) Stn,dpt,1 ≤ Stn,dpt−1,1 ∀ tn, ∀ dpt > c (3.21) Similarly, the row-wise neighbour between the draw point and row/tunnel r must be open to open draw point tn, dpt: Stn,dpt,1 ≤ Stn+1,dpt,1 ∀ tn < r, ∀ dpt (3.22) Stn,dpt,1 ≤ Stn−1,dpt,1 ∀ tn > r, ∀ dpt (3.23) 84 Chapter 3. Sequence Optimization in Block Cave Mining Adapting to Mining Layout The formulation proposed by Malkin has two apparent drawbacks. The first is that the centre point must be specified. This can be overcome by running the formulation using each draw point as the centre point (nDP T times) and evaluating the objective function for each centre point. This may be feasible for a smaller mine but may be too time consuming for the larger mines. Alternatively, a heuristic could be developed to select good centre points. The second drawback is that since the axes are set by the centre point, the axes for the diamond are fixed for the duration of the mine. In reality, a better schedule might result from allowing the axes to vary over time as is possible under current methods. This formulation requires two approximately perpendicular axes which can be realized in the tunnels and the across direction discussed in the previous section. We use the convention that in the mining data, the “rows” are the tunnels and the “columns” are the across-tunnel pairs of draw points. Using the naming convention (tunnel, across), the draw points zigzag down the tunnels as shown earlier in Figure 3.12. If we choose the centre point 4, 21 we get the axes for the diamond as shown by the dark line in Figure 3.15. It might seem that the results from centre point 4, 21 and centre point 4, 22 would be the same, but this is not necessarily true. For example, if draw point 4, 22 has a low value, a cave generated from centre point 4, 21 may not include 4, 22 and draw points in adjacent tunnels. However, draw points on these tunnels may be included if low valued 4, 22 were chosen as the centre point, resulting in a different cave. Malkin Model Formulation The Malkin formulation is as follows Data : nT N LS, F dpttn , Ldpttn , T, m, ptn,dpt , wtn,dpt , λ TnC,DptC the coordinates of the centre point, where the diamond axes intersect Decision Variables : Stn,dpt,t Total Number of (binary) Decision Variables: nDP T S × T 85 Chapter 3. Sequence Optimization in Block Cave Mining Figure 3.15: Example of Centre Point 4,21 Objective : max V alue = nTX N LS Ldpt Xtn tn=1 dpt=F dpttn wtn,dpt X Stn,dpt,t e−λ(t−1) t Constraints : Column Connected The column-wise neighbour between the draw point and column c must be open in order to open draw point i, j: t X u=1 Stn,dpt,u ≤ t X Stn,dpt+1,u ∀ tn, ∀ dpt < DptC, ∀ t (3.24) u=1 86 Chapter 3. Sequence Optimization in Block Cave Mining t X Stn,dpt,u ≤ u=1 t X Stn,dpt−1,u ∀ tn, ∀ dpt > DptC, ∀ t (3.25) u=1 Row Connected The row-wise neighbour between the draw point and row r must be open to open draw point i, j: t X t X Stn,dpt,u ≤ u=1 Stn+1,dpt+1,u ∀ tn < T nC, ∀ dpt odd, ∀ t u=1 (3.26) t X Stn,dpt,u ≤ u=1 t X Stn,dpt−1,u ∀ tn < T nC, ∀ dpt even, ∀ t u=1 (3.27) t X Stn,dpt,u ≤ u=1 t X t X Stn,dpt+1,u ∀ tn > T nC, ∀ dpt odd, ∀ t (3.28) u=1 Stn,dpt,u ≤ u=1 t X Stn−1,dpt−1,u ∀ tn > T nC, ∀ dpt even, ∀ t u=1 (3.29) Start Once (3.11) Global Capacity (3.12) Total Number of Constraints: nDP T S + T + nDP T S × T × 2 Figure 3.16 illustrates the Column Connected Constraints in the tunnel setting. The arrows point to the draw point that must be open before the indicated draw point can open. The dark line shows the across-tunnel axis, and it is evident that all the arrows point towards this axis. The draw points directly on the axis, in this example those numbered 21, do not have constraints. Figure 3.17 shows the Row Connected constraint with the second axis at tunnel 4, shown in black. Again, the arrows point to the draw point that must be open before the indicated draw point can open. It is evident that all the arrows point towards the within-tunnel axis, tunnel 4. The draw points directly on the axis do not have constraints. The combined picture is presented in Figure 3.18 87 Chapter 3. Sequence Optimization in Block Cave Mining Figure 3.16: Example of Centre Point 4,21 and Direction of the Column Connected Constraint((3.24)-(3.25)). Each arrow represents a constraint and points to the draw point that must be opened before the draw point at the tail of the arrow can open. 88 Chapter 3. Sequence Optimization in Block Cave Mining Figure 3.17: Example of Centre Point 4,21 and Direction of the Row Connected Constraint((3.26)-(3.29)). The arrows show the draw point that must be open before the draw point at the tail of the arrow can open. 89 Chapter 3. Sequence Optimization in Block Cave Mining Figure 3.18: Example of Centre Point 4,21 and Arrows for All Constraints 90 Chapter 3. Sequence Optimization in Block Cave Mining Theoretical Properties The basic formulation of the one-period model as proposed by Malkin and Wolsey, (3.20) — (3.23), gives an integral solution. A linear programming problem of the form max{cx : Ax ≤ b, x ∈ Rn+ } will have an optimal solution that is integral if the matrix A is totally unimodular Wolsey (1998). It is straightforward to show that A is unimodular and we conclude that the polyhedron is integral. The multiperiod formulation presented is more complicated for several reasons. First, the simple neighbour constraints are extended to the multiperiod setting. Secondly, the Start Once and Global Capacity constraints are added. Nominally, the multiperiod constraints say that a given draw point can not open in period t unless the neighbour closer to the axis started in period t or prior. These could be written in the form Stn,dpt,t ≤ t X Stn,dpt+1,u ∀ tn, ∀ dpt < DptC. (3.30) u=1 Instead we have chosen to use an aggregated form of the constraints (3.24) — (3.29) which we will show is a tighter formulation. In the aggregated form, the constraints say a given draw point can only have started in period t or earlier if the neighbour closer to the centre point started in period t or prior, for all periods 1 ≤ t ≤ T . To demonstrate the aggregated constraints are stronger than the disaggregated constraints we look at the following example for a single draw point. Period one In period one both constraints are the same Stn,dpt,1 ≤ Stn,dpt+1,1 Period two In period two the disaggregated constraint is Stn,dpt,2 ≤ Stn,dpt+1,1 + Stn,dpt+1,2 and the aggregated constraint is Stn,dpt,1 + Stn,dpt,2 ≤ Stn,dpt+1,1 + Stn,dpt+1,2 This constraint is stronger than the disaggregated one as there are two non-negative variables on the left-hand side instead of the single one 91 Chapter 3. Sequence Optimization in Block Cave Mining in the disaggregated constraint. There is also no possibility of having Stn,dpt,1 = 1 and Stn,dpt+1,2 = 1 since the period one constraint would require that if Stn,dpt,1 = 1 then Stn,dpt+1,1 = 1 and then the Start Once constraint would not allow Stn,dpt+1,2 = 1 The following example shows a fractional solution that satisfies the disaggregated constraints (3.30) but violates the aggregated constraints (3.24): 1 1 Stn,dpt,1 = , Stn,dpt,2 = 2 2 3 Stn,dpt+1,1 = , Stn,dpt+1,2 = 0. 4 The disaggregated constraints are satisfied: Stn,dpt,1 = 1 3 ≤ Stn,dpt+1,1 = 2 4 1 3 ≤ Stn,dpt+1,1 + Stn,dpt+1,2 = , 2 4 but the aggregated constraints are violated Stn,dpt,2 = Stn,dpt,1 = 1 3 ≤ Stn,dpt+1,1 = 2 4 3 Stn,dpt,1 + Stn,dpt,2 = 1 6≤ Stn,dpt+1,1 + Stn,dpt+1,2 = . 4 The multiperiod setting also requires the addition of the Start Once constraint. This constraint is required to ensure that the valuable draw points are not opened multiple times in order to increase the objective value. An alternative may be to only restrict the centre point, draw point T nC, DptC to open once. T X ST nC,DptC,u ≤ 1 t=1 Combining thisP constraint with the aggregated neighbour constraints, should limit the sum Tt=1 Stn,dpt,u to one for all draw points thus ensuring each opens at most once. This constraint was not tested computationally. Finally, the Global Capacity constraint is needed in its multiperiod form to reflect the non-uniformity of draw point durations. If all draw points had uniform duration then this constraint and the overall formulation could be 92 Chapter 3. Sequence Optimization in Block Cave Mining simplified. Withj uniform duration p one is only interested in periods t = kp k T for k = 1, 2, . . . p and the Global Capacity constraint becomes nTX N LS Ldpt Xtn Stn,dpt,t ≤ m ∀ t tn=1 dpt=F dpttn Initial inspection suggests that the multiperiod formulation may be integral with uniform duration of all draw points, but that with varying durations the Global Capacity constraint becomes an obstacle to integrality. Results This formulation provides solutions with one, contiguous open cave if the mine footprint is convex. The shape may not always be a diamond, but it is based on two axes. Sample results from data set P2 are shown below in Figure 3.19. The tunnels are the zig-zag lines running approximately parallel to the horizontal axis and the axes for the formulation are shown in black and intersect at the chosen centre point. However, when the footprint is not convex, multiple caves can form. This can be seen in Figure 3.20 which shows results on data set P4. The draw points on the horizontal axis of the diamond are only constrained to have a neighbour closer to the vertical axis of the diamond. In the case of a gap in the footprint, the neighbour to the edge draw point is on the other side of the gap, and it opened in period 1. This allows draw points along the horizontal axis of the diamond to open and form an additional cave. This can be fixed by bending the axes around the footprint. We have chosen to bend the axes around the opening in the mine footprint, then continue in a straight line tangent to the opening. An example of this is in Figure 3.21. This method finds a solution, and finds it relatively quickly for some centre points in data set P2. Results from runs from data set P2 and data set P4 are shown in Table 3.3 for comparison. The comment Original Axes refers to the setting of the axes from the center point even if they cross a empty section of the mining footprint. The comment Bent Axes means that the axes are bent around any empty sections. There are no such comments for data set P2 as there are no empty sections in this data set. While some runs from data set P2 complete in less than a minute, most runs from data set P4 still take over an hour. This table also indicates that there are some variations in objective values from different centre points. Most of the cases 93 Chapter 3. Sequence Optimization in Block Cave Mining Figure 3.19: Data Set P2, Period 1 Opening where the axes are bent took longer to solve than those with the original axes, but not in all cases. Since the solutions without bent axes resulted in multiple caves, the extra time is the price for a feasible solution. Similarly, in most but not all cases the bending of the axes increased the objective value. We are still faced with the problem of finding a good centre point. Since data set P2 has only 236 draw points, the model was repeatedly run using each draw point as the centre point of the formulation. The results of these runs show the variation in objective value as well as differences in running time over the footprint. Table 3.4 and Figure 3.22 below show some of the results from these runs. 94 Chapter 3. Sequence Optimization in Block Cave Mining Data Set P2 P2 P2 P2 P2 Run Time (sec) 89.8 55.2 4,217.9 11.6 25,763.7 Centre Tunnel 5 6 1 5 10 Centre Dpt 12 12 20 10 21 Total NPV 452,746 451,207 412,537 451,894 424,758 comments min Total NPV fastest run time longest run time P4 P4 P4 P4 P4 P4 P4 P4 P4 P4 3,178.6 5,280 4,176.4 3,105.6 3,076.6 4,028.3 4,488 6,184.9 2,978.6 3,605 13 13 14 14 14 14 28 31 31 31 90 90 65 65 95 95 135 57 103 103 606,162 609,400 608,322 608,648 604,280 607,859 606,454 607,713 609,901 609,416 Original Axes Bent axes Original Axes Bent axes Original Axes Bent Axes Original Axes Original Axes Original Axes Bent Axes max Total NPV Table 3.3: Table of Running Times of Malkin Model from Data Sets P2 and P4 95 Chapter 3. Sequence Optimization in Block Cave Mining Data Set P2 P2 P2 P2 P2 P2 P2 P2 P2 P2 P2 P2 P2 P2 P2 P2 P2 P2 P2 P2 P2 P2 P2 P2 P2 P2 P2 P2 P2 P2 P2 .. . Run Time (sec) 544.234 387.891 196.031 17.188 120.422 517.094 879.984 173.406 178.844 168.516 371.829 423.516 4714.88 4217.95 1052.12 447.437 117.297 91.828 95.094 96.234 69.625 75.36 266.407 128.313 664.437 178.843 201.047 241.547 1398.74 2149.45 557.906 .. . Centre Tunnel 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 .. . Centre Dpt 7 8 9 10 11 12 13 14 15 16 17 18 19 20 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 5 .. . Total NPV 435,414 436,801 443,198 443,198 443,232 444,724 443,340 442,610 437,418 436,455 430,991 429,554 413,992 412,537 431,796 432,407 439,908 441,360 445,872 446,274 446,743 447,803 447,159 446,376 442,037 441,616 434,607 434,960 422,581 421,452 436,816 .. . Table 3.4: Partial Table of Running Times for Data Set P2 96 Chapter 3. Sequence Optimization in Block Cave Mining Figure 3.20: Data Set P4, Centre Point (14,65) 97 Chapter 3. Sequence Optimization in Block Cave Mining Figure 3.21: Horizontal Axis Bent Around Hole in Footprint, Centre Point (14,65) 98 Chapter 3. Sequence Optimization in Block Cave Mining Figure 3.22: Histogram of Total NPV There is a range of about 9% from the lowest NPV of 412,537 to the highest NPV of 452,746 for data set P2, so it is worth some effort to find the centre point yielding the highest value. The distribution is shown in Figure 3.22. Note, these runs were done with the CPLEX setting of MIPGAP 0.01. This means that computation is stopped when a solution within 1% of optimality is reached. For this range of objective values, 1% is approximately 4,400 which explains some of the variation in the results. A surface plot of objective values over the footprint of the mine is shown in Figure 3.23 for data set P2. For this data set, highest objective values are for centre points in the centre of the mine, dropping off quite steeply at the mine edges. While there is some variation at the centre, there is no single, small cluster of high objective value points. Figure 3.24 is a preliminary attempt to heuristically locate the optimal centre point. Each of the three graphs is a contour plots over the footprint of the mine for data set P2. The draw points are the intersection of the grid lines. The tunnels run horizontally across the plot and the across numbering runs vertically. The plot on the right is for the TOTAL NPV given the draw point is the centre point for the Malkin model. The dark shape in the 99 Chapter 3. Sequence Optimization in Block Cave Mining Figure 3.23: Surface Plot of Total NPV for Data Set P2 centre indicates the highest value draw points. The two plots on the left are representations of the value of the individual draw points. The top one is the ratio wpii , the total value divided by the duration. The lower one is the ratio wi modified to represent the discount factor and is exp(−λp . This figure does i) not support the conjecture that individual draw point value predict TOTAL NPV if used as a centre point. There are three clusters of draw points with the highest wpii ratio. The top group does not at all overlap with the high TOTAL NPV draw points, and there are perhaps one or two from the other wi two groups that overlap. There are also three clusters of high exp(−λp ratio i) draw points, but there is only a single draw point (tunnel 4, draw point 10) that is also in the cluster of high TOTAL NPV draw points. 100 Chapter 3. Sequence Optimization in Block Cave Mining Figure 3.24: Projections of the Surface Plots for Draw Point Values of wi exp(−λpi ) and Total NPV for Data Set P2 wi pi , 101 Chapter 3. Sequence Optimization in Block Cave Mining Variation: Wandering Axes formulation One potential drawback to the formulation proposed by Malkin and Wolsey is that the axes for the diamonds are fixed over the entire duration of the mine. It is not clear whether this is an necessary restriction on the mine growth or not, and it may overly constrain the problem. The following variation of the above formulation allows the axes to deviate from the tunnel and draw point chosen as the centre point. The idea is that if the axes are allowed to deviate, or wander, then the formulation may lead to better objective values; we have named it the Wandering Axes formulation. The variation is executed by picking a centre point, then deciding where the axes will run. When this has been done, the model runs as above with the constraints requiring neighbours between the draw point and the axes at the location of the draw point to be open before the draw point can start. The centre point is heuristically chosen as the draw point with the highest average value of a cluster of neighbours, and the progress of the axes are also influenced by the average value of neighbours. The neighbours are defined as those within a diamond shape containing approximately m draw points around the centre point. For simplicity we looked at diamond shapes of size n(n−1) n−1 + n−1 2 2 when n is an odd integer. The diamond extends 2 draw points along the tunnel, in either direction from the centre point. In the adjacent tunnel there are n − 2 draw points in the diamond. Each adjacent tunnel is included with the number of draw points decreasing by two until n−1 2 tunnels away there is a single draw point in the diamond. Since not every centre draw point will have a full diamond of neighbours ( for example the draw points at the edge of the mining footprint ) the average value over the number of existing neighbours is used. The following procedure was used. 1. For each draw point the average neighbour value was calculated. The total value of the draw point wi was used. The centre point was chosen as the draw point with the highest neighbour value. For data set P2 a diamond of size 25 was used and the chosen centre point was tunnel 6, draw point 12. 2. The axes started at the chosen centre point T nC, DptC and stayed with the centre tunnel and draw point in the adjacent tunnels and draw points for one step. From here the axis following the centre draw point alternates between choosing the direction of highest neighbour value in one tunnel, holding that same path in the next tunnel, then again deciding the highest value route in the next tunnel. That is, in 102 Chapter 3. Sequence Optimization in Block Cave Mining the next adjacent tunnels T nC ± 2, the neighbour values of three draw points (T nC + 2, DptC − 1), (T nC + 2, DptC), (T nC + 2, DptC + 1) are compared, and the axis moves to the draw point with the highest neighbour value. The axis stays the same in the adjacent tunnel, but in the next, a similar choice is made. Similarly, for the along tunnel axis a decision is made at DptC + 2 between (T nC − 1, DptC + 2), (T nC, DptC + 2), (T nC + 1, DptC + 2) for the draw point with the highest neighbour value, and also at DptC − 2 between (T nC − 1, DptC − 2), (T nC, DptC − 2), (T nC + 1, DptC − 2). An example is shown in Figure 3.25 Figure 3.25: Table of Neighbour Values for Data Set P2. Rows are tunnels and columns are draw point numbers. The boxes show the three values considered in each decision and the bold values indicated the selected axes. In the data set P2 this criterion chose tunnel 6, draw point 12 as the centre point. The results from this run of the model are shown in Table 3.5. Runs from the basic Malkin model with centre points 5, 12 and 6, 12 are included for comparison. Data Set P2 P2 P2 Run Time (sec) 54.6 89.782 55.156 Centre Tunnel 6 5 6 Centre Dpt 12 12 12 Total NPV 452,530 452,746 451,207 Cplex Optimality gap 0.10% 0.02% 0.36% comments Wandering Axes basic Malkin basic Malkin Table 3.5: Comparison of Results of Wandering Axes Variation The results show that the Wandering Axes variation did not give a higher objective value than the basic Malkin model, in that the values of 452,746 for the basic model and 452,530 for the Wandering Axes variation are within the 103 Chapter 3. Sequence Optimization in Block Cave Mining optimality gap of CPLEX. The Wandering Axes variation was successful as a single run heuristic of approximating in a single run the maximum objective value obtained by running the basic formulation over all draw points. Figure 3.26: Results from the Wandering Axes Model Variation on Data Set P2 Figure 3.26 shows the resulting plan for the mine. The black lines show the axes. While the along tunnel axis does “wander” there is little movement in the across-tunnel axis. Figure 3.27 shows the resulting plan for the mine from the basic Malkin model centred at tunnel 6 draw point 12. The black lines show the axes and Figure 3.28 compares wandering axis and basic Malkin model results side by side. 104 Chapter 3. Sequence Optimization in Block Cave Mining Figure 3.27: Results from the Malkin Model Centered at Tunnel 6, Draw Point 12 (a) (b) (Wandering Axes) (Basic Malkin) Figure 3.28: Same Figures Side by Side 105 Chapter 3. Sequence Optimization in Block Cave Mining This model runs quickly and does produce a single, contiguous cave. The best value for data set P2 is 452,746. The Wandering Axes heuristic did not produce a higher objective value, but did closely approach the best value over the data set. Results from this data set suggest that the extra work required to determine the changes in the axes was not justified by improved objective value. The formulation does work so could be evaluated on other data sets if there was reason to believe an improvement is achievable. The Malkin model guarantees a single cave in each period which is an improvement over the Basic model. The main cost of this model is the determination of the fixed centre point. For smaller data sets it may prove feasible to run the model for all draw points. In Section 3.5 the performance of this model is compared to others we developed for this problem. Contribution The single period formulation of the Malkin model was suggested by Malkin and Wolsey after discussions with Maurice Queyranne. They have not published this work. Nor, to my knowledge, have they done further work on it or Sequence Optimization. The extension of the formulation to multiple periods and the implementation on data sets is new work I did for this thesis and is part of the contribution of this thesis to the research literature. 106 Chapter 3. Sequence Optimization in Block Cave Mining 3.4.5 Formulations Based on 4 Vertices A third approach to formulating this problem is rooted in the operational preference for a diamond-shaped group of draw points in each period (Deiring (2006)). The mining engineers favour a cave with a diamond shape as it has the desired area to perimeter ratio to have good rock behaviour. Ideally a square cave is favoured over a long and thin one. Diering (2006) defines a diamond as “On a given front [of the cave], the line of draw points is relatively straight”. We define a diamond as a shape with 4 vertices such that the diagonals intersect at right angles. The formulations that follow are based on this definition. 4 Vertices In the 4 Vertices model the general idea is that a diamond is defined by its 4 vertices which induce its 4 edges and interior points. The points in a diamond can be described by the number of open neighbours. In the mine layout a typical draw point has two within-tunnel neighbours and one across-tunnel neighbour. Draw points at the edges of the footprint may only have a subset of these three neighbours. Clearly a point in the interior of the diamond has all three neighbours open. Points on an edge have two neighbours open and a vertex has either one or two open neighbours. Diamond Example 1 Diamond Example 2 Figure 3.29: Example of Diamonds with Vertices Having Only One Open Neighbour 107 Chapter 3. Sequence Optimization in Block Cave Mining Figure 3.29 shows examples of diamonds in which the vertices have only one neighbour. The zig zag lines represent the tunnels and the draw points occur at the bends in the zig-zag. The vertical lines connect the draw points to their across-tunnel neighbours. The vertices are marked with asterisks, draw points with two open neighbours are marked with open circles, and those with three open neighbours are marked with open squares. With this picture in mind, an attempt was made to write a formulation to open draw points if they either have at least two open neighbours or are one of the four vertices. The size of the diamond is constrained by the Global Capacity constraint, and each draw point can only be opened once. There were difficulties in writing the formulation for the constraint to achieve that a draw point requires two neighbours open, or be one of the 4 vertices in order to open. The attempts quickly identified two problems with this model concept. 1. The vertex definition of having only one open neighbour does not hold for all possible diamonds. 2. Interior points are not forced to have all three neighbours open and holes can form within the diamond The first problem is illustrated in Figure 3.30. The draw point in the top vertex position has 2 open neighbours and so can open without being designated as a vertex. As a result, another draw point, not in this diamond, can be assigned the vertex and initiate another diamond, not connected to the one shown. The second problem arises from the constraint that a draw point must have at least two open neighbours to open. Ideally there would be one constraint for edges, allowing only two open neighbours, and one for interior draw points requiring all three. We were not able to write a set of constraints that differentiated between edge and interior points. Instead, a more restrictive definition of the diamond was sought. 108 Chapter 3. Sequence Optimization in Block Cave Mining Figure 3.30: Example of Diamond with a (Top) Vertex Having Two Open Neighbours 109 Chapter 3. Sequence Optimization in Block Cave Mining 2 Vertices and Cones The 2 Vertices and Cones formulation is an adaptation of the 4 Vertices formulation to specifically address the target of one contiguous cave. In this formulation the diamond is identified by a “top” vertex (Vertex1), a “bottom” vertex (Vertex4), and the intersection of a upward-pointing cone defined by Vertex1, and an downward-pointing cone defined by Vertex4. If a draw point is Vertex1 in period t then the draw points in the cone above it are forced open to form BottomCone for period t. Similarly if a draw point is Vertex4 in period t then the draw points in the cone below it are forced open to form TopCone for period t. The draw points in the intersection of the two cones must either have opened in an earlier period i < t or in that period t. Figure 3.31: Three Examples of Pairs of Cones and Their Intersections. Only when a vertex from each cone are in both cones (middle case) does a diamond result The intersections of the two cones will only give the desired diamond shape if both vertices are in both cones, see Figure 3.31. A constraint is desired that will ensure that both vertices are in both cones. The definition of the cones defines the possible shapes that can occur. More specifically, the definition of the cones fixes the interior angles at each vertex. Two attempts were made at this model. In the first, unsuccessful, attempt, constraints were written to define which draw points could be included in a cone, but allow the model to determine the actual shape of the cones during execution. The cone definition began with the corresponding 110 Chapter 3. Sequence Optimization in Block Cave Mining vertex, then neighbours adjacent to the vertex could be included, then their neighbours etc. Despite all attempts, the formulations either generated multiple caves in a single period or over multiple periods, and proved to have many of the problems of the 4 Vertices model; the results are not presented here. The second attempt was based on a much more rigid definition of the cones. In this 2Cone model there is a defined cone that is forced open for each vertex. Based on discussions with Gemcom (Deiring (2006)), we proceeded with cones that grow by two draw points per incremental tunnel. If the intertunnel spacing equals the inter draw point spacing this will result in a base angle of approximately 45 degrees. 2Cone Formulation If the vertices are limited to the footprint of the mine, then it will never be possible to include all the draw points in the defined cone and open them. For example, if the top vertex is in the top tunnel at draw point dpt, then at most its two within-tunnel neighbours dpt − 1 and dpt + 1 will be in the cone, and dpt − 2 can not open unless the vertex moves to dpt − 1 in the next period. To overcome this difficulty, the space for the vertices is extended beyond the mine footprint with extra tunnels. In general, enough tunnels should be added below the first tunnel and above the last tunnel so that the maximum tunnel length in the footprint can be in a cone if it occurred any tunnel. An extended footprint is created by adding eT N LS tunnels above and below the mine footprint. For a cave that grows at two draw points per tunnel, the eT N LS is one half the longest tunnel length. An example is shown in Figure 3.32. Data : nT N LS, T, m, ptn,dpt , wtn,dpt , λ eTNLS the number of tunnels added above and below the mine footprint to create the extended footprint. For a cave that grows by two draw points per tunnel, the value is half of largest acrosstunnel value for all tunnels, eT N LS = d0.5 ∗ maxtn (Ldpttn )e. Decision Variables : Stn,dpt,t binary variable that is 1 in the period the draw point is in both cones in a given period. Note the definition of this variable is altered slightly for this model only. 111 Chapter 3. Sequence Optimization in Block Cave Mining Figure 3.32: An Example of the Extended Footprint for the 2Cone Model. The outline of the mine footprint is shown, and the 8 tunnels (horizontal lines) shown in solid lines. Since the longest tunnel is 20 draw points long, eT N LS = 10 and 10 additional tunnels above and below the mine footprint are shown in dashed lines. Two example cones are shown, one originating at a vertex in the extended footprint. Vertex1tn,dpt,t binary variable that is 1 if the draw point is the bottom (lowest tunnel) vertex in a given period. This variable is defined over the extended footprint 1 ≤ tn ≤ nT N LS + eT N LS, 1 ≤ dpt ≤ 2 ∗ eT N LS Vertex4tn,dpt,t binary variable that is 1 if the draw point is the top (highest tunnel) vertex in a given period. This variable is defined over the extended footprint 1 + eT N LS ≤ tn ≤ nT N LS + 2 ∗ eT N LS, 1 ≤ dpt ≤ 2 ∗ eT N LS BottomConetn,dpt,t binary variable that is 1 if the draw point is the the cone projecting up from Vertex1 in a given period TopConetn,dpt,t binary variable that is 1 if the draw point is in the the cone projecting down from Vertex4 in a given period 112 Chapter 3. Sequence Optimization in Block Cave Mining Total Number of (binary) Decision Variables: nDP T S × T × 5 Objective : max V alue = Ldpt Xtn nTX N LS wtn,dpt X Stn,dpt,t e−λ(t−1) t tn=1 dpt=F dpttn Constraints : Start Once (3.11) Global Capacity (3.12) One Vertex4 : in each period nT N LS+2∗eT N LS X N LS 2∗eT X tn=1+eT N LS V ertex4tn,dpt,t ≤ 1 ∀ t (3.31) dpt=1 One Vertex1 : in each period nT N LS+eT N LS X N LS 2∗eT X tn=1 V ertex1tn,dpt,t ≤ 1 ∀ t (3.32) dpt=1 Define Top Cone : Starting at Vertex4, the cone grows by two draw points per tunnel. T opConetn,dpt,t = nT N LS+eT X N LS dpt+(tun−tn) X tun=tn V ertex4tun+eT N LS,d,t d=dpt−(tun−tn) ∀ tn, ∀ dpt, ∀ t (3.33) Define Bottom Cone BottomConetn,dpt,t = tn X dpt+(tn−tun) X V ertex1tun+eT N LS,d,t tun=1−eT N LS d=dpt−(tn−tun) ∀ tn, ∀ dpt, ∀ t (3.34) 113 Chapter 3. Sequence Optimization in Block Cave Mining Open Intersection of Cones draw points in both cones must be in the pattern In Both t−1 X Stn,dpt,u ≥ T opConetn,dpt,t +BottomConetn,dpt,t −1 ∀ tn, ∀ dpt, ∀ t u=1 (3.35) In Top Stn,dpt,t ≤ T opConetn,dpt,t ∀ tn, ∀ dpt, ∀ t (3.36) In Bottom Stn,dpt,t ≤ BottomConetn,dpt,t ∀ tn, ∀ dpt, ∀ t Previous Cave in Current Cave t−1 X Stn,dpt,u ≤ BottomConetn,dpt,t ∀ tn, ∀ dpt, ∀ t ≥ 2 (3.37) (3.38) u=1 t−1 X Stn,dpt,u ≤ T opConetn,dpt,t ∀ tn, ∀ dpt, ∀ t ≥ 2 (3.39) u=1 Total Number of Constraints : nDP T S+T ×3+nDP T S×T ×7−nDP T S×2 ———————————————————————— Results The 2Cone model was successfully run on data set P2. The model ran in 22 seconds and returned a single cave in each period with a total NPV of 356,187. This model significantly faster than the Malkin model, which took 89 seconds on the maximum value run alone. The value of the 2Cone result was approx 79% of value of the maximum value Malkin run. The initial cave was centred in the middle of the footprint, on tunnels 6 and 7, close to the centre point for the maximum value Malkin run. In Section 3.5 the performance of this model is compared to others we developed for this problem on two trial data sets. Contribution The 2Cone formulation presented in this section is a new idea I developed for this thesis work. The formulation and implementation is new work on the problem of Sequence Optimization and is part of the contribution of this thesis to the research literature. 114 Chapter 3. Sequence Optimization in Block Cave Mining 3.4.6 Column Generation Since other approaches did not seem to guarantee a single contiguous cave that could move freely over the mining footprint, a completely different approach was proposed by Maurice Queyranne. This approach uses the well known technique called column generation Wolsey (1998). Column generation breaks the problem down into a Master problem which selects a solution from a set of presented options, and a Subproblem that generates additional feasible options that will improve the objective function. The presented options are denoted as ‘columns’ and the Subproblem generates more ‘columns’ for the Master problem, hence the name column generation. The appeal of this technique for solving Sequence Optimization is that while the Master problem must be run as a linear/integer program, the Subproblem is not limited to the linear/integer programming structure. Note that these “columns” are not the same columns used to describe the mining footprint. Clearly, if we could enumerate all the possible opening patterns at the start of the procedure the Master problem would be sufficient and the problem solved in one pass. Since each opening pattern is already an acceptable, single contiguous cave, the constraints in the Master problem would guarantee the other requirements such as Start Once Constraint, Global Capacity Constraint, and any others such as the lead/lag time of adjacent draw point openings. However, there are far too many possible patterns for this to work. ColGen Model Formulation The Master problem will contain a variable that indicates which columns are chosen. While this will be a binary integer variable, the linear program relaxation of the Master problem is solved in order to generate dual values/shadow prices for each of the constraints. These values are used by the Subproblem to generate a new column that is not only feasible, but also that will improve the objective function of the Master problem. After convergence is reached, the Master problem is run as an integer program to determine the final solution. Column generation has been successfully used in large scheduling problems such as course registration at colleges and airline scheduling (Lübbecke (2005)). In these settings the constraints on feasible schedules are complex, and there are a large number of feasible schedules, which is similar to our problem. The overall objective is to determine in which period to start each draw 115 Chapter 3. Sequence Optimization in Block Cave Mining point in order to maximize total NPV subject to a global capacity constraint and shape constraints on the open mine in each period. Decision Variables : Stn,dpt,t binary variable that is 1 in the period the draw point is opened. Constraints : Start Once (3.11) Global Capacity (3.12) Open Mine shape The open draw points in each period must form a contiguous shape that is roughly a diamond. Master Problem Our motivation in applying the column generation method is to remove the Open Mine Shape constraint from the rest of the problem as it has proven difficult to formulate in the linear and integer programming setting. The Master problem will choose one mine opening pattern for each period t ≤ T to form a feasible opening sequence that maximizes NPV. We construct the Master problem as Data : nT N LS, F dpttn , Ldpttn , T, m, ptn,dpt , wtn,dpt , λ Kt the number of feasible mine opening patterns for period t, 1 ≤ t ≤ T Ltn,dpt,t,k a binary array of feasible mine opening patterns (tn, dpt) for period t, 1 ≤ tn ≤ nT N LS, dpt = F dpttn . . . Ldpttn , 1 ≤ k ≤ Kt , 1 ≤ t ≤ T . These patterns are created in the Subproblem and are included as data in the Master problem. Decision Variables :Stn,dpt,t Xt,j binary variable that is 1 if mine opening pattern L∗,∗,t,j is chosen in period t, 1 ≤ t ≤ T Objective : max V alue = nTX N LS Ldpt Xtn tn=1 dpt=F dpttn wtn,dpt X Stn,dpt,t e−λ(t−1) t 116 Chapter 3. Sequence Optimization in Block Cave Mining Constraints : Start Once (3.11) Global Capacity (3.12) Link Selected Mines Together: First period Stn,dpt,1 = allpatterns X Ltn,dpt,1,k X1,k ∀ tn, ∀ dpt (3.40) k=1 Link Selected Mines Together: All other periods Stn,dpt,t = allpatterns X Ltn,dpt,t,k Xt,k k=1 − allpatterns X Ltn,dpt,t−1,k Xt−1,k k=1 ∀ tn, ∀ dpt, ∀ t > 1 (3.41) Select At Most One Pattern Per Period allpatterns X Xt,j ≤ 1 ∀ t (3.42) j=1 Improving Objective Function When the Master problem is run as a linear program, a shadow price will be generated for each constraint. In choosing new patterns for each period to add to L (the role of the Subproblem), the Master objective function will be improved if the reduced cost of the added patterns improves the objective. In this case, the objective is max V alue = nTX N LS Ldpt Xtn tn=1 dpt=F dpttn wtn,dpt X Stn,dpt,t e−λ(t−1) t for the given discount rate λ. For each period t, we seek a new pattern such that the reduced cost is greater than zero. The reduced cost is a function of the objective coefficient of X for period t, in this case zero, and the shadow prices of the constraints on X for period t. Let us assign the following dual variables to the constraints of interest. 117 Chapter 3. Sequence Optimization in Block Cave Mining Link Selected Mines Together: First Period (3.40) dual variable πtn,dpt,1 ∀ tn, ∀ dpt Link Selected Mines Together: All Other Periods (3.41) dual variable πtn,dpt,t ∀ t > 1, ∀ tn, ∀ dpt Select At Most One Pattern Per Period (3.42) dual variable µt ∀ t The dual formulation related to the X variables becomes Xt,k : − XX Ltn,dpt,t,k πtn,dpt,t + tn dpt XT,k : − XX Ltn,dpt,t,k πtn,dpt,t+1 + µ1 ∀ t < T tn dpt XX Ltn,dpt,T,k πtn,dpt,T + µT tn dpt Simplifying, the dual formulation related to the X variables becomes Xt,k : − XX Ltn,dpt,t,k (πtn,dpt,t − πtn,dpt,t+1 ) + µt ∀ t ≤ T tn dpt XT,k : − XX Ltn,dpt,T,k πtn,dpt,T + µT tn dpt An optimal solution to the (relaxed) Master problem is attained if all possible columns have negative reduced costs, that is, there is no way to generate a new column with positive reduced cost. Thus the Subproblem is to generate an opening pattern with positive reduced costs for each period to be added to the columns X of the Master problem. In each period t, we are looking to add a new pattern Ltn,dpt,t,i such that XX − Ltn,dpt,t,i (πtn,dpt,t − πtn,dpt,t+1 ) − µt ≥ 0 tn dpt Equivalently, we are looking for an “acceptable” Ltn,dpt,t,i for each period satisfying XX Ltn,dpt,t,i (πtn,dpt,t+1 − πtn,dpt,t ) − µt ≤ 0 tn dpt 118 Chapter 3. Sequence Optimization in Block Cave Mining Subproblem The Subproblem is to generate new “acceptable” patterns that satisfy X − (πtn,dpt,t − πtn,dpt,t+1 ) − µt ≥ 0 {tn,dpt} ∈ pattern for each period, if they exist. Two Subproblem routines were developed, both based on the two intersecting cone concept. The first is an integer program formulation, and the second is an enumeration routine. It was hoped that since the Global Capacity Constraint was included in the Master problem the dual values would transfer over information on the size of the pattern. The Subproblem is run for each period t ≤ T . Data : nT N LS, F dpttn , Ldpttn , eT N LS Dualtn,dpt the dual value of each draw point for the chosen period t, Dualtn,dpt = πtn,dpt,t − πtn,dpt,t+1 , 1 ≤ tn ≤ nT N LS, F dpttn ≤ dpt ≤ Ldpttn The integer programming formulation of the Subproblem for a single period is based on constructing cones with a fixed angle at the vertex. The cone starts with the vertex then grows by two draw points per tunnel. Decision Variables : Pattn,dpt binary variable that is 1 if draw point is the pattern Vertex1tn,dpt binary variable that is 1 if the draw point is the bottom (lowest tunnel) vertex in a given period. This variable is defined over the extended footprint 1 + eT N LS ≤ tn ≤ nT N LS + 2 ∗ eT N LS, 1 ≤ dpt ≤ 2 ∗ eT N LS. Vertex4tn,dpt binary variable that is 1 if the draw point is the top (highest tunnel) vertex in a given period. This variable is defined over the extended footprint 1 + eT N LS ≤ tn ≤ nT N LS + 2 ∗ eT N LS, 1 ≤ dpt ≤ 2 ∗ eT N LS BottomConetn,dpt binary variable that is 1 if the draw point is in the cone projecting up from Vertex1 in a given period TopConetn,dpt binary variable that is 1 if the draw point is in the cone projecting down from Vertex4 in a given period 119 Chapter 3. Sequence Optimization in Block Cave Mining Objective : max V alue = Ldpt Xtn nTX N LS P attn,dpt Dualtn,dpt tn=1 dpt=F dpttn Constraints : One Vertex4 nT N LS+2∗eT N LS X N LS 2∗eT X tn=1+eT N LS V ertex4tn,dpt ≤ 1 ∀ tn, ∀ dpt (3.43) dpt=1 One Vertex1 nT N LS+eT N LS X N LS 2∗eT X tn=1 V ertex1tn,dpt ≤ 1 ∀ tn, ∀ dpt (3.44) dpt=1 Define TopCone Starting at the vertex, the cone grows by two draw points per tunnel. T opConetn,dpt = nT N LS+eT X N LS dpt+(tun−tn) X tun=tn V ertex4tun+eT N LS,d d=dpt−(tun−tn) ∀ tn, ∀ dpt (3.45) Define BottomCone BottomConetn,dpt = tn X dpt+(tn−tun) X V ertex1tun+eT N LS,d tun=1−eT N LS d=dpt−(tn−tun) ∀ tn, ∀ dpt (3.46) Open Intersection of Cones draw points in both cones must be in pattern In Both P attn,dpt ≥ T opConetn,dpt + BottomConetn,dpt − 1 ∀ tn, ∀ dpt (3.47) 120 Chapter 3. Sequence Optimization in Block Cave Mining In Top P attn,dpt ≤ T opConetn,dpt ∀ tn, ∀ dpt (3.48) In Bottom P attn,dpt ≤ BottomConetn,dpt ∀ tn, ∀ dpt (3.49) The enumeration routine of the Subproblem for a single period is also based on the two cone idea. The algorithm can be summarized as follows: For each possible pair of Top Vertex, Bottom Vertex draw points • Determine the mine opening pattern for the vertex pair • Determine the size (number of draw points), and value of the mine opening • If the current pattern has a higher value than the previous contender continue, else go to next vertex pair • If the current pattern is the correct size, the current pattern replaces the contender • Go to the next vertex pair The mine opening pattern is determined by the vertex pair and, as above, starts at the top vertex and grows by two draw points per tunnel until the mid point, then reduces by two draw points per tunnel until the bottom vertex is reached. If we denote tunnel and across coordinates of the top and bottom vertices as (Ttn , Tdpt ) and (Btn , Bdpt ), for each tunnel t, Ttn ≤ t ≤ Btn , the mine opening consists of draw points {max(Tdpt − (t − Ttn ), Bdpt − (Btn − t)), . . . , min(Tdpt + (t − Ttn ), Bdpt + (Btn − t))}. Initial Results The first run was done using the data set P2, size 236 draw points. This data set was chosen both because it is small and because earlier work has shown that this data set tends to violate the one contiguous mine constraint. A value of 20 was chosen for m in the Master problem and the model run for 10 periods. The duration of the draw points range from 2 to 8 periods with very few of duration less than 4. When the iterations converged, and the master model was finally run as an integer program, there was a huge gap between the IP and LP objective values. Inspection of the generated set of patterns showed some patterns 121 Chapter 3. Sequence Optimization in Block Cave Mining smaller than the maximum active size m, and some that were clearly too large, i.e., larger than m in period 1. Neither of these groups of patterns were truly feasible patterns. Closer inspection of the intermediate results of the Master problem showed very fractional solutions that pieced together a collection of patterns to maximize the LP relaxation of the problem. Clearly the Global Capacity constraint did not transfer enough information from the Master problem to the Subproblem and it was concluded that both a minimum and maximum size constraint were needed for the Subproblem. Determining Limits on Subproblem Mine Size For period 1 the mine should be of size m. For all other periods, the size is a function of the activity in the previous periods and it is not possible to determine exact size unless all draw points have the same duration. Lower size limits for each period can be determined by relaxing the Open Mine Shape assumption and instead assuming the draw points will open in order of increasing duration. Number the draw points di such that pdi ≤ pdi+1 . Start with the m shortest duration draw points, {d1 , . . . , dm }, in period 1, then open the next one, dm+1 in period pd1 + 1 when the first active draw point closes. Continue to open draw points as the active draw points close, keeping track of the size of the mine in each period. This generates a lower bound on the number of open draw points or size of the mine in each period. Similarly, the upper size limits for each period can be determined by opening the draw points in order of decreasing duration. Since these estimation methods do not consider the feasibility the pattern size they are conservative and, if there is a wide spread of duration values, are not expected to be tight bounds. Generating Sequences of Patterns The preliminary runs also revealed that while the algorithm generated feasible mine patterns in a given period, it did not guarantee the generation of a feasible sequence of mine patterns over all periods. The linear relaxation of the model could use fractional solutions to incorporate the newly generated patterns and increase the objective value, but the final run as an integer program often reverted back to the initial set of columns to obtain a feasible solution. We considered two solutions. In the first all mine patterns generated by the Subproblem can be selected by the Master problem for any period. This increases the number of feasible solutions as a mine pattern of size m could be 122 Chapter 3. Sequence Optimization in Block Cave Mining selected in all periods. In the second, a set of mine patterns is generated by the Subproblem for each period. This solution was strongly recommended by members of the thesis committee, so it was chosen. Since it is not possible to input the pattern of the previous period pattern to the Subproblem, there is no guarantee that the model will return a set of patterns that can be linked from period to period; there may be no feasible solution to the Master problem. To accommodate this, the definition of a feasible solution to the Subproblem is extended to include mine patterns in all other periods that together form a feasible sequence with the generated pattern. Once the Subproblem generates a mine pattern for period t, an algorithm generates a feasible pattern, or precursor for the period t−1 that is subset of the pattern t. The algorithm is repeated for each period t − 1 . . . 2 generating a sequence of feasible, linked patterns that could lead to the opening of the pattern generated by the Subproblem. If a feasible sequence for periods 1 . . . t − 1 is generated, then another algorithm generates a feasible, linked sequence of patterns for periods t + 1 . . . T and all the mine patterns are added to the Master problem. If no feasible sequence is generated for period 1 . . . t − 1, then the pattern generated by the Subproblem is not added to the Master problem. The algorithms for generating patterns for periods 1 . . . t − 1, precursor patterns, and for period t + 1 . . . T , following patterns, use enumeration to try to maximize the dual value of each pattern. For the precursors, all pairs of vertices and projected cones are evaluated to see if they are contained in the current pattern, and meet the Global Capacity constraint. The pattern with the highest dual value is selected as the precursor, and the algorithm is repeated to find a precursor to this pattern. Similarly, for each of the following patterns, an estimate of the maximum size of the pattern is generated from the durations of the current pattern and all pairs of vertices and corresponding projected cones that cover the current pattern are evaluated to find the one with the maximum dual value. Generating the Initial Columns The column generation process begins by running the LP relaxation of the Master problem, which requires a feasible solution. A simple initial solution is the same set of m draw points, in a feasible shape, in each of the periods. As an alternative, the algorithms developed for mine sequence generation, discussed above, were adapted to generate a feasible sequence. First the Subproblem is run to generate the mine pattern for period one, then the algorithm is used to generate a sequence of following patterns. 123 Chapter 3. Sequence Optimization in Block Cave Mining Since dual values are not available for the Subproblem or following pattern algorithm, draw point values are substituted in all periods. We called this the greedy/myopic algorithm as it opens the adjacent draw points of highest value in each period. It is discussed below, in Section 3.4.7 Contribution Column Generation is a well used technique in solving large problems, including parallel machine scheduling. In applying Column Generation to Sequence Optimization I developed the mine sequence algorithms necessary to produce a feasible sequence of caves. The algorithms and implementation are new work on the problem of Sequence Optimization and is part of the contribution of this thesis to the research literature. 3.4.7 Greedy/Myopic Algorithm The Greedy/Myopic algorithm is a direct result of the work done on the ColGen model. In order to begin the ColGen algorithm and generate the first set of dual prices, a feasible set of columns must be included in the Master problem. For the first few attempts the initial set of columns was generated by hand. With the modification of the ColGen algorithm to generate mine openings that preceded and followed the mine opening generated by the Subproblem, the opportunity to generate a feasible initial set of columns presented itself. The Greedy algorithm uses the ColGen Subproblem routine to generate a mine opening pattern for the first period, then uses the routine written to generate mine openings for the rest of the periods. Draw point weights are substituted for the reduced costs in the routines. As a result, the highest valuation mine opening is selected for each period resulting in a myopic or greedy algorithm. This is essentially the same method used to generate the initial columns for the ColGen algorithm. Results Since both the ColGen model and the Greedy algorithm use the same cave shape definition as the 2Cone model, the generated caves are similar to those from the 2Cone model. In Section 3.5 the performance of both the ColGen model and the Greedy algorithm on its own are compared to others we developed for this problem. 124 Chapter 3. Sequence Optimization in Block Cave Mining Contribution The idea of a greedy algorithm using a myopic view is well established. In this work, I have created an algorithm that develops a feasible sequence of caves using a greedy approach. The algorithm and implementation is new work on the problem of Sequence Optimization and is part of the contribution of this thesis to the research literature. 125 Chapter 3. Sequence Optimization in Block Cave Mining 3.5 Application of Models to Test Data Sets This section describes the application of the developed models to two new data sets. The models were developed using the four previously described data sets, so there was a risk that the models were tailored to these data sets. Thus we sought out new data sets from Gemcom. They provided us with two data sets coming from a project with Rio Tinto, one of the worlds leading mining and exploration companies. Unfortunately this data is confidential and so is not publicly available. Evaluating the models on two new data sets was an attempt to objectively evaluate and compare the models. The objective was to evaluate the efficiency and practical applicability of the models to aid in the Sequence Optimization process. In addition, the execution speed of the models is compared. All models discussed in the previous sections were run and the solutions and time to solve were compared. Some additional runs were also conducted with constraints on the starting opening pattern constraints to answer additional planning considerations. 3.5.1 Data Sets We now describe the two data sets. RT1 is a mine with a small footprint consisting of 332 draw points to be opened over 10 years. Planners believe that the small size of the footprint reduces the possible combination of sequences and makes it easier to determine a good Sequence Optimization by traditional methods. Running the new models on this data set allows us to demonstrate their potential on a data set not considered difficult to plan for using existing techniques. This mine is already under production which provided a basis of comparison for Rio Tinto. The mining sequence was published graphically (Calder (2000)) but without values so a value comparison is not possible. By contrast, RT2 is a mine with a huge footprint of 6510 draw points to be opened over 20 periods. The size of the data set clearly presents challenges to any solution technique, and so it is a good test of the new models. Data Set RT1 The data set RT1 consists of 332 draw points along 19 vertical tunnels, the longest of which is 30 draw points long; see Figure 3.33 below. 126 Chapter 3. Sequence Optimization in Block Cave Mining Figure 3.33: Data Set RT1 Provided Data: Table 3.6 is a sample of the data available for each draw point as supplied by Gemcom: Draw Point Name: A naming convention used by Gemcom to identify each draw point. The name consist of an “L” followed by a tunnel then either a “W” or an “E” followed by a draw point number X, Y: are the location coordinates of the draw point Tot Dol: this is the total value of the column of ore above the draw point. It is assumed that this is already discounted to represent the net present value (NPV); the details behind the generation of this data were not provided to us as they were considered confidential by Gemcom. Tonnage: the total tonnage (units) of the column of ore above the draw point. Required Data: Duration and Value: Duration and NPV value of each draw point is required for the models. The duration of each draw point is calculated by dividing the amount of ore to be mined by the assumed constant draw rate. 127 Chapter 3. Sequence Optimization in Block Cave Mining Draw Point Name L01W01 L01W02 L01W03 L01W04 L01W05 L01W06 .. . X 13348.78 13348.78 13348.78 13348.78 13348.78 13348.78 .. . Y -24016.7 -24033.7 -24050.7 -24067.7 -24084.7 -24101.7 .. . Tot Dol 11061.73 8689.515 6138.938 7823.849 17987.16 16093.53 .. . Tonnage 452,008 432,104 443,745 444,517 455,378 439,121 .. . Table 3.6: Available Data for Data Set RT1 In this project, the data in column Tonnage was divided by the assumed draw rate which was provided at 120 tons/day. textbf120 tons/day (120 tons/day *365 days/year = 43,800 tons/year) was used. To determine the duration of each draw point the following formula was used Duration = round(Tonnage/43,800 tons per year) The round function converts the result to the nearest integer value (up or down). As a result, the durations ranged from 7 to 17 years and the distribution of durations is shown in Figure 3.34. The total value of each draw point, discounted to the period it is opened, is also required as input to the models. The results in this thesis used the Tot Dol column from the provided data for this input. In this data set there are 37 draw points with negative value that must be opened to complete the mine. As can be seen in Figure 3.35, these draw points with negative value are grouped in tunnels 6, 7 and 8. The highest value draw points are in a cluster just two tunnels over. Maximum active draw points: The downstream processing capacity is captured in the models with a maximum active draw points constraint for each period. The downstream processing capacity was provided in the form of an annual production rate which increases from 3 to 12.41 Mt/y over the first four years of operation. In order to determine the number of active draw points the annual production rate was converted into equivalent number of draw points. Table 3.7 shows the estimated number of draw points in each period, for the constant production rates supplied, 120 tons/day. Discount rate: A discount rate of 10% was supplied by Gemcom Other data: The coordinates of the draw points were converted into tunnel and draw point numbers consistent with the grid system used in the 128 Chapter 3. Sequence Optimization in Block Cave Mining Figure 3.34: Histogram of Draw Point Durations at 120 tons/day for Data Set RT1 Constant Draw Rate 120 tons/day Year 1 Year 2 Year 3 Year 4 Subsequent years Annual production rate Mt 3 6 9 12.41 12.41 Maximum active draw points 70 138 206 284 284 Table 3.7: Maximum Active Draw Points for Data Set RT1 models. In this data set the tunnels were parallel to the “Y” coordinate and the across-tunnel draw points were perpendicular to the tunnels. The data were 129 Chapter 3. Sequence Optimization in Block Cave Mining Figure 3.35: Distribution of Draw Point Values over the Mining Footprint for Data Set RT1 rotated through 90◦ to align the tunnels across the grid to be consistent with our prior usage. Under this transformation, the tunnel numbering is in the reverse order from the tunnel assigned by Gemcom, and tunnel 1 is a combination of L20E and L19W. The “W” draw points take the even numbering and the “E” become the odd. Thus L19W01 becomes tunnel 1, draw point 1 and L20E01 becomes tunnel 1, draw point 2. The Gemcom draw point numbering convention begins with 01 at the start of each tunnel, but due to the shape of the footprint, the ends of the tunnels are not aligned. The new draw point numberings correct for this. 130 Chapter 3. Sequence Optimization in Block Cave Mining Data Set RT2 The data set RT2 consists of 6510 draw points along 56 vertical tunnels, the longest of which is 196 draw points long; see Figure 3.36 below. Figure 3.36: Data Set RT2 Provided Data: Table 3.8 is a sample of the data available for each draw point as supplied by Gemcom: 131 Chapter 3. Sequence Optimization in Block Cave Mining DPTNAME P01 5621W P01 5524W P01 5524E P01 5523W P01 5523E .. . XCOORD 495460.3 495408.8 495417.9 495419.1 495428.3 .. . YCOORD 3683259 3683232 3683226 3683249 3683243 .. . Tot Dol -214.0677338 -200.9185181 -435.5866394 158.174881 10.73392391 .. . TON14M 116,799 75,699 101,510 74,210 86,200 .. . Table 3.8: Available Data for Data Set RT2 DPT Name: A naming convention used by Gemcom to identify each draw point. The name begins with a letter and two digits followed by an underscore then 4 digits and either “W” or “E”. The meaning of the segment before the underscore is not clear but could refer to a section of the mine, but it appears that the underscore is followed by a tunnel number and a draw point number. While the tunnel numbers appear to be consistent across sections, the draw point numbering is not. XCOORD, YCOORD: are the location coordinates of the draw point Tot Dol: this is the total value of the column of ore above the draw point. TON14M: the total tonnage (units) of the column of ore above the draw point. Required Data: Duration and Value: Duration and NPV value of each draw point is required for the models. The duration of each draw point is calculated by dividing the amount of ore to be mined by the assumed constant draw rate. In this project, the data in column TON14M was divided by the assumed draw rate of 270 tons/day. This value for the draw rate was chosen in conjunction with the maximum active draw points to ensure that the entire footprint could be mined in 20 periods. 270 tons/day (270 tons/day *365 days/year = 98,550 tons/year) was used. To determine the duration of each draw point the following formula was used: Duration = round(TON14M/98,550 tons per year) As a result, the durations ranged from 1 to 5 years, with a single draw point at 5 years, and the distribution of durations is shown in Figure 3.37. 132 Chapter 3. Sequence Optimization in Block Cave Mining Figure 3.37: Histogram of Durations for Data Set RT2 The total value of each draw point, discounted to the period it is opened is also required as input to the models. The results in this thesis used the Tot Dol column from the provided data and the calculated duration to calculate this input. The total value is calculated using the PV(present value) function in Microsoft Excel. It is the PV of Duration installments of (Tot Dol/Duration) payments at the discount rate. Maximum active draw points: The downstream processing capacity is captured in the models with a maximum active draw points constraint for each period. The maximum active draw points values were chosen along with the daily draw rate to ensure that all draw points could be mined in 20 periods and are shown below in Table 3.9. Discount rate: A discount rate of 10% was supplied by Gemcom Other data: The coordinates of the draw points were converted into tunnel and draw point numbers consistent with the grid system used in the models. In order to establish the grid, the data was projected onto a new coordinate set. The new axis R runs parallel to the tunnels and is obtained by rotating the original coordinates through 50◦ . This is the rotation required 133 Chapter 3. Sequence Optimization in Block Cave Mining Constant Draw Rate 270 tons/day Year 1 Year 2 Year 3 Year 4 Subsequent years Annual production rate Mt 28,382,400 54,793,800 82,387,800 109,587,600 109,587,600 Maximum active draw points 288 566 836 1112 1112 Table 3.9: Maximum Active Draw Points for Data Set RT2 to align the first tunnel along the horizontal axis. The new S axis runs in the across draw point axis and is perpendicular to the R axis. It is obtained by rotating the original coordinates through 60◦ . The procedure followed was as follows 1. Choose a point as the “origin”, it will have new coordinates (0,0). Denote the original coordinates of this point as (Xo ,Yo ) 2. calculate the distance from each draw point to the origin using the X and Y coordinates 3. calculate the new R coordinate as distance × sin{90 − 50 − arcsin{(Y − Yo )/distance}} 4. calculate the new S coordinate as distance×sin{arcsin{(Y −Yo )/distance}− 60} 134 Chapter 3. Sequence Optimization in Block Cave Mining 3.5.2 Models Five sequence optimization models were run on data sets RT1 and RT2. Four of the models were written in the linear/integer programming framework using the language AMPL and solved using CPLEX solver Version 9.1. The fifth model, a greedy algorithm, was written in Visual Basic. All five models used the same input data, and the same layout of draw points. The models were run on a Intel Pentium (R) D CPU 3.20 GHz machine with 1.99 GB of RAM. Basic Model: This is the original “across and along” model, i.e., contiguity along and across each tunnel, described in Section 3.4.3. The constraints ensure that in each period, if a tunnel contains open or active draw points that they form a single contiguous segment. There are also constraints to enforce a similar contiguity constraint in the across-tunnel dimension. Recall from Section 3.4.3 that earlier runs have shown that enforcing contiguity in only two directions does not guarantee a single contiguous cave in each period. Maximize: NPV Subject to : Start Once : (3.11) Global Capacity : (3.12) Contiguous along tunnels: (3.13) - (3.15) Contiguous across tunnels: (3.16) - (3.19) Table 3.10 shows the number of binary variables and constraints for the Basic model. In addition to the set of variables that indicate which draw points are open in each period, there are two sets of binary variables which denote the first open draw point along and across each tunnel. There are four contiguity constraints for each draw point in each period plus two for each tunnel in each period. These form the majority of the constraints for the model. For Data Set RT1 there are a total of almost 14 000 constraints, and this grows to almost 530 000 for Data Set RT2. Malkin Model: This model is a formulation based on the ideas of Malkin and Wolsey, Malkin (2006). For a convex mine footprint, the model, described in Section 3.4.4, guarantees a single cave in each period by defining a fixed centre point to the cave and only allowing draw points to open if they are connected by open or active draw points to the centre point. A 135 Chapter 3. Sequence Optimization in Block Cave Mining nDP T S × T × 3 Binary Variables Constraints Start Once Global Capacity Contiguous Along Contiguous Across Total Constraints nDP T S T nDP T S × T × 2 + nT N LS × T nDP T S × T × 2 + nT N LS × T nDP T S + T + nDP T S × T × 4 + nT N LS × T × 2 RT1 9 960 RT2 390 600 332 10 6 970 6 970 13 922 6 510 20 261 520 261 520 529 570 Table 3.10: Number of Variables and Each Type of Constraints for Basic Model Runs on Data Set RT1 and RT2 fixed set of axes are constructed from the centre draw point to help with the formulation. The fixed centre point model was used as the increased data preparation for the “wandering axes” variation did not result in improved results in the developmental runs in Section 3.4.4. Maximize: NPV Subject to : Start Once : (3.11) Global Capacity : (3.12) Column Connected: (3.24)-(3.25) Row Connected: (3.26)-(3.29) Table 3.11 shows that the Malkin model is considerably smaller than the Basic number, both in number of binary variables and number of constraints. The Malkin model has one third of the variables, and approximately one half of the constraints of the Basic model. For Data Set RT1 there are approximately 7 000 constraints, and approximately 130 000 constraints for Data Set RT2. 2Cone Model: The 2 vertices 2 cone model, described in Section 3.4.5, is an attempt to get a single cave that has a diamond shape. The cave is the intersection of two cones or triangles each determined by the decision variable representing the vertex of the cone. Figure 3.31 illustrates the concept. 136 Chapter 3. Sequence Optimization in Block Cave Mining nDP T S × T RT1 3 320 RT2 130 200 nDP T S T nDP T S × T nDP T S × T nDP T S + T + nDP T S × T × 2 332 10 3 320 3 320 6 982 6 510 20 130 200 130 200 266 930 Binary Variables Constraints Start Once Global Capacity Column Connected Row Connected Total Constraints Table 3.11: Number of Variables and Each Type of Constraints for Malkin Model Runs on Data Set RT1 and RT2 Maximize: NPV Subject to : Start Once : (3.11) Global Capacity : (3.12) One top vertex in each period: (3.31) One lower vertex in each period: (3.32) Cone Definition : (3.33) - (3.34) Open Intersection of Cones : (3.35) - (3.37) Previous Cave Contained in Current Cave: (3.38) - (3.39) The 2Cone model is the largest of the three integer programming models, as can be seen in Table 3.12. There are five sets of binary variables, including the two for the vertices over the extended footprint. This model has approximately 21 000 binary variables for Data Set RT1, and just over 1 000 000 for Data Set RT2. A large number of constraints are needed to define the cones and their intersection resulting in a model with approximately twice as many constraints as the Basic model. ColGen: The column generation technique, described in Section 3.4.6, uses a LP model to choose from a set of cave patterns to form a maximum value, feasible sequence over all periods. A separate algorithm, not limited by the LP/IP framework, uses the dual prices generated in the LP to add new candidates to the set of cave patterns. In this case, the algorithm was based on the 2 Cone concept and found the highest value candidates by enumeration. 137 Chapter 3. Sequence Optimization in Block Cave Mining Binary Variables Constraints Start Once Global Capacity One of Each Vertex Cone Definition Cone Intersection Previous Cone Total Constraints nDP T S × T × 5+ eT N LS 2 × T × 2 nDP T S T T ×2 nDP T S × T × 2 nDP T S × T × 3 nDP T S × (T − 1) × 2 nDP T S + T × 3+ nDP T S × T × 7 − nDP T S × 2 RT1 21 100 RT2 1 035 160 332 10 20 640 960 976 938 6 510 20 40 260 400 390 600 247 380 904 950 6 9 5 22 Table 3.12: Number of Variables and Each Type of Constraints for 2Cone Model Runs on Data Set RT1 and RT2 Greedy : The greedy model uses the cave generation algorithms from the ColGen model to generate an initial set of cave patterns that form a feasible sequence and is described in Section 3.4.7. Since there were no dual values available, the total value of the draw point is used. 138 Chapter 3. Sequence Optimization in Block Cave Mining 3.5.3 Results Results on Data Set RT1 Table 3.13 compares the total NPV and number of caves from running the models discussed above on Data Set RT1. As the Basic model has the most relaxed constraints, it is expected to give the highest NPV so will be used as a benchmark for the other models. Basic Malkin 2Cone ColGen Greedy Total NPV $ 3,507,820 $ 3,457,480 $ 3,314,810 $ 3,333,938 $ 3,333,520 # of Caves >2 1 1 1 1 % of max NPV 100 98.6 94.5 95 95 Table 3.13: Results for Data Set RT1 The highest NPV is from the Basic model, however inspection of the results shows that this model does not actually result in a single contiguous cave. Note that the Malkin result is (3,457,480/3,507,820 = 0.986) 98.6% of the Basic model result and the 2Cone result is (3,314,810/3,507,820 = 0.945) 95% of the Basic model result. Somewhat surprising is the success of the Greedy model to achieve such a good result;(3,333,520/3,507,820 = 0.950) 95.0% of the Basic model result. Subsequent iterations of the ColGen model, from the Greedy initial solution, did result in improvement of the objective value, but it is almost negligible. (3,333,938/3,507,820 = 0.950) 95.0% of the Basic model result was achieved with ColGen model. Table 3.14 reports the differences in active cave size over the ten periods for the different models. There is very little difference between the models with the exception of the Basic model. There are 37 draw points with negative value that must be mined for structural reasons. Some of the models may delay these until the last period and therefore discount this loss. This may explain why the Basic model only opens 132 of the possible 140 draw points in the second period. As can be seen from the results in Table 3.15, the models took from 4 seconds to almost an hour to run. Listed are both the times for the complete set of Malkin models, with a run with each of the 332 draw points as the fixed centre point, and the single run that produced the highest value pattern. The Basic model ran far faster than the 2Cone and Malkin models, 139 Chapter 3. Sequence Optimization in Block Cave Mining Number of draw points active in each period Basic Malkin 2Cone ColGen Greedy 1 70 70 70 70 70 2 132 140 140 140 140 3 205 206 198 204 205 4 284 284 284 284 284 5 284 284 284 284 284 6 284 284 284 284 284 7 284 284 284 284 284 8 284 284 283 276 279 9 283 284 283 284 279 10 254 241 223 232 232 Table 3.14: Cave Development Results for Data Set RT1 Total NPV Solution Time Basic $ 3,507,820 4 seconds 2Cone $ 3,314,810 3 050 seconds ≈ 50 minutes Malkin: all draw points 2 416 seconds ≈ 40 minutes Malkin (maximum value) $ 3,457,480 2 seconds “Malkin: all draw points” is the aggregation of 332 runs of the Malkin model, once for each draw point as the fixed centre point. ‘Malkin (maximum value)” is the single run that produced the highest value Table 3.15: Solution Times: however the resulting multiple caves are undesirable. In this case, the 2Cone model took 10 minutes longer to run than the Malkin model, but produced a higher NPV. We now look at the mine shapes generated by solving these models on data set RT1. The figures show how well the results conform to the objective of a single cave as well as the preferred diamond shape. We also compare the progression of the caves over time to see how the results from the different models vary. 140 Chapter 3. Sequence Optimization in Block Cave Mining Basic: The following sequence was determined using the Basic model. As has been found in the past, this model does not guarantee a single contiguous cave. The results show that there is a high value group of draw points on the right side that is chosen in the first period. It also shows that there is a valuable group on the left side that is also chosen in the first period, despite not being adjacent to the group on the right. Figure 3.38: Results of Basic Model on Data Set RT1 - Period 1 141 Chapter 3. Sequence Optimization in Block Cave Mining Figure 3.39: Results of Basic Model on Data Set RT1 - Period 2 Figure 3.40: Results of Basic Model on Data Set RT1 - Period 3 142 Chapter 3. Sequence Optimization in Block Cave Mining Figure 3.41: Results of Basic Model on Data Set RT1 - Period 4 Figure 3.42: Results of Basic Model on Data Set RT1 - Period 10 143 Chapter 3. Sequence Optimization in Block Cave Mining Malkin: The Malkin model was run repeatedly, changing the centre point over all draw points to find the centre point yielding the highest NPV. The run using this centre point was chosen for reporting and the pattern is shown below. As with the Basic model, in the first period draw points on right are chosen with a movement towards the left in the second period. The caves for years 2 and 3 extend along the across-tunnel axis and do not look diamond like. We can see that the constraints are met, in each period there are no gaps between any open draw point and either axis. In the first period, the cave is a diamond shape, but the axes of the diamond do not align with the designated axes. In the second and third period, the lack of symmetry around the axes results in a non-diamond shape. Figure 3.43: Best Results of Malkin Model for Data Set RT1 - Period 1 144 Chapter 3. Sequence Optimization in Block Cave Mining Figure 3.44: Best Results of Malkin Model for Data Set RT1 - Period 2 Figure 3.45: Best Results of Malkin Model for Data Set RT1 - Period 3 145 Chapter 3. Sequence Optimization in Block Cave Mining Figure 3.46: Best Results of Malkin Model for Data Set RT1 - Period 4 Figure 3.47: Best Results of Malkin Model for Data Set RT1 - Period 10 146 Chapter 3. Sequence Optimization in Block Cave Mining 2Cone: The 2Cone model generated a single cave in a diamond pattern. The first period cave is very similar to that chosen in the Basic model and the highest value Malkin model. After choosing the high value draw points in the first period, the model moved to the right before moving left through the negative value draw points. Figure 3.48: Results of 2Cone Model for Data Set RT1 - Period 1 147 Chapter 3. Sequence Optimization in Block Cave Mining Figure 3.49: Results of 2Cone Model for Data Set RT1 - Period 2 Figure 3.50: Results of 2Cone Model for Data Set RT1 - Period 3 148 Chapter 3. Sequence Optimization in Block Cave Mining Figure 3.51: Results of 2Cone Model for Data Set RT1 - Period 4 Figure 3.52: Results of 2Cone Model for Data Set RT1 - Period 10 149 Chapter 3. Sequence Optimization in Block Cave Mining ColGen: The ColGen Model also generated a single cave in a diamond pattern, but the first period cave was not as centrally located as the other models. The first period cave includes some high value draw points but extends over to the right hand edge. Over the remaining periods the cave is extended left, then ends with draw points at the bottom of the layout. Figure 3.53: Results of ColGen Model for Data Set RT1 - Period 1 150 Chapter 3. Sequence Optimization in Block Cave Mining Figure 3.54: Results of ColGen Model for Data Set RT1 - Period 2 Figure 3.55: Results of ColGen Model for Data Set RT1 - Period 3 151 Chapter 3. Sequence Optimization in Block Cave Mining Figure 3.56: Results of ColGen Model for Data Set RT1 - Period 4 Figure 3.57: Results of ColGen Model for Data Set RT1 - Period 10 152 Chapter 3. Sequence Optimization in Block Cave Mining Greedy: As noted above, there is very little difference between the Greedy algorithm results, and the improvements made by subsequent iterations of the ColGen model. Differences do not appear until year 4 and only on the edges of the cave. Figure 3.58: Results of Greedy Model for Data Set RT1 - Period 1 153 Chapter 3. Sequence Optimization in Block Cave Mining Figure 3.59: Results of Greedy Model for Data Set RT1 - Period 2 Figure 3.60: Results of Greedy Model for Data Set RT1 - Period 3 154 Chapter 3. Sequence Optimization in Block Cave Mining Figure 3.61: Results of Greedy Model for Data Set RT1 - Period 4 Figure 3.62: Results of Greedy Model for Data Set RT1 - Period 10 155 Chapter 3. Sequence Optimization in Block Cave Mining Comparison of first two periods over all models This section is included to facilitate comparison of the results for the first two periods across all models. As can be seen below, all models start with the draw points in tunnels 13 through 19; tunnels numbered from left to right (Figure 3.33). The results are very similar, with the exception of those for the Basic model which moves to tunnels 1 through 10 in the second period and results in multiple caves. While all models created a diamond-shaped cave in the first period, the 2Cone, ColGen and Greedy models are more successful at generating diamond-shaped caves in the second period. Combining this with the NPV value results suggests the 2Cone model performs better than the other models. Figure 3.63: Opening Patterns for First Two Periods for Basic Model 156 Chapter 3. Sequence Optimization in Block Cave Mining Figure 3.64: Opening Patterns for First Two Periods for Malkin Model Figure 3.65: Opening Patterns for First Two Periods for 2Cone Model 157 Chapter 3. Sequence Optimization in Block Cave Mining Figure 3.66: Opening Patterns for First Two Periods for ColGen Model Figure 3.67: Opening Patterns for First Two Periods for Greedy Model 158 Chapter 3. Sequence Optimization in Block Cave Mining Sensitivity: This data set was chosen for this project as it was thought by planners that its relatively small size provided few opening sequences and would not lead to much sensitivity to sequence patterns. To test sensitivity to sequence, some alternate scenarios, intended to be different from those generated above, were run as follows: 1. Start at LHS (2Cone): Force mining to start on the left-hand side of footprint, i.e., Tunnel 5 if counted from the LHS. This was done by forcing all the draw points in tunnel 5 to start in period 1 or period 2. All the results above started with the high value draw points in the centre or right of centre of the layout. The choice of tunnel 5 is distinct from the centre point of the model results above, and allows for a diamond shaped first period cave, without running into the edge of the footprint. 2. Start on Edge (First or Last Tunnel) (2Cone): Force mining to start on an edge (Tunnel 1 or Tunnel 19) instead of the centre of the footprint. As mentioned in Section 3.4.4, some mining sequences originate at the edge of a mining footprint so that the equipment used to remove the mined rock does not cross an active mining front. If the sequence does not start on a footprint edge, often the mined rock must be carried on a level below the mining to cross the mining front. Sequences that start on a footprint edge are of interest to understand the tradeoffs between these two approaches. 3. Start on Edge (Any) (2Cone): Force mining to start on an edge (Tunnel 1 or Tunnel 19 or “top” or “bottom”) 4. Look for Minimum NPV (2Cone): run the model to minimize NPV. Since the easiest way to do this is delay starting all draw points until the end, or not open any at all, a new constraint was added to require at that least “maximum number of active draw points - 4” must be active in each period. This sequence can provide a measure of the worst case scenario of a badly chosen sequence, and demonstrate the potential value of the models on this data set. 5. Lowest value Malkin: The Malkin model solution with the fixed centre point yielding the lowest NPV. Again, this is a worst case scenario to show effect of chosen sequence on value. 159 Chapter 3. Sequence Optimization in Block Cave Mining 6. Force mining to start in the middle (2Cone). This was done by forcing draw points in tunnels 10 and 11 (numbered from LHS) to start in period 1 or period 2. It was chosen to reproduce the published sequence (Calder (2000)) for this data set. Figure 3.68 below shows the sequence that was published for the data set. Figure 3.68: Published Sequence Optimization for Sample Data Set Table 3.16 reports the results of the sensitivity runs. The 2Cone model from the earlier runs was used as a benchmark to compare the values of these alternate runs, based on its high NPV value and diamond-shaped caves. The results in Table 3.16 do show that the total NPV of caves produced on this data set is sensitive to the opening sequence with many of the alternates 10% below the 2Cone value from the original runs. The results also show a variation in solution times. The individual Malkin run is faster than most of the 2Cone results. The more constrained the 2Cone model is, the faster the run completed. While the longest 2Cone solution time of 8 hours is not fast by any measure, it would be possible to run this model overnight, especially for the scale of costs involved with mine planning. There is not much variation in the number of active draw points in each cave between the runs; some runs did not open a full 70 draw points in the first period. See Table 3.17 The caves from the first two years are presented below in Figure 3.69 through Figure 3.73 to confirm that the objectives of each scenario were met. 160 Chapter 3. Sequence Optimization in Block Cave Mining Total NPV Basic 2Cone Start at LHS (2Cone) Start on Edge (First or last tunnel)(2Cone) Start on Any Edge (2Cone) Minimize (2Cone) Malkin (low value) Start in centre (2Cone) % of 2Cone NPV 105.8 100 91.3 88.4 4 sec ≈ 0.001 h 3 050 sec ≈ 0.8 h 694 sec ≈ 0.2 h 2 598 sec ≈ 0.7 $ 3,304,930 99.7 1 564 sec ≈ 0.4 h $ 2,809,050 $ 3,063,040 $ 3,136,930 84.7 92.4 94.6 29 132 sec ≈ 8 h 303 sec ≈ 0.08 h 149 sec ≈ 0.04 h $ $ $ $ 3,507,820 3,314,810 3,025,550 2,932,900 Solution time Table 3.16: Results of Sensitivity Runs: Basic 2Cone Start at LHS Start on Edge (First or last tunnel) Start on Edge (Any) Minimize Malkin (low value) Start in centre 1 70 70 66 67 Number of draw 2 3 4 132 205 284 140 198 284 135 202 284 137 202 284 points active in each period 5 6 7 8 9 284 284 284 284 283 284 284 284 283 283 284 284 284 284 283 284 284 284 284 279 10 254 223 234 267 69 138 205 284 284 284 284 282 283 237 67 70 137 140 202 206 280 284 280 284 280 284 280 284 260 284 238 284 251 268 70 137 205 283 283 283 283 279 263 205 Table 3.17: Cave Sizes in Sensitivity Runs: 161 Chapter 3. Sequence Optimization in Block Cave Mining Figure 3.69: Results of 2Cone Model Starting at Tunnel 5 Figure 3.70: Results of 2Cone Model Starting on Edge Tunnel 162 Chapter 3. Sequence Optimization in Block Cave Mining Figure 3.71: Results of 2Cone Model Starting on Any Edge Figure 3.72: Results of 2Cone Model Minimized 163 Chapter 3. Sequence Optimization in Block Cave Mining Figure 3.73: Results of Minimum Malkin 164 Chapter 3. Sequence Optimization in Block Cave Mining Start at centre This run was an attempt to reproduce the published sequence for this data set. Figure 3.74 below is an overlay of the sequence for the first 3 periods with the published sequence. While the published sequence starts in the middle of tunnel 5, as opposed to the bottom in the model result, both move to open the draw points in the tunnels to the right before moving towards the tunnels on the left. The value of the published sequence is not available for comparison. Figure 3.74: Overlay of Published Results and First 3 periods of Generated Results The following figures show the generated results without the overlay of the published results. As mentioned, the first period cave is in the middle tunnels, but touches the lower edge of the mining footprint. In the second period it moves to the left, and after that grows in both directions. 165 Chapter 3. Sequence Optimization in Block Cave Mining 166 Chapter 3. Sequence Optimization in Block Cave Mining Figure 3.75: Results of 2Cone Model Starting at Centre 167 Chapter 3. Sequence Optimization in Block Cave Mining Results on Data Set RT2 Table 3.18 reports the results of the models discussed above. Due to the length of time of each run it was not feasible to run the Malkin model on all 6510 possible centre points. Instead a sampling was taken based on the draw points selected in the first period of the Basic model. The 2Cone model did not run to completion on this data set. After two days of calculations the results were not returned and it appeared the maximum memory was exceeded. Similarly, there are no results for the ColGen model as there were no results after several days of running and no intermediate results were returned. Total NPV Basic Malkin (29,80) Malkin (23,85) Malkin (26,95) Malkin (35,70) Malkin (25,100) Malkin (39,50) 2Cone ColGen Greedy $ $ $ $ $ $ $ 11,062,600 10,934,200 10,877,400 10,861,900 10,831,700 10,803,000 10,706,000 $ 10,441,300 Number of Caves 2 1 1 1 1 1 1 % of Basic NPV 100 98.8 98.3 98.2 97.9 97.6 96.8 1 94.4 Solve Time (sec) 55 719.6 5 372.8 4 299.7 3 649.2 33 128.1 4 867.4 25 849.4 infinite infinite > 6 days ≈ ≈ ≈ ≈ ≈ ≈ ≈ 19 hr 1.5 hr 1 hr 1 hr 9 hr 1 hr 7 hr Table 3.18: Results for Data Set RT2: The highest NPV is from the Basic model. However, inspection of the results, see Figure 3.76 shows that this model produces two contiguous caves. The Malkin run with the highest result is (10,934,200/11,062,600= 0.988) 98.8% of the Basic model result and the Malkin run with the lowest result is (10,706,000/11,062,600= 0.968) 96.8% of the Basic model result, demonstrating that the choice of centre point can improve objective value. As with data set RT1, the Greedy model achieved a very good result achieving (10,441,300/11,062,600= 0.944) 94.4% of the Basic model result. The models took considerably longer to run on this data set than on Data Set RT1, as expected. The Basic model took almost a day to complete, and the Greedy model almost a week. Since the ColGen model initializes with the Greedy model, it is not surprising that it did not complete. The Greedy algorithm evaluates the caves formed by each pair of draw points, 168 Chapter 3. Sequence Optimization in Block Cave Mining roughly 42 million combinations for data set RT2. Clearly there is some room to improve this algorithm and as a result, improve the speed of the ColGen model. While individual Malkin runs appear much quicker, many runs would be needed to find the best centre point. Table 3.19 shows the number of active draw points in each period for the Basic, the highest and lowest total value from the sample of Malkin runs, and the Greedy algorithm. There are 201 draw points with negative value that must be mined for structural reasons. Some of the models may delay these until the last period and therefore discount this loss. This may explain why the models wait until the last period to open the remaining draw points. We now look at the mine shapes generated by solving these models on data set RT2. Inspection of the figures provides a count of caves formed. We also compare the shape and progression of the caves over time to see how the results from the different models vary. 169 Chapter 3. Sequence Optimization in Block Cave Mining Number of draw points active in each period 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 Basic Malkin Maximum of sample (29,80) Malkin Minimum of sample (5,80) Greedy 288 566 835 1090 1099 1109 1095 1103 1080 1064 1077 1155 1139 1184 1108 792 22 3 1 163 288 566 831 1087 1086 1117 1102 1154 1126 1073 1030 1113 1113 1172 1126 846 21 2 0 116 278 543 794 1067 1067 1075 1067 1081 1082 1092 1048 1055 1137 1159 1190 1068 42 0 0 124 288 566 836 1117 1106 1102 1086 1112 1108 1089 1095 1095 1104 1134 1114 906 189 42 35 33 Table 3.19: Cave Development Results for Data Set RT2: 170 Chapter 3. Sequence Optimization in Block Cave Mining Basic: The following sequence was determined using the Basic model. As has been found on other data sets, this model did not return a single contiguous cave in the first period. The results show that there is are two high value groups of draw points on the right side that are chosen in the first period. It also shows that there is a valuable group on the bottom left side that is not adjacent to those chosen in the first period but which the model seeks in the fourth period. 171 Chapter 3. Sequence Optimization in Block Cave Mining Figure 3.76: Results from Basic Model on Data Set RT2 172 Chapter 3. Sequence Optimization in Block Cave Mining Malkin: As each run of the Malkin model took between 1 and 9 hours it was not feasible to run the model for each of the 6510 draw points as the starting point. Instead 6 runs were done, each starting in one of the two caves selected in the first period of the Basic model run. Of the runs done, the maximum total NPV was centre point (29,80) with NPV $10,934,200, which was 2.6% higher than the minimum at centre point (39,50) with NPV $10,706,000. While there is not much difference in the total NPV between the Malkin runs, the two centred on (35,70) and (39,50) took much longer time to run than the others. These two runs were the only two centred in the right hand cave generated in first period of the Basic model. Starting Tun Dpt Total NPV 29 23 26 35 25 39 $ $ $ $ $ $ 80 85 95 70 100 50 10,934,200 10,877,400 10,861,900 10,831,700 10,803,000 10,706,000 Number of Caves Solve Time (sec) 1 1 1 1 1 1 5 372.8 4 299.7 3 649.2 33 128.1 4 867.4 25 849.4 ≈ 1.5 hours ≈ 1 hour ≈ 1 hour ≈ 9 hours ≈ 1 hour ≈ 7 hours Table 3.20: Results from Selected Runs of Malkin Model on RT2 Due to similarities in the caves generated, graphical results from only the runs centred at (29,80) and (39,50) are shown below. 173 Chapter 3. Sequence Optimization in Block Cave Mining Tunnel 29, Draw point 80 This run was started in the left hand cave selected in the first period of the Basic model run. As with the Basic model, the cave moves to the right, then heads towards the bottom left of the footprint in period five. 174 Chapter 3. Sequence Optimization in Block Cave Mining Figure 3.77: Results from Malkin Model (29,80) for Data Set RT2 175 Chapter 3. Sequence Optimization in Block Cave Mining Tunnel 39, Draw point 50 This run was centred on the upper right corner of the right hand cave selected in the first period of the Basic run on this data set. As with the other Malkin run from the right hand side, Malkin (35,70), the cave moves to the left, almost reaching the edge by period five. The cave then moves to the bottom left of the footprint, which is notably different from all runs other than Malkin (35,70). Figure 3.78: Results from Malkin Model (39,50) 176 Chapter 3. Sequence Optimization in Block Cave Mining Greedy: The Greedy model returned a diamond-shaped sequence that begins in the same area as the Basic model. The cave grows to the left hand edge in the early periods but the tight shape restrictions prevent it from reaching the bottom left until much later in the run. 177 Chapter 3. Sequence Optimization in Block Cave Mining Figure 3.79: Results from Greedy Algorithm on Data Set RT2 178 Chapter 3. Sequence Optimization in Block Cave Mining 3.5.4 Conclusions The new models demonstrated their contribution on the two data sets. On both data sets the models not only identified high value opening sequences, but also demonstrated sensitivity of total NPV to opening sequences. While the Basic model did not produce a feasible opening in the early periods, it did provide an upper bound on the NPV as well as indicate the sections of the footprint with high value. The Malkin model ran well on both data sets, and provided feasible opening sequences as well as an opportunity to test the sensitivity of the results to different initial opening constraints. The 2Cone model was effective on the smaller data set, but proved to be too large for the larger data set. On the smaller data set, the 2Cone model also proved an effective tool for runs to investigate the sensitivity of the NPV to restrictions on the opening sequence. The Greedy algorithm produced results that were close to those of the other models for Data Set RT1. For Data Set RT2 the total NPV was close to the maximum, but the generated caves were noticeably different from those from the Basic and Malkin model runs. Finally, the ColGen algorithm produced results similar to those of the 2Cone model for Data Set RT1, but, like the 2Cone model, was not able to manage the large number of draw points in Data Set RT2. From these two data sets, it seems that no one model provides the solutions that a practitioner seeks, but instead each contribute to the understanding of the problem. The Basic model provides an upper bound on the total NPV and indicates high value starting areas. Unfortunately, the frequent generation of multiple caves makes the solution impractical. In large data sets, running the Malkin model on all centre points may be infeasible due to time constraints. However, individual runs can be useful in investigating sensitivity to centre point or other constraints. The 2Cone model also proved useful, and ran in a reasonable time on the smaller data sets. The Greedy algorithm performed reasonably well but as it was not faster than the others it seems less useful than the other models. 3.6 Conclusions In the preceding sections we presented five model formulations to assist the mine planner with the sequence optimization process for planning a block cave mine. The models each have strengths and weaknesses which will be discussed below. This is a complex problem, and this work has not found a simple answer. Instead, a combination of the presented models is proposed to provide insight to the planner. 179 Chapter 3. Sequence Optimization in Block Cave Mining Model Number of (Binary) Decision Variables Number of Constraints Basic nDP T S × T × 3 nDP T S + T + nDP T S × T × 4 +nT N LS × T × 2 Malkin nDP T S × T nDP T S + T + nDP T S × T × 2 2Cone nDP T S × T × 5+ eT N LS 2 × T × 2 nDP T S + T × 3 + nDP T S × T × 7 −nDP T S × 2 Table 3.21: Integer Programming Formulation Comparison All models share the common goal of creating a single, contiguous cave, the size of which is regulated over time by a production capacity. The models differ on how they achieve the single cave, as this proved a challenging constraint to implement. All models require the same data preparation, namely the numbering of draw points on a two-dimensional grid with one axis aligned with the planned tunnels, and a second axis perpendicular to the first. The first axis is referred to as “along-tunnel” and the second as “across-tunnel”. Three of the models, Basic, Malkin and 2Cone, use general integer programming methods for the formulation. The ColGen model, based on the column generation technique and a Greedy algorithm round out the set of models. Table 3.21 compares size of the three models based on integer programming formulations. The Malkin model is the smallest both in number of decision variables, and number of constraints. Basic Model The Basic model is the most general of the models. The cave is allowed to grow as long as, in each period, the opened draw points form a single, contiguous group along each tunnel, and a single contiguous group in the across-tunnel direction. With this first model, there is no attempt to create the desired diamond shape, only a single, contiguous cave in each period. The single, contiguous group is created by specifying a direction to travel along each axis, then, for each row, naming an open draw point following a closed draw point the “First opened” and allowing only one in each period. 180 Chapter 3. Sequence Optimization in Block Cave Mining Thus, in each row of each axis there is only one transition from closed to open draw points, creating a contiguous group along each row. Running this model revealed that it does not guarantee a single cave in each period. Apparently, for some data sets, the constraints are needed along more axes than the two defined for this work. Despite the failure of the formulation to produce a feasible sequence, there is value to this model. The formulation may return a single cave in each period, or even a sequence that can be modified into something useful by the mine planner. If multiple caves are returned, the result can be used as an upper (unreachable) bound against which the sequences created by the skilled mine planner, and other models can be evaluated. The caves generated from the early periods can also be used to highlight high value areas which may be used as input into other models presented. The Basic model ran to completion on the two data sets tested. It ran in 4 seconds on the smaller data set, but this rose to 19 hours for the larger one. To some, 19 hours may seem unexpectedly long, but in practice it may be acceptable to run it over the weekend if it does not need to be repeated often. Malkin Model The Malkin model was based on the 2D Integral Formulation proposed by P. Malkin and L. Wolsey. In this formulation, a centre point is specified which locates a dominant axis in each of the across-tunnel and along-tunnel directions. In each period, a draw point can open if those between it and each dominant axis are open. The formulation is the smallest of those presented as it is simple and does not require any decision variables in addition to the set describing which draw points are open. For a convex mine footprint, this formulation guarantees a single, contiguous cave in each period. While the dominant axes can suggest a diamond shape, there is nothing in this model that evenly distributes open draw points on either side of the axes. In trials on our data set, this leads to a variety of shapes from a square with the origin of the axes at one corner, to irregular shapes with long, thin protrusions along an axis. Additional work would be required by the mine planner to obtain a usable sequence of caves, and this work would change the total NPV of the sequence. The obvious challenge with this model is the task of choosing a centre point. In this work we were able to run the model on all centre points in less than an hour for a smaller data set. On the larger data set, individual runs took over an hour making it infeasible to run all 6510 centre points. For the 181 Chapter 3. Sequence Optimization in Block Cave Mining larger data set we ran the model on a selection of centre points chosen from the caves formed by the Basic model in the first period. The best result from this selection of runs was evaluated against the upper bound formed by the Basic model result. On the smaller test data set, the best Malkin model achieved 98.6% of the Total NPV of the Basic model, while the best of the selection of runs on the larger data set achieved 98.8% of the Total NPV of the Basic model. The test data sets showed that the strengths of this model are that it can return sequences for a single cave, and they can be close to the upper bound created by the Basic model. This is in contrast to some of the other models that did not return a solution on the larger data set. This model does provide the mine planner with a way to compare a variety of starting points in the mine footprint. For example, production method constraints may require the cave to begin along an edge of the footprint. Running the model on a selection of points around the perimeter could identify which edge is most favorable. In theory all of the integer programming models can perform this task, but this model is the smallest and individual runs were faster than Basic model runs on our trial data sets. 2Cone Model In the 2Cone model the cave for each period is defined by the intersection of two cones, or triangles. Each cone has a draw point as its vertex, then grows from the vertex by two draw points per tunnel. Other than the constraints that ensure previous caves are contained within the current one, the vertices are unrestricted across the mining footprint. This formulation will provide a single cave in each period on a convex data footprint. The shape of each cave is much closer to a diamond shape than those produced by the Basic model and Malkin model, especially in the earlier periods. The definition of each of the cones uses two sets of decision variables, one for the vertex, and one for the cone. As a result, this was the largest of the integer programming formulations. On the smaller data set, this model took slightly longer to run than running the Malkin model on all centre points. The value of the result was lower than that of the best Malkin model run, 94.5% of the Basic model result, suggesting that the shape constraint was more restricting than the fixed dominant axes constraint. This model did not complete on the larger data set. The strength of this model is in the creation of diamond shaped caves. This formulation has specified that the cave grows at two draw points per 182 Chapter 3. Sequence Optimization in Block Cave Mining tunnel, a restriction which may not suit all, if any, mining footprints. We would suggest that this model, or a similar one with a different diamond proportions, could be used in conjunction with one of the above models. The above models can identify draw points that should be in the initial cave. If the cave sequence is not well shaped, the 2Cone model could be used to suggest improvements. With an additional constraint specifying a set of draw points to open in the first period, the 2Cone model can be run to create well shaped caves over all periods. Consultation with the expert mine planner can suggest other diamond proportions that are acceptable. ColGen Model The ColGen model uses the Column Generation technique to separate the task of generating feasible caves from that of picking the cave in each period that maximizes the total NPV. The later task still utilizes Integer Programming techniques, but we were free to use other methods to generate the caves. In this work, we used an algorithm based on the 2Cone model to generate the feasible caves for each period. The technique first requires a set of caves for all periods, then the Master program selects the sequence that maximizes the total NPV. The shadow prices generated by the Master program are then used to generate new caves that will improve the NPV, and so the model iterates until no more improvement is possible. Creating the set of caves for each period proved complicated as there was not always sufficient information in the shadow prices to create reasonably sized caves that overlapped in consecutive periods. By reasonably sized we mean close to but not exceeding the Global Capacity Constraint. The algorithms added to estimate limits on cave sizes in each period, and generate feasible preceding and succeeding sequences greatly increased the size of the model, and the time to run. On the small trial data set, the result from the ColGen model was within half a percent of the 2Cone model result. This is not surprising as both models were based on the same cave shape definition. As with the 2Cone model, the larger trial data set proved to be too large, and a solution was not returned. Greedy Algorithm The Greedy model is the algorithm that was developed to generate the initial cave sequence for use in the ColGen model. The algorithm enumerates all 183 Chapter 3. Sequence Optimization in Block Cave Mining caves of maximum size, and chooses the one with the highest total of draw point weights for the first period cave. For each subsequent period, it counts how many draw points have completed, then enumerates all caves of the new size that fully overlap with the cave from the previous period. The cave with the highest incremental draw point weight becomes the cave for that period and the process is repeated. In this work, caves matching the shape of those in the 2Cone model were generated. This model ran to completion on both the trial data sets. On the smaller data set, the model generated a sequence that achieved $3,333,520 which is essentially the same as the ColGen result of $3,333,938. On this data set, the additional iterations of the ColGen model did not significantly improve on the starting sequence generated by the Greedy algorithm. As mentioned above, the Greedy algorithm / ColGen model result is also very close to that from the 2Cone model. On the larger data set, the Greedy algorithm ran for over 6 days, but did return a solution, unlike the ColGen and 2Cone models. The value of the solution was 94.5% of the Basic model solution, not leaving much room for improvement had the ColGen model run to completion. In these two trial data sets, the Greedy algorithm performed strongly. If cave shapes that can not be well defined in the integer programming framework are desired, we suggest running this model, modified to the desired shape, and comparing the result to the Basic model benchmark. If there is a large gap, the ColGen model may add some value, but if not, the extra computation may not be worthwhile. In theory, the Greedy algorithm (and the ColGen model) can be adapted to generate other cave shapes. The results from these trials show that even with the simple, single cave definition, the Greedy algorithm can be slow. This suggests additional, or more complex cave definitions may perform even worse. Depending on the size of the data footprint, and the complexity of the problem, the mine planner can decide if using this model is worth the time to reach a solution. Future Work This work is a first attempt to use models to generate opening sequences for planning block cave mining. In the process, we learned, as often happens with real world problems, that writing the rules, or constraints of this task is complex. We have shown that it is possible to write models that can generate sequences, but before any further work is done we recommend consulting expert mine planners to learn their opinion on how useful the models are. 184 Chapter 3. Sequence Optimization in Block Cave Mining The slow speed of the Greedy algorithm directly impacts the speed of the ColGen model. The Greedy algorithm is very simplistic and evaluates the cave made by each pair of draw points, approximately nDP T S 2 of them. Future work could be focused on this algorithm to make it more efficient and faster. Much further down the road, if useful sequences can be generated, work could be done to combine this program with PC-BC, Gemcom’s mine planning software. In real life, much of the data is stochastic (ore price, discount rate, etc.) and so having a model that takes this into account would be valuable. The models in this work are all in the integer programming setting. There are many other optimization techniques available, and this problem could be tried on them. For example constraint programming (Marriott 1998) or genetic algorithms and other global optimization algorithms (Weise 2008). 185 References Ahuja, Ravindra K., Thomas L. Magnanti, and James B. Orlin. 1993. Network flows: theory, algorithms, and applications. Prentice Hall. Americal Museum of Natural History. (2008). Mining a Kimberlite Pipe. In The Nature of Diamonds. Retrieved June 25, 2008 from http://www.amnh.org/exhibitions/diamonds/mining.html. Barber, J., Thomas, L., & Casten, T. 2000. Freeport Indonesia’s Deep Ore Zone Mine. Proceedings MassMin 2000. 291. Calder, K., Townsend, P., & Russell, F. 2000. The Palabora Underground Mine Project. Proceedings MassMin 2000. 219–225. Chen, F. 2000. Optimal Policies for Multi-Echelon Inventory Problems with Batch Ordering. Operations Research. 48 no. 3, 376–389. Dantzig, G. B., & Fulkerson, D. R. 1956. On the Max-Flow Min-Cut Theorem of Networks. Linear Inequalities and Related Systems 215–221. Diering, T. 2000. PC-BC: A Block Cave Design and Draw Control System. Proceedings MassMin 2000. 469–484. Diering, T. 2006. Personal communication. Garey, M.R., & Johnson, D.S. 1979. Computers and intractability: a guide to the theory of NP-completeness. W.S. Freeman.. GEMCOM. 2006. Company website. http://www.gemcomsoftware.com/ Gertsch, R. E. & Bullock, R.L. 1998. Techniques in Underground Mining. Society for Mining, Metallurgy, and Exploration (U.S.) Giannini, L. M., Caccetta, L., Kelsey, P., & Carras, S. 1991. PITOPTIM: A New High Speed Network Flow Technique for Optimum Pit Design Facilitating Rapid Sensitivity Analysis. AusIMM Proc. 2 57–62. Howard, R. 1960. Dynamic Programming and Markov Processes. (M.I.T. Press). Lee, J., Leung, J. & Margot, F. 2003. Min-up / Min-down Polytopes. Tepper School of Business Paper 272, http://repository.cmu.edu/tepper/272 186 Chapter 3. Sequence Optimization in Block Cave Mining Lerchs, H. & Grossman, I. F. 1965. Optimum Design for Open Pit Mines. CIM Bulletin. 58 (January) 47–54. Lübbecke, M & Desrosiers, J. 2005. Selected Topics in Column Generation. Operations Research [serial on the Internet]. 53 no. 6, 1007-1023. Malkin, P. & Wolsey, L. 2006. Possible Formulations of the Mining Problem. Manuscript. Malkin, P. & Wolsey, L. 2006. Personal Correspondence. Marriott, K. & Stuckey, P. J. 1988. Programming with Constraints: An Introduction. MIT Press, Cambridge MA Newman, A., Rubio, E., Weintraub, A. & Eurek, K. 2010. A Review of Operations Research in Mine Planning. Interfaces 40 no. 3, 222-245. Newman, P. 1987.“Indirect Utility Function” The New Palgrave:A Dictionary of Economics. Eds John Eatwell, Murray Milgate and Peter Newman. Palgrave Macmillan. 1987. Parkinson, A., McCormick, S.T. 2005. Optimal Replenishment with Two Delivery Sizes. (MSOM Conference Presentation). Pinedo, M. 1995. Scheduling: Theory, Algorithms, and Systems. Prentice Hall.. Puterman, M. 1994. Markov Decision Processes. (John Wiley and Sons). Rajan, D. & Takriti, S. 2005. Minimum Up/Down Polytopes of the Unit Commitment Problem with Start-Up Costs. IBM Research Report Resolution Coper Mining. (2009). Block Cave Mining. Block Cave Mining. Retrieved November 12, 2009, from http://www.resolutioncopper.com/res/ourapproach/BlockCaveMining.pdf Rubio, E. 2002. Long Term Planning of Block Cave Operations using Mathematical Programming Tools. Masters Thesis, UBC. Van Vyve, M. 2006. Personal communication. Veinott, A. 1965. The Optimal Inventory Policy for Batch Ordering. Operations Research. 13 1103–1145. Weintraub, A., Pereira, M. & Schultz, X. 2008. A Priori and A Posteriori Aggregation Procedures to Reduce Model Size in MIP Mine Planning Models. Electronic Notes Discrete Math 30 297- 302. Weise, T. 2008. Global Optimization Algorithms - Theory and Applications -. http://www.it-weise.de. Wolsey, L. 1998. Integer Programming. John Wiley and Sons.. 187 Chapter 3. Sequence Optimization in Block Cave Mining Yegulalp, T. M., & Arias, A. J. 1992. A fast algorithm to solve the ultimate pit limit problem. Proc. 23rd International APCOM Symposium 391– 397. 188
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Essays on sequence optimization in block cave mining...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Essays on sequence optimization in block cave mining and inventory policies with two delivery sizes Parkinson, Anita Frances 2012
pdf
Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
Page Metadata
Item Metadata
Title | Essays on sequence optimization in block cave mining and inventory policies with two delivery sizes |
Creator |
Parkinson, Anita Frances |
Publisher | University of British Columbia |
Date Issued | 2012 |
Description | Chapter 1 is an introductory chapter for this thesis work. It sets the scene by describing the motivation and industrial setting for each project. In Chapter 2, "Optimal Inventory Replenishment with Two Delivery Sizes", we consider a periodic review inventory system where a retailer can order in multiples of a fixed quantity Q₁, or multiples of Q₂ = 2Q₁, where the per unit material cost is less for ordering Q₂. We extend results of Veinott, and of Fangruo Chen, to show that an optimal replenishment policy has a reorder point R, as well as a second parameter controlling when the last order should be for Q₁ instead of Q₂ under a linear cost structure. In Chapter 3, "Sequence Optimization in Block Cave Mining" we investigate the use of integer programming models to aid the practitioner in the early planning stages of a Block Cave Mine. Given the footprint of the ore body divided into draw points or grid squares, sequence optimization determines which draw points to open in which period to meet the physical constraints of the mining process and maximize the total net present value of the mine. Traditionally done by trial and error by experts in the field, this is a first attempt to use modelling techniques to automate and optimize the process. We develop three integer programming models and discuss the challenges of formulating the problem in this framework. Two additional models are developed for comparison, one using the Column Generation technique and one using a greedy or myopic algorithm. All models are run on two data sets provided by our industrial partner, and the performance and results are compared. This work demonstrates that integer programming models can generate opening sequences but, like many "real life" problems, this one is complicated. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2012-08-21 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0073023 |
URI | http://hdl.handle.net/2429/42988 |
Degree |
Doctor of Philosophy - PhD |
Program |
Business Administration |
Affiliation |
Business, Sauder School of |
Degree Grantor | University of British Columbia |
GraduationDate | 2012-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2012_fall_anita_parkinson.pdf [ 12.81MB ]
- Metadata
- JSON: 24-1.0073023.json
- JSON-LD: 24-1.0073023-ld.json
- RDF/XML (Pretty): 24-1.0073023-rdf.xml
- RDF/JSON: 24-1.0073023-rdf.json
- Turtle: 24-1.0073023-turtle.txt
- N-Triples: 24-1.0073023-rdf-ntriples.txt
- Original Record: 24-1.0073023-source.json
- Full Text
- 24-1.0073023-fulltext.txt
- Citation
- 24-1.0073023.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
data-media="{[{embed.selectedMedia}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0073023/manifest