Close
About
FAQ
Home
Collections
Login
USC Login
Register
0
Selected
Invert selection
Deselect all
Deselect all
Click here to refresh results
Click here to refresh results
USC
/
Digital Library
/
University of Southern California Dissertations and Theses
/
Routing and inventory model for emergency response to minimize unmet demand
(USC Thesis Other)
Routing and inventory model for emergency response to minimize unmet demand
PDF
Download
Share
Open document
Flip pages
Contact Us
Contact Us
Copy asset link
Request this asset
Transcript (if available)
Content
ROUTING AND INVENTORY MODEL FOR EMERGENCY RESPONSE TO MINIMIZE UNMET DEMAND by Zhihong Shen __________________________________________________________________ A Dissertation Presented to the FACULTY OF THE GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulfillment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (INDUSTRIAL AND SYSTEMS ENGINEERING) August 2008 Copyright 2008 Zhihong Shen Acknowledgements A journey is easier to travel together. It is a great pleasure to thank the many people who accom- panied and supported me and made this thesis possible. Over the course of the last 4 years, I have been privileged and grateful to work closely with my Ph.D. advisors: Prof. Maged Dessouky and Prof. Fernando Ord´ o˜ nez. Without the guidance, advice, encouragement and high quality academic training they provided, this work is impossi- ble. Their profound knowledge, their enthusiasm and integral view on research, has made deep impression on me and shaped my perspective on research fundamentally. Besides of being ex- cellent advisors, they were as close as great supportive and respectful friends to me. I will never forget those discussions we made during the most difficult times during my Ph.D. study, which brought me so much comfort and laughter and helped me regain my strength and confidence to accomplish this research. I am indebted to my qualifying and thesis committee members: Prof. Sheldon Ross, Prof. Randy Hall and Prof. Sven Koenig for their detailed review, constructive criticism and excellent advice during different stages of this research. I warmly thank Prof. Mohamed I. Dessouky and our department chair Prof. James E. Moore, II, for their constant support and encouragement for my research. ii I owe my most sincere gratitude to the many people who have taught me and advised me during different stages of my life and without their inspiration, I could not reach the point I am now in the academic pursuing. They include, my high school teachers: Mr. Zhenjiang Wang, Mr. Chaozhi Tang; my college teachers: Prof. Daoheng Yu, Prof. Junxiong Tang; my graduate school teachers and advisors: Dr. Olivier Peyran, Prof. Kim-Chuan Toh; my advisor in Temasek Lab: Dr. Kai-Yew Lum and my supervisor in asprecise Pte. Ltd.: Dr. Qizhang Liu. I wish to thank my peer students in ISE department: Dr. Hongzhong Jia, Dr. Worawan Suteewang, Dr. Jie Liu, Mr. Pavan Murali, Ms. Hui Wang, Mr. Haifeng Ji, Mr. Indrajeet Dixit, Ms. Yen-Ming Lee, and Mr. Diego Prats, for sharing not only the many heat discussions on the academic problems but also lots of fun outside the office and classroom. I am very grateful to my best friends in Los Angeles: Dr. Bei Wang, Ms. Min Xie, Dr. Timmy Lee, Ms. Chen Wang and Mr. Hua Hui; and my best friends over the years: Ms. Yue Sun, Ms. Juanjuan Qi, Dr. Kai Li, Mr. Yangfan Zhou, Ms. Shanshan Chen, Ms. Danni Guo and Mr. Seng-Piow Loh. I am grateful for their help in getting me through the difficult times, their accompaniment when I was frustrated, and all the emotional support, entertainment and caring they provided. My debt to Mr. Hao Xu, is beyond words. He was the single person that I was at ease with in both good times and bad times. The understanding and encouragement, the fun and love he brought to me is unforgettable. I am grateful to him for a love that has grown and been wonderful. Last but not least, the special gratitude is due to my family: my mother, Ms. Hesheng Pu; my father, Mr. Yuanchen Shen; my brother; Mr. Zhiqian Shen; my sister-in-law, Ms. Jingyuan Li and my niece, Xinmeng Shen. Their unstopping love and trust has always been the most important thing to me in this world. My father passed away while I was doing the Ph.D. study. It was him iii who passed me the curiosity and desire to the knowledge and wisdom through the blood. This dissertation is dedicated to him and my mother. iv Table of Contents Acknowledgments ii List of Tables vii List of Figures viii Abstract x Chapter 1 Introduction 1 1.1 Perishable Inventory Management Problem 7 1.2 Vehicle Routing Problem 9 1.3 Contribution 10 Chapter 2 Literature Review 13 2.1 Perishable Inventory Management Problem 13 2.1.1 Policies 14 2.1.2 Ordering Policy 16 2.1.2.1 Assumptions 16 2.1.2.2 Modeling Methods/Solution Approaches 17 2.1.3 Gap 20 2.2 Vehicle Routing Problem 20 2.2.1 Stochastic Vehicle Routing Problem 21 2.2.1.1 Classification 22 2.2.1.2 SVRPs with Stochastic Customers and Demands 23 2.2.1.3 SVRPs with Stochastic Travel Time and Service Time 25 2.2.2 Routing with Profit 26 2.2.3 Gap 28 Chapter 3 Perishable Inventory Management Problem with Minimum Inventory Volume Constraint 29 3.1 Problem Statement 29 3.2 EMQ Model with Perishability 32 3.3 Modified EMQ Model 34 3.4 Constraint on the Inventory Cycle Length Parameter 40 3.5 Total Cost Evaluation and Boundary Conditions 44 3.5.1 Notation 44 3.6 Local and Global Optimality 64 3.6.1 Local Optimality 65 3.6.2 Global Optimality 65 3.7 Exact Solution Algorithm and Complexity Analysis 69 3.8 Computational Experiments 70 3.8.1 Parameter Estimation 70 3.8.2 Model Comparison 73 3.8.3 Sensitivity Analysis 75 v Chapter 4 Stochastic Routing Problem to Minimize Unmet Demand 81 4.1 Problem Statement 81 4.2 Two Stage Solution Process 83 4.3 A Priori Routing Problem 84 4.3.1 Notation 85 4.3.2 Deterministic Model 87 4.3.3 Stochastic Model 92 4.4 Recourse Strategies 95 4.4.1 Notation 95 4.4.2 LP Recourse 96 4.4.3 Knapsack Recourse 97 4.4.4 Reroute Recourse 100 4.5 Solution Approaches 100 4.5.1 Tabu Search for A Priori RoutingProblem 100 4.5.1.1 Object-Oriented Tabu Search Framework 101 4.5.1.2 Tabu Search Procedure 105 4.5.1.3 Neighborhood Structure and Tabu Moves 105 4.5.1.4 Tabu Tenure, Aspiration and Stopping Criteria 108 4.5.1.5 Post-processing for a Complete Route Construction 108 4.5.2 Approximation Heuristic for Knapsack Recourse 109 4.6 Computational Experiments 111 4.6.1 Numerical Setup and Simulation Processes 112 4.6.2 Evaluation of Tabu Search Quality 114 4.6.3 Evalution of A Priori Route 118 4.6.3.1 The Deterministic Routes vs the CCP Routes 118 4.6.3.2 The Role of the Safety Factor 120 4.6.4 Comparison of Recourse Strategies 123 Chapter 5 Conclusion and Future Work 125 5.1 Conclusion 125 5.2 Future Work 127 Bibliography 129 vi List Of Tables 3.1 Estimated numerical value of the parameters used in the experiment. . . . . . . . 72 4.1 Percentage of unmet demand over total demand for the deterministic model. . . . 115 4.2 Remaining supply at depot for the deterministic model. . . . . . . . . . . . . . . 115 4.3 Tabu Heuristic Result with CPLEX Bounding for 10 Customers and 3 Vehicles. . 117 4.4 Tabu Heuristic Result with CPLEX Bounding for 20 Customers and 4 Vehicles. . 117 4.5 Tabu Heuristic Result with CPLEX Bounding for 50 Customers and 10 Vehicles. 118 4.6 Average percentage of the unmet demand over the total demand for DP and CCP models (lognormal distribution) (%) . . . . . . . . . . . . . . . . . . . . . . . . 119 4.7 Unmet Demand Percentage Comparison Between 3 Recourse Strategies for 50 Customers and 10 Vehicles. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124 vii List Of Figures 1.1 Distribution Processes for Emergency Response. . . . . . . . . . . . . . . . . . 3 1.2 Problem Decomposition for Emergency Response. . . . . . . . . . . . . . . . . 4 3.1 Goal of the Perishable Inventory Management Research for Emergency Response. 31 3.2 Perishable Inventory System with I min = 0 or small I min . . . . . . . . . . . . . . . 33 3.3 Illustration for the Modified EMQ Model. . . . . . . . . . . . . . . . . . . . . . 35 3.4 All 5 Cases of Possible Scenarios. . . . . . . . . . . . . . . . . . . . . . . . . . 38 3.5 Classification Criteria for 5 cases. . . . . . . . . . . . . . . . . . . . . . . . . . 39 3.6 Illustration for the Stability Condition. . . . . . . . . . . . . . . . . . . . . . . . 40 3.7 Graph Illustration for an Extreme Boundary Case. . . . . . . . . . . . . . . . . . 41 3.8 Graph Illustration for Case 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 3.9 Graph Illustration for Case 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 3.10 Graph Illustration for Case 3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 3.11 Graph Illustration for the Boundary Value of Case 3. . . . . . . . . . . . . . . . 54 3.12 Graph Illustration for Case 4. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 3.13 Graph Illustration for the Boundary Value Case 4. . . . . . . . . . . . . . . . . . 58 3.14 Graph Illustration for Case 5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 3.15 Graph Illustration for Boundary Value of Case 5. . . . . . . . . . . . . . . . . . 63 viii 3.16 Graph Illustration for the total cost w.r.t. Q . . . . . . . . . . . . . . . . . . . . . 64 3.17 Graph Illustration for the Local Optimum Point . . . . . . . . . . . . . . . . . . 65 3.18 Graph Illustration for the Global Optimum Proof - General Case . . . . . . . . . 67 3.19 Graph Illustration for the Global Optimal Proof - Lower Bound. . . . . . . . . . 68 3.20 Graph Illustration for the Global Optimal Proof - Upper Bound. . . . . . . . . . 68 3.21 Complete Exact Algorithm for Solving the Modified EMQ Model . . . . . . . . 69 3.22 Total cost comparison of a standard model and our proposed model . . . . . . . . 73 3.23 Profit v.s. Y with Q max = 1600 at different I min level . . . . . . . . . . . . . . . . 76 3.24 Profit v.s. Y which allows salvage at most 300 million pills . . . . . . . . . . . . 77 3.25 Profit v.s. Y with Q max = 1600 at different unit government-paid price . . . . . . 78 3.26 Illustration of the different slopes in the profit v.s. Y plot . . . . . . . . . . . . . 79 4.1 A priori route model at the planning stage. . . . . . . . . . . . . . . . . . . . . . 86 4.2 The architecture of the Tabu search implementation. . . . . . . . . . . . . . . . . 102 4.3 The architecture of the optimization engine for tabu search. . . . . . . . . . . . . 104 4.4 Tabu Search Procedure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 4.5 Illustration of the DEM Moves. . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 4.6 Illustration of the λ-interchange Moves. . . . . . . . . . . . . . . . . . . . . . . 108 4.7 Flowchart of Simulation Process. . . . . . . . . . . . . . . . . . . . . . . . . . . 114 4.8 Simulation comparison for the DP and CCP models. . . . . . . . . . . . . . . . . 119 4.9 Simulation comparison for different safety factors. . . . . . . . . . . . . . . . . 122 ix Abstract Rapid and efficient wide-scale distribution of medical supplies plays a critical role in assuring the effectiveness in managing the risks of large-scale emergencies such as a bio-terrorism attack. Important issues in the design of such an efficient distribution network involve deciding how to route distribution vehicles and how to manage these inventories. The high uncertainty, limited supply and overwhelming demand associated with a large-scale emergency may result in signif- icant unmet demand, which is a direct representation of the most undesirable consequence - loss of life. Solving appropriate vehicle routing and inventory management problems in a coordinated manner can ensure the design of a logistic network capable of efficiently distributing medical sup- plies to decrease the potential fatalities in responding to a large-scale emergency. In this work, we develop models and solution approaches to solve a perishable inventory management problem and a vehicle routing problem in the context of emergency response to minimize unmet demand. To effectively manage the huge volume of perishable medical supplies to guarantee their freshness, we present a modified Economic Manufacturing Quantity (EMQ) model for perish- able items with a minimum inventory volume constraint. Minimizing the cost of maintaining such a system can be formulated as a non-convex non-differentiable unconstrained optimization problem. An exact algorithm with polynomial complexity is developed to solve this problem. To efficiently distribute the medical supplies for large-scale emergencies, a two-stage solution x approach is proposed by solving a stochastic routing problem in the first planning stage and a deterministic scheduling problem in the second operational stage. We formulate a mixed integer model which incorporates the routing with profits and the traditional complete routing for the first time to address the planning stage problem. A chance-constrained approach is applied to handle the uncertainty in both demand and travel time in the stochastic planning stage model. Three recourse strategies are implemented and compared for the operational stage. We develop a tabu heuristic and approximated knapsack heuristic to solve models in both stages. Numerical experiments are conducted to evaluate our models and solution approaches based on simulated large-scale emergencies. xi Chapter 1 Introduction Large-scale emergencies (or major emergencies) are defined as any event or occurrence that over- whelms local emergency responders, which severely impacts the operation of normal life, and has the potential to cause substantial casualties and property damage. Examples are natural dis- asters (earthquake, hurricane, flooding, etc.) and terrorist attacks, like September 11 th , 2001. To manage the risks and consequences of large-scale emergencies and mitigate their impacts to the population, wide-scale distribution plans of medical supplies must be developed. Careful and systematic pre-planning, as well as efficient and professional execution in re- sponding to a large-scale emergency, can save many lives. Rational policies and procedures ap- plied to emergency response could maximize the effectiveness of the scarce resources available in relation to the overwhelming demands. A key ingredient in an effective response to an emergency is the prompt availability of necessary supplies at emergency sites. It is challenging for practi- tioners to manage the huge volume of medical supplies to guarantee the readiness of the delivery and the freshness of the stock, as well as to deliver these massive supplies in a short time period to dispersed demand areas. Operations research models can play an important role in addressing and optimizing the logistical problems in this complex distribution process. Larson et al. (2005, 1 2006) conducted a detailed analysis based on well-known and recent large-scale emergencies. They emphasized the need for quantitative, model-oriented methods provided in the operations research field to evaluate and guide the operational strategies and actions in response to major emergencies. In certain emergency settings, medication or antidotes must be applied within a specified time limit from the occurrence of the event to maximize their effectiveness, and in some situations, to save lives. Traditional pharmaceutical supply chains are not adequate to push a huge volume of medical supplies to the affected population within a short time frame. Usually the local caches are only enough to serve the first responders. The majority of the medical supplies used would be from the national repository, which is called the SNS (Strategic National Stockpile) and is separated into pre-positioned push packages and a larger inventory managed by manufacturers (which is called Vender Managed Inventories - VMIs and it is about 80% of the total SNS in volume). The SNS is administrated by the Centers for Disease Control and Prevention in the United States. For example, to address emergencies of infectious diseases, the SNS contains 300 million doses of smallpox vaccines and enough antibiotic to treat 10 million people for anthrax. We illustrate the distribution process in Fig. 1.1. If there is a public health emergency severe enough to cause local supplies to run out, the SNS supplies can be accessed. Initial medical supplies from the push packages will be delivered to any state in the U.S. within 12 hours of federal and local authorization. This initial delivery can be followed by supplies from the VMIs if needed. After the stockpile reaches the regional central depot, it is each state’s responsibility to have plans to receive and distribute SNS medicine and medical supplies to local communities as quickly as possible. The state investigates the efficient strategies and best practices in identifying potential points of dispensing (PODs) sites where the population in the nearby neighborhood 2 comes to seek medical services, as well as in constructing suitable routes to distribute the supplies from the regional central depot to the selected PODs. SNS Regional Central Depot Selected Points of Dispensing Population Eligible Staging Area Push Packages Medical Services VMI (80%) Federal Stockpile Figure 1.1: Distribution Processes for Emergency Response. As the primary goal of this work is to establish an effective and rapid distribution system for the medical supplies to respond to a large-scale emergency, it is natural to decompose the whole problem into three separate however inter-connected key decisions, which can be integrated to provide a complete medical supply chain management solution to large scale emergencies (see Fig. 1.2). The medical supplies are the items of interests in the system. We are concerned with three key questions in the distribution process: 1. How to manage these medical supplies? 2. Where to store them? 3. How to distribute them? These three questions need to be addressed considering the specific challenges posed by large- scale emergency situations, such as: • a massive quantity of the medical supplies is required in the stockpile at all time; 3 • the medical supplies are perishable items; • medicines are needed to be delivered within a preferred time frame to maximize its effects; • the supply is limited within a short time frame compared with the overwhelming demand; • there is high uncertainty associated with the emergency scenarios, especially in the demand and the road condition (travel time). Where to store? (Coverage) How to distribute? (Coverage + Fast) How to manage? (Huge Volume + Perishable) Medical Supplies Facility Location Vehicle Routing Inventory Model Capacity Demand Pts Supplies Part II Part I Figure 1.2: Problem Decomposition for Emergency Response. A perishable inventory management model, a facility location model and a routing model are for- mulated to address these three questions separately to respond to the challenges presented above. The solutions to each problem can be coordinated into an overall efficient management strategy for the distribution process in the emergency response context. The minimum amount that must be kept in the inventory system is enforced by the SNS policy, which is the key constraint in the 4 inventory model and also provides the capacity constraint to the location model as well as the total supply constraint to the routing model. The results of the location model are the selected PODs, which will feed into the routing model as the demand points input. The facility location model and solution approaches are the topic of another Ph.D. thesis work (Jia, 2006). In this work, we are concerned with the inventory management problem faced by manufacturers to maintain a significant part of the SNS as VMIs (Part I of this work) and the distribution process after the SNS supplies have been received at the local level (Part II of this work). We strive to establish efficient inventory management policies and dispatching plans to help save lives through a timely delivery plan. In a large scale emergency setting, the initial pre-positioned supplies will typically not be suf- ficient to cover all demands in the recommended times. The unmet demand in an emergency situ- ation can directly result in loss of life, an impact that outweighs any operational costs. Therefore, the overriding objective in our work is to mitigate the impact of the emergency by minimizing the amount of unmet demand. It requires both an efficient inventory management scheme for the massive supplies to guarantee the readiness and freshness of the stockpile for continuous service and an effective plan to distribute the huge quantity of supplies as much as possible within the preferred time line. Currently, the integrated inventory-distribution problem has attracted more attention in the research community and has gained acceptance in many supply chain environments. This class of problems is “characterized by the simultaneous vehicle routing and inventory decisions that are present in a dynamic framework such that earlier decisions influence later decisions” (Baita et al. 1998). It routes supplies from a central depot to dispersed customer sites which have their own inventory systems to maintain. A constraint links the demand rate, the inventory level, and 5 the amount to be delivered by the routing plan at each customer site; and it is modeled to min- imize both routing and inventory carrying cost (Dror and Ball, 1987; Trudeau and Dror, 1992; Viswanathan and Mathur, 1997). Different groups of researchers extended it to the case of a multi-depot problem (Hwang, 1999), the production-distribution-inventory problem (Fumero and Vercellis, 1999), and the case of incorporating backorders (Abdemaguid et al., 2006), etc. However, these models in the current literature cannot directly translate to address our inventory- distribution problem to respond to large-scale emergencies in two aspects. First, in the spatial aspect, most problems in the literature assumes a small-scale independent inventory system run- ning at each distributed demand point; while our problem is concerned with the efficient inventory control plan for a huge volume of perishable items at a central manufacturer’s warehouse. Sec- ond, in the temporal aspect, the inventory decision and the routing decision are faced about at the same time and hence can be made simultaneously to further reduce costs in the traditional supply chain context; however, in our problem, the inventory management problem is the focus well before any emergencies have happened at the manufacturer’s warehouse and the routing plan is encountered both before and after a specific emergency at the regional level. The issues in both spatial and temporal aspects make it difficult to combine the inventory and routing problems into a single model in response to large-scale emergencies. Hence we do the decomposition and formulate separately as illustrated in Fig. 1.2. We formulate the problem of managing the scarce (relative to the overwhelming emergency demand) and huge volume (relative to the regular market demand) medical supplies as a modi- fied Economical Manufacturing Quantity (EMQ) model for perishable items, and the problem of disbursing supplies as a vehicle routing problem (VRP) to minimize unmet demand. These two problems and a facility location problem (Jia, 2006) are closely coupled and they are integrated 6 to provide a complete medical supply chain management solution to large-scale emergencies. The differences between the models developed in this research and the available models in the existing literature lie in that the developed models take into account the unique characteristics of large-scale emergencies (e.g., overwhelming demand, high degree of uncertainties of travel time and demand requirements, and low probability of occurrence) and hence could better mitigate the impact of the emergencies. 1.1 Perishable Inventory Management Problem The readiness and freshness of the continuous massive supply from the VMIs is crucial to the distribution process under the stress from both the massive volume requirement and the short delivery time preference. The primary goal of the vendor stockpile inventory management for the SNS is to ensure the availability of the huge volume of medical supplies at any point of time. At the same time, it is also imperative to achieve the maximal economical value of the stockpile inventory by exercising an inventory policy which addresses both the regular market demand and the SNS minimum inventory volume requirement. As we discussed before, the majority of the SNS are from the VMIs and firms are paid by the government for the production and storage of the drugs. For example, in a potential anthrax attack, the federal government aims to disperse an- tibiotics to 10 million people. This stockpile represents enough Cipro, a common antibiotic with a 9 year lifespan that works against anthrax and other infections, to meet regular market demand for years. Considering the obvious fact that it is more efficient to sell some of the stockpile rather than holding it until it expires, current SNS policy allows firms to sell the drugs to the private market (where Cipro works as a handy antibiotic for non-anthrax related illnesses) when they are 7 6 months prior to expiration rather than let the drugs spoil. However, the size of the stockpile is large compared with the expected market demand which leads to waste. Given more time and flexibility there is the opportunity to capture a significant amount of salvage value for the unsold stockpile. From the manufacture’s perspective, if it can apply a more sophisticated inventory holding policy which allows for the constant usage of the stockpile to meet the regular market demand and refill with new production at the same time to maintain the minimum stockpile re- quirement, then the firm can save on the total cost in maintaining the stockpile inventory, hence making it possible to further reduce the price charged to the government. The two main challenges for this inventory management problem are the massive minimum inventory requirement and the perishable property of the medical supplies. In the large-scale emergency context, the minimum inventory requirement (which is the medical supplies required for the SNS) is huge and in a scale which is comparable with the total regular market consumption within its shelf life. Any trivial extension of current inventory models is no longer valid since the minimum inventory stockpile cannot be timely consumed by the regular demand without significant spoilage. A more sophisticated inventory ordering policy is needed to guarantee the freshness and readiness of the required minimum volume in stock incorporating with the regular market demand as well as to reduce the maintenance cost of such an inventory system. We propose a modified EMQ model with a minimum inventory volume constraint for perishable items to address these issues. With the aid of inventory plots, we can formulate the problem of minimizing the total cost of maintaining the proposed inventory system as a non-convex non- differentiable unconstrained optimization problem with respect to the production batch size (Q). We propose an exact solution procedure to obtain the optimal production batch size (Q ∗ ) with polynomial complexity. Numerical experiments are conducted to demonstrate the advantage of 8 our proposed model over a standard model and a sensitivity analysis is performed to illustrate the tradeoffs between different government controlled parameters in the system. 1.2 Vehicle Routing Problem To respond to a large scale emergency, a vehicle routing/dispatching model provides a method to plan for the first responder’s reaction to the emergency, and in the event an emergency occurs, a methodology to quickly adjust the dispatching policy to provide an efficient response. In a large scale emergency, we consider that the objective of the routing problem is to minimize the unmet demand that can occur because resources and first responders are overwhelmed. Unmet demand is an overriding objective in an emergency situation since it can result in loss of life, which outweighs traditional VRP objectives. An additional important aspect in meeting the demand in a large-scale emergency is the response time. Since the administration of the medication to the population within a preferred time-frame makes an appreciable health difference in a large-scale emergency, it is important to associate each demand with a time window, to address situations in which a late delivery directly leads to loss of life, such as in an anthrax attack. We decompose our problem into a two-stage problem. In the first planning stage, we formu- late the problem into a mixed integer programming (MIP) model well in advance of any possible emergency to generate pre-planned routes, which can be used for mock trial runs and training pur- poses. The planning stage model combines both the routing model with profits and the traditional complete routing model for the first time. We represent the randomness in the vehicle routing problem such as uncertain demand requirements and travel times assuming a known probability distribution. A chance-constrained model is applied to handle the uncertainty which guarantees 9 that the time and demand constraints are not violated more than a prescribed probability. More- over, because of the massive service requirements, the demand at a location is not necessarily satisfied by a single truckload. As such, the planning stage routing problem allows for split de- livery (i.e., a point can be visited more than once if the demand exceeds the load capacity of available vehicles) to the demand points from multiple vehicles. In the second operational stage, at the time of the emergency when more information is revealed, we need to quickly respond to the event and generate the delivery requirements with the planned routes as well as making ad- justments on the routes if necessary. We refer to the second stage model as the recourse strategy. We propose three different recourse strategies: adjust the quantity in each vehicle (LP recourse strategy); adjust quantity and skip customers (Knapsack recourse strategy) and complete reroute (re-planning strategy). We propose a tabu heuristic for the planning stage model and an approx- imation heuristic for the knapsack recourse strategy. Numerical experiments are conducted to illustrate the effectiveness and efficiency of the models and the proposed solution approaches. 1.3 Contribution In this work, we investigate and analyze the special requirements for the distribution process in responding to a large-scale emergency. To better mitigate the risks and consequences of large- scale emergencies, the process has been decomposed into three closely inter-connected problems. In this thesis, we focus on two of them: a perishable inventory management problem and a routing problem, which can be solved in a coordinated manner to minimize the unmet demand, which is, to decrease the potential fatalities in a large-scale emergency. The contribution of this thesis are: 10 • propose a modified EMQ model for perishable items with a minimum inventory volume constraint (first time in the literature); • demonstrate the problem of minimizing the total cost of maintaining the proposed inventory system as an unconstrained non-convex non-differentiable optimization problem; • present an exact solution procedure with polynomial complexity to solve the modified EMQ model; • estimate parameters in the proposed EMQ model for a potential anthrax scenario and con- duct numerical experiments to illustrate the advantage of our model over a standard model; • perform sensitivity analysis on government controlled parameters in the proposed EMQ model and provide insights from perspectives of both parties (the federal government and the firms); • separate the distribution of medical supplies from a regional central depot into a stochastic planning stage problem and a deterministic operational stage problem; • formulate a mixed integer model to first time combine the routing with profits and the traditional vehicle routing model for the planning stage routing problem; • apply the chance-constrained technique to handle the stochastic elements in both demand and travel time in the planning stage routing problem; • propose three recourse strategies for the operational stage scheduling problem; • develop a tabu search heuristic and approximated knapsack heuristic for models in both stages; 11 • conduct numerical experiments to demonstrate the effectiveness and efficiency of the rout- ing models and the heuristics in both stages; • use simulation results to study the relationship between the deterministic model, the chance- constrained model and the robust optimization model for the planning stage routing prob- lem. The rest of this thesis is organized as follows: section 2 reviews the literature on the perishable inventory management problem and the routing problem. In section 3, we state the requirement and assumptions imposed by the massive perishable inventory management problem and formally present the model followed by an exact solution procedure, its complexity analysis and numerical experiments. In section 4, we analyze the characteristics of large-scale emergency and its special requirements on the routing problem; then we propose a two stage solution process with detailed models and solution approaches on both stages; numerical experiments and analysis on the re- sults conclude this section. In section 5, we finish with conclusions and propose future research directions. 12 Chapter 2 Literature Review In this section, we review well-known problems in the literature which are relevant to our problem at hand. This section is organized as follows: we first review the perishable inventory manage- ment problem; then the routing problem under two broad categories: classical stochastic vehicle routing problems (SVRPs) and routing with profits, are reviewed. 2.1 Perishable Inventory Management Problem In most inventory systems, it is assumed that stock items can be stored indefinitely to meet future demands. However, the effects of perishability cannot be ignored for certain types of inventories, which undergo changes in storage so that they may become partially or entirely unsuitable for consumption as time passes. Typical examples are fresh produce, blood cells, chemicals, photo- graphic films, drugs and other pharmaceuticals. There are two types of perishable items: fixed lifetime and random lifetime. The fixed lifetime category includes those items whose lifetime is known a priori and is independent of all other parameters of the system. On the other hand, some product lifetime is a random variable with some specified probability distribution (a typical case is 13 exponential decay). Some items experience “age-independent” perishability, which may be phys- ically fit for use but have lost economic value. Examples include fashion garments, Christmas trees/cards. For those items which experience some physical decay are called “age-dependent” perishable goods. In some cases, the rate of deterioration or decay is low and negligible; its effect can be ignored so that such items are still treated as fixed lifetime goods in the literature, e.g. blood cells, drugs, etc. However, in many other cases, the rate of deterioration plays a major role and its impact must be explicitly considered. It can be further classified into constant decay rate and probabilistic decay rate (random lifetime). Raafat (1991) provided a thorough review on the continuously deteriorating inventory models. 2.1.1 Policies For the perishable inventory management problem, policies in four different aspects have at- tracted attention from the research community. • Ordering Policy: This is the most investigated policy in the literature for the perishable inventory problem. It focuses on answering the following two questions: when to order and how much to order. Nahmias (1982) provided a well known review on this topic. Based on different assumptions on the system setup and parameters, researchers modeled and solved this problem with many different approaches. We will discuss it in more detail in the next subsection. • Issuing Policy:Issuing policy concerns the sequence to remove items from a stockpile of finitely many units of varying ages. The most general approaches include issuing either the oldest units first (FIFO) or the newest units first (LIFO). Pierskalla and Roach (1972) 14 studied the optimal issuing policies for perishable goods with random supply and demand and fixed life-time under several possible objective functions. In most literature, the FIFO policy is widely adopted. • Disposal Policy: At the occasion when the items in stock move very slow, an overstock situation may happen that the potential sales value of excess stock minus the expected inventory maintenance cost does not match the salvage value; the strategic disposal of part of the inventory is desirable. The topic on when and how much to dispose as well as the optimum number of items to keep in the inventory under stochastic demand and perishing has been studied by Rosenfield (1989, 1992). In Vaughan’s (1994) work, he treated the age at which the vendor outdates units as a decision variable and analyzed the interaction between the ordering and outdating policies. • Pricing Policy: The pricing of a product is a well-established research topic in revenue management and accepted wide attention and interest among the research community. Here we are only interested in the pricing policy which is closely coupled with the ordering pol- icy in the perishable inventory management content (multi-period newsvendor problem). The price is a decision variable and the potential/forecasting demand is price-sensitive. The pricing decision and ordering quantity decision can be made either sequentially or simulta- neously by different modeling approaches (Gallego and Van Ryzin, 1994; Abad, 1996 and 2003; Burnetas and Smith, 2000; Chun, 2003). 15 2.1.2 Ordering Policy The literature on the ordering policy for the perishable inventory management has flourished over the past several decades. Research has been done with different assumptions and modeling methods/solution techniques. We will discuss the models based on the broad classification with the assumptions and the modeling methods/solution approaches. 2.1.2.1 Assumptions 1. Lifetime: As we discussed earlier in this section, the perishability property assumes either items are not suitable for use after a fixed amount of time (Goh et al., 1993; Ravichandran, 1995; Willian and Patuwo, 1999, 2004; Hsu and Lowe, 2001) or the usability continuously deteriorate over time (Giri and Chaudhuri, 1998; Yan and Cheng, 1998). 2. Periodic vs. continuous review: As a periodic review policy is easy and cheap to adopt in practice, it is the dominate methodology studied in the perishable inventory management liter- ature (Fries, 1975; Nandakumar and Morton, 1993; Hsu, 2000; Hsu and Lowe, 2001). However, continuous review policy attracts the attention among different groups of researchers (Weiss, 1980; Ravichandran, 1995; Liu and Lian, 1999) because of its amenability to theoretical analysis on the optimality condition for both Poisson and general demand renewal processes. 3. Demand process: Most of the perishable inventory literature treats demand as random variables. The early works assumed the demand in different periods to be independent and identi- cally distributed and a general probability distribution function is used to study the optimal policy (Fries, 1975; Chazan and Gal, 1977; Nandakumar, 1993). Later on, researchers usually assumed Poisson distribution on the arrival rate and demand size in order to use Markovian model to for- mulate the problem (Weiss, 1980; Graves, 1982; Ravichandran, 1995; Tekin et al., 2001). Goh 16 et al. (1993) and Kopach et al. (2006) assumed multiple independent demand level based on the age and emergency level and addressed the effects of different substitution policies for the blood bank application. Sobel and Zhang (2001) considered the combination of both deterministic and stochastic demand. 4. Purchase vs. production: The delivery rate of products into stock can be assumed to be instantaneous or not (infinite or finite production rate). Graves (1982) assumed a continuous production system and aimed to decide the optimal production rate. On the other hand, Yan and Cheng (1998) made the assumption that the production rate is given and made the decision on when to stop and restart the production. 5. Lead time: For the purchase case, zero lead time is assumed by default. Ravichandran (1995) relaxed this assumption and Williams and Patumo (1999, 2004) specifically studied how the positive lead time would affect the optimal ordering quantity level. 6. Lost sale and backorder: It is common to allow either lost sale or backorder in most systems. Yan and Cheng (1998) addressed the situation that some of the demand is lost and the rest is backlogged. 7. Holding cost: Traditional inventory models assume a constant holding cost per unit per unit time. Giri and Chaudhuri (1998) and Hsu (2000) treated the holding cost as a nonlinear function of the age of items in stock and the on-hand inventory level to address the effects of continuous deterioration. Ferguson (2007) studied how to estimate the holding cost parameters. 2.1.2.2 Modeling Methods/Solution Approaches: 1. Zero first order derivative point (stationary point) over total cost function: Based on different sets of assumptions, it is usually straightforward to write the governing equation on 17 the inventory level over time and then obtain the inventory carrying cost, along with the fixed ordering cost, the viable ordering cost, the shortage cost, and the outdating cost. This can be modeled as an unconstrained nonlinear system. If the total relevant cost is a continuous and second-order differentiable function over the decision variable of current system (in most cases, the ordering quantity or reorder level), then we can obtain the first-order stationary point as the optimum. (Ravichandran, 1995; Giri and Chaudhuri, 1998; Liu and Lian, 1999). 2. Heuristics/Approximations: For stochastic demand circumstances, exact optimal policies are not only difficult to compute, but also demanding to implement due to the requirement to keep track of the age distribution of the stock. Nahmias (1982) provided a good summary on the early works of approximated optimal policies and heuristics to obtain the optimal policy parameters. Nandakumar and Morton (1993) developed heuristics from “near myopic” bounds for the fixed-life perishable items and demonstrated the accuracy of these heuristics with less than 0.1% average error over a wide range of problems. Goh et al. (1993) applied different approximation methods to study a two-stage perishable inventory models which have items either “young” or “old” and compared two managerial policies that whether or not allow the substitution with “young” items when the “old” ones are out of stock. 3. Markovian model: The queueing model with impatient customers has been used to analo- gize perishable inventory systems. The queue corresponds to the inventory stockpile, service process to the demand, arrival of customers to the replenishment of inventory, and the time a customer will stay in queue before leaving due to impatience corresponds to the shelf-life. Early works (Chazan and Gal, 1977; Graves, 1982) usually wrote the descriptive transition probabil- ities then obtained the stationary distribution to evaluate the performance measures such as the 18 expected outdating amount. Weiss (1980) and Liu and Lian (1999) also use the performance measures to construct a cost function from which the optimal policy parameter can be computed. 4. Dynamic programming (DP): Since dynamic economic lot size problem was first pro- posed by Wagner and Whitin (1958), which reviews the inventory periodically and the demand is deterministic in every review period, dynamic programming technique has been widely adopted in solving many variations and extensions by different groups of researchers in the non-perishable inventory context. Hsu (2000) provided a brief review on this topic. There is a well-known zero- inventory property under which no inventory is carried into a production period that is necessary for an optimal solution. Hsu (2000) demonstrated that this zero-inventory property may not hold for any optimal solution with perishable items and proposed a new dynamic programming re- cursion based on an interval division property that solved the problem in polynomial time. Hsu and Lowe (2001) further extended the dynamic programming approach to solve the perishable economic lot size problem in which backorder is allowed. Sobel and Zhang (2001) used the DP technique to address a perishable inventory system with both stochastic and deterministic demand co-existing at the same time. 5. Fuzzy theory: Recently, rapid growth of fuzzy set theory brings it into many different application fields. Katagiri and Ishii (2002) introduced it into the perishable inventory control model with a fuzzy shortage cost and a fuzzy outdating cost; hence the expected profit function is represented with a fuzzy set. The effect of the fuzziness on the obtained ordering quantity is investigated. 19 2.1.3 Gap The perishable inventory management literature is abundant in discussing policies, modeling techniques and solution approaches based on different sets of assumptions which strive to make reasonable and reliable approximation of real-life problems. To the best of our knowledge, no work has been conducted in the literature with the minimum volume constraint on the inven- tory size throughout the planning horizon. This extension is trivial when the required minimum volume is not significant compared with the amount consumed by a regular market demand rate within the shelf life as it can be timely and completely refreshed by the regular demand. How- ever, when the minimum volume is huge and in a scale which is comparable with the total regular market consumption within the shelf life, such as the medical supplies required for the SNS in the large-scale emergency context, strategic inventory ordering and disposing policy is needed to guarantee the freshness and readiness of the required minimum volume in stock at all time as well as the low cost of maintaining such an inventory system. In this work, we address this modeling gap by formulating a perishable inventory management model with minimum volume inventory constraint and providing an exact solution approach to this model in chapter 3. 2.2 Vehicle Routing Problem The Vehicle Routing Problem is defined on a given topological graph. The classical standard VRP generates a set of routes which visit each customer exactly once. It aims to minimize the total travel time and/or the operational cost. When some uncertainty in the parameters is introduced into the problem, it is called a Stochastic VRP (SVRP). A comprehensive overview of the Vehicle Routing Problem and its stochastic variant can be found in Toth and Vigo (2002) and 20 Gendreau et al. (1996). Routing problems with profits are a generalization of the vehicle routing problems, where it is not necessary to visit all customers. A profit is associated with each demand node and a uniform deadline is applied to every node. The problem aims to maximize the total collected profits subject to the time restriction (constrained by the travel cost). This is called the Orienteering Problem (OP) when only a single vehicle is present and is motivated by an outdoor sport. The multi-vehicle version of the OP is referred to as the Team Orienteering Problem (TOP) in the literature. In comparison to the classical routing problems, the routing problems with profits have an added element of complexity. They require selecting a subset of nodes to visit while determining an order on the tour. If all the customers can be visited within the time constraints, then it becomes the classical routing problem and the extra complexity vanishes. We review these two well-known bodies of literature in the following subsections. 2.2.1 Stochastic Vehicle Routing Problem The Vehicle Routing Problem is defined on a given graph G=(V,A), where V ={v 1 ,v 2 ,...,v n } is a set of vertices andA ⊆ (v i ,v j ) : i6= j,v i ,v j ∈ V is the arc set. An optimal set of routes, composed of a cyclic linkage of arcs starting and ending at the depot, is selected to serve a given set of customers at vertices. This problem was first introduced by Dantzig and Ramser in 1959 to solve a real-world application concerning the delivery of gasoline to service stations. Toth and Vigo (2002) provides a comprehensive discussion on problem formulations, solution techniques, important variants and applications. Other general surveys on the deterministic VRP also can be found in Laporte (1992). The Stochastic VRP (SVRP) introduces some element of uncertainty in the parameters within the system in question. The SVRP differs from its deterministic counter- part: several fundamental properties of the deterministic VRPs no longer hold in the stochastic 21 case, and solution methodologies are considerably more complicated. A general review on the SVRP appeared in Gendreau et al. (1996). In this subsection, the general SVRP problem is briefly introduced and classified according to type of uncertainty, modeling methods and solution techniques. Two specific topics, SVRP with stochastic customers and demands and SVRP with stochastic travel time and service time, follow the general discussion. 2.2.1.1 Classification The Stochastic Vehicle Routing Problem can be broadly classified based on the following criteria: • Uncertainty in the problem: The uncertainty can be present in different parts of the Ve- hicle Routing Problem. It can be divided into VRP with stochastic customers (VRPSC), VRP with stochastic demands (VRPSD), VRP with stochastic travel time (VRPSTT) and VRP with stochastic service time (VRPSST). • Modeling method: The modeling method can also be the criterion to classify SVRPs. Stochastic VRPs can be cast into a stochastic programming framework, which is further di- vided into chance constrained program (CCP) and stochastic program with recourse (SPR). Another approach to model the SVRP is to view it as a Markov decision process. Some more recent modeling approaches include neurodynamic programming/reinforcement learn- ing methodology. • Solution techniques: Different solution techniques are the direct result of different mod- eling methods and the nature of the model. They usually broadly fall into two categories: exact methods and heuristic methods. The exact solution methods successfully applied to 22 the SVRP include branch and bound, branch and cut, integer L-shape method and gener- alized dynamic programming. Under some mild assumptions, several classes of chance- constrained SVRPs can be transformed into equivalent deterministic VRPs (Stewart and Golden, 1983; Laporte et al., 1989). Numerous heuristics have been used to solve the SVRP, e.g. the savings algorithm, tabu search, etc. The heuristic solution techniques usually can be further divided as constructive heuristics, improvement heuristics and meta-heuristics. 2.2.1.2 SVRPs with Stochastic Customers and Demands The VRP has stochastic demands (VRPSD) when the demands at the individual delivery (pickup) locations behave as random variables. The first proposed algorithm for the VRPSD by Tillman (1969) was based on Clarke and Wright (1964) savings algorithm. Another early major contribu- tion on the VRPSD comes from Stewart and Golden (1983), who applied the chance constrained programming and recourse methods to model the problem. Later on, Dror (1986) illustrated the impact that the direction of a designed route can have on the expected cost. A major contribu- tion to the study of the VRPSD comes from Bertsimas (1988, 1992). This work illustrated the a priori method with different recourse policies (re-optimization is allowed) to solve the VRPSD and derived several bounds, asymptotic results and other theoretical properties. A number of greedy heuristics were proposed and analyzed in an asymptotic sense. Bertsimas and Simchi- Levi (1996) surveyed the development in the VRPSD with an emphasis on the insights gained and on the algorithms proposed. A family of stochastic location-routing problems was studied by Laporte (1989). In this problem, the set of depots is not known a priori. Instead, depot sites with given operating costs must be determined from a candidate set, simultaneously with the optimal delivery routes through a set of customers having random supplies. The chance constrained model 23 was proposed and a branch and bound algorithm was applied. Besides the conventional stochastic programming framework, a Markov decision process for single stage and multistage stochastic models were introduced to investigate the VRPSD in Dror et al. (1989) and Dror (1993). More recently, a re-optimization type routing policy and a rollout algorithm for the VRPSD were intro- duced by Secomandi (2001). The VRP with stochastic customers (VRPSC), in which customers with deterministic de- mands and a probability p i of being present, and the VRP with stochastic customers and demands (VRPSCD), which combines the VRPSC and VRPSD, first appeared in the literature of J´ ez´ equel (1985), Jaillet (1987) and Jaillet and Odoni (1988). Bertsimas (1988) gave a more systematic analysis and presented several properties, bounds and heuristics. Unlike most of the work on the VRPSC which assume binary demands on each customer, Waters (1989), considered non-binary demands and compared different operating policies. Gendreau et al. (1995, 1996) proposed the first exact solution, L-shaped method, and a meta-heuristic, tabu search, for the VRPSCD. Another research direction of the VRP with stochastic demands or customers incorporates a dynamically changing environment, where the demands vary over time. This is motivated by some real world applications (e.g., mobile repair) and it aims to minimize waiting time, not travel cost. Bertsimas and Van Ryzin (1991, 1993) pioneered this work and named it dynamic traveling repairman problem (DTRP). They applied queuing theory to solve it and analyzed several dis- patching policies for both light and heavy traffic. Based on their results, Papastavrou (1996) and Swihart and Papastavrou (1999) defined and analyzed a new routing policy using a branching pro- cess. The multiple scenario approach (MSA) (Bent and Hentenryck, 2004) and waiting strategies (Branke et al., 2005) have also been introduced for this problem recently. 24 2.2.1.3 SVRPs with Stochastic Travel Time and Service Time Compared with stochastic customers and demands, the research on the stochastic travel time and service time problem of the TSP and VRP have received less attention. The VRP with stochas- tic travel time (VRPSTT) describes the uncertain environment of the road traffic condition. Kao (1978) first proposed heuristics based on dynamic programming and implicit enumeration for the TSP with stochastic travel time (TSPSTT). Carraway, Morin and Moskowitz (1989) used a gener- alized dynamic programming methodology to solve the TSPSTT. Laporte et al. (1992) performed systematic research on the VRP with stochastic service and travel time (VRPSSTT). They pro- posed three models for the VRPSSTT: chance constrained model, 3-index recourse model, and 2-index recourse model. They presented a general branch-and-cut algorithm for all 3 models. The VRPSSTT model was applied to a banking context and an adaptation of the savings algorithm was used in the work of Lambert et al. (1993). Jula et al. (2005) developed a procedure to estimate the arrival time to the nodes in the presence of hard time windows. In addition, they used these estimates embedded in a dynamic programming algorithm to determine the optimal routes. Hadjiconstantinou and Roberts (2002) formulated a VRP with stochastic service times to model a company who receives calls (failure reports) from customers and dispatches technicians to customer sites. They used a two stage recourse model and a paired tree search algorithm to solve it. 25 2.2.2 Routing with Profit Hayes and Norman (1984) first introduced the orienteering problem with real data for the 1981 Karrimor Mountain Marathon and proposed a dynamic programming solution approach. Tsili- girides (1984) deduced a common mathematical model, proposed heuristics and ran the numer- ical tests on up to 33 locations. Later on, a number of heuristics have been proposed by various groups of researchers, e.g., Golden et al. (1987, 1988), Ramesh et al. (1991), and Chao et al. (1996a). All these heuristics were tested using the small size benchmark problems provided in Tsiligirides (1984) with improving solution quality. More recently, many types of sophisticated meta-heuristic techniques have been applied on the OP, for example: artificial neural networks (Wang et al., 1995), tabu search (Gendreau et al., 1996), genetic algorithms (Tasgetiren and Smith, 2000), and ant colony optimization (Liang and Smith, 2001). These meta-heuristics can provide comparable results on much larger scale problems (up to 300 nodes); among them, the tabu search has been shown to provide an optimality gap that is less than 1% in a number of experiments. An- other focus of the OP research are exact solution methods. Laporte and Martello (1990) used a branch-and-bound approach with a knapsack bound; while Ramesh et al. (1992) applied a short- est spanning tree relaxation style bounding. Leifer and Rosenwein (1994) discussed an LP-based bounding procedure. Later on, Fischetti et al. (1997, 1998) and Gendreau et al. (1998a) proposed a branch-and-cut method to solve instances with up to 500 nodes. Feillet et al. (2005) provided a more detailed review on this topic. Two works have considered stochastic elements in the OP. Tang and Miller-Hooks (2005a) ad- dressed an orienteering problem with stochastic service time. The problem was formulated as a chance-constrained stochastic programming problem. An exact branch-and-cut algorithm and 26 a “construct-and-adjust” heuristic were proposed to solve it. Ilhan et al. (2008) focused on the uncertainty in demands in the OP and maximized the probability of collecting more than a pre-defined target level of the total return. They proposed an exact solution method based on a parametric formulation and developed a bi-objective GA heuristic. The routing problem with profits was extended to a multi-vehicle context by Butt and Cav- alier (1994). They were concerned with the athletes recruitment problem and proposed a greedy algorithm for solving it. Dunn (1992) applied dynamic programming to obtain a quick and useful schedule for a military operational situation. Chao et al. (1996b) used a straightforward extension of their heuristic for the OP (Chao et al., 1996a) to the TOP context and solved up to 102 nodes. Later on, Butt and Ryan (1999) proposed a branch-and-price exact solution procedure based on a set-partitioning formulation, and made efficient use of constraint branching and column gener- ation. In recent years, Tang and Miller-Hooks (2005b) introduced tabu heuristics for the TOP; Archetti et al. (2007) proposed two variants of a generalized tabu search algorithm and a variable neighborhood search algorithm and their computational experiments on benchmark problems (Chao et al., 1996b) generated best known solutions so far. All the works mentioned above for the TOP consider the constraint on the temporal dimension only and assume an uncapacitated fleet. Goel and Gruhn (2005) studied a more complicated problem arising in air-cargo transport. They were concerned with a multiple pickup and delivery problem with a heterogeneous fleet. Load acceptance, compatibility constraints, time windows constraints, capacity constraints, as well as a timely manner of handling dynamic requests were under consideration at the same time. A variable neighborhood search method (Mladenovic and Hansen, 1997) was developed to solve the problem on real-life scale test cases. 27 2.2.3 Gap Both the SVRPs and the routing with profits are well investigated respectively by many researches over the past several decades. However, to the best of our knowledge, there is no work to combine these two well-known problems together and model them in a single problem in current literature. In the distribution process of a large-scale emergency, a tight deadline is imposed to guarantee the effectiveness of the medicines to be delivered, which mimics the routing problem with profits and strives to minimize the potential fatality rate; on the other hand, a complete route is desirable since it provides a good starting point to guide a efficient recourse plan and is also considered as an effective strategy to handle the significant uncertainty associated with the emergency response. Our work fills this modeling gap and proposes a heuristic solution approach to this model in chapter 4. 28 Chapter 3 Perishable Inventory Management Problem with Minimum Inventory Volume Constraint 3.1 Problem Statement As part of the national emergency preparedness plan, the federal government keeps a Strategic National Stockpile (SNS), which includes a combination of vaccines, prophylaxis, and medical supplies to respond to a bio-terrorist attack on the United States. The SNS has stockpiled enough medicine to protect people in several large cities at the same time. About 20% of the SNS is ready to deploy as push packages and another 80% is in the form of VMIs. The federal government pays firms for producing and keeping inventory of these huge volume of stockpiles. The availability of the medicine within a short time line and the freshness of the massive volume stockpile places a unique challenge for the inventory management system for emergency response. Note that the federal government requires a minimum inventory level of the VMIs stockpile and these drugs have a limited shelf-life. The traditional EMQ (Economic Manufactur- ing Quantity) model can be readily extended to address the perishability property of the stockpile 29 with zero minimum inventory requirement by properly upper-bounding the EMQ cycle to guar- antee the inventory to be consumed within its shelf-life. In the case that the minimum inventory required in the stockpile is significantly smaller than the total regular market consumption within the shelf-life time of the drug, a tighter limit on the EMQ cycle which incorporates the considera- tion on the freshness of the minimum inventory as well will suffice. However, when the minimum inventory is in a scale comparable with the total regular market consumption within the shelf-life, a trivial extension to any current perishable inventory policies is no longer adequate. It is there- fore imperative to develop a new inventory policy specially geared to the perishable VMI system for the SNS, which not only satisfies the minimum inventory requirement but also minimizes the operational cost of maintaining such a system by incorporating the regular market demand. For example, in a potential anthrax attack, the stockpile contains enough antidote for 10 mil- lion people. This stockpile represents enough Cipro, a common antibiotic with a 9 year lifespan that works against anthrax and other infections, to meet regular market demand for years. Cur- rently, the SNS policy allows manufacturers to sell the pills 6 months prior to expiration rather than let the drugs spoil; however, considering the fact that the size of the stockpile is huge com- pared with the regular market demand while the drugs are so close to their expiration date, the potential salvage value is low. From the manufacture’s perspective, if it can apply a more so- phisticated inventory holding policy which allows the constant usage of the stockpile to meet the regular market demand and refill with new production at the same time to maintain the minimum stockpile requirement, then the firm will save on the total cost in maintaining the stockpile in- ventory, hence making it possible to further reduce the price charged to the government. From the government’s perspective, if it allows firms to sell the pills earlier than 6 months prior to the expiration date, there is an opportunity to capture a significant amount of salvage value for the 30 unsold stockpile. The unique challenge of this problem lies in efficiently maintaining a constant minimum level of perishable inventory. If items are non-perishable (or the minimum stockpile can be consumed by the regular market demand well before its expiration), we can apply a FIFO issuing policy which keeps a regular EMQ model on top of the minimum inventory level and sat- isfies both the extra minimum stockpile constraint and the regular market demand. However, the issue arises when the minimum stockpile is large and cannot be consumed by the market demand before it expires. Hence, in this research, we aim to propose a single inventory system which combines two types of demand: the regular market demand and the emergency demand (the per- ishable VMI for SNS) as shown in Fig. 3.1 to meet both of them and minimize the operational cost from the manufacturer’s perspective. VMI (Imin) (Emergency Demand) Regular Market Demand One System (One Model) Perishable How ?? Figure 3.1: Goal of the Perishable Inventory Management Research for Emergency Response. In this chapter, we first discuss the assumptions and policies we adopt and then study a straightforward extension of the EMQ model on perishable items with zero or a “small” min- imum inventory requirement and quantify the “small” amount. Next, we propose a modified EMQ model for perishable items with a minimum inventory volume constraint to address the special requirements for our emergency response context. Then the necessity of enforcing a max- imum inventory cycle constraint on an important system parameter to guarantee the freshness of the stockpile is demonstrated. We investigate the effect of specified policies on the inventory 31 plot of this inventory system with massive perishable items for different production batch sizes (Q). We still maintain the idea of an EMQ model because of its simplicity. We decompose the cost of maintaining such a system into four components: inventory holding costs, fixed ordering costs, purchasing costs and salvage costs. With the aid of the inventory plot, we can express the four parts of the total cost as a non-convex and non-smooth function of Q. To minimize the total relevant cost with respect to Q, we formulate an unconstrained optimization problem and present an exact solution procedure to obtain the optimal production batch size (Q ∗ ). At last, we conduct two different numerical experiments with an anthrax attack example. The first experi- ment demonstrates the advantage of our proposed model over a standard model, which runs two separate systems to meet the regular market demand and emergency demand respectively. In the second experiment, we perform sensitivity analysis on those government controlled system pa- rameters to provide some insights for both parties (the firm and the government) as they negotiate the parameter settings. 3.2 EMQ Model with Perishability In this section, we first present the assumptions and policies we adopt, then discuss a simple extension of an EMQ model for perishable items. In this work, we assume a single fixed-life perishable item is produced, consumed and stored for an infinite continuous time horizon. We denote T s as the shelf-life of the inventory (in years). The regular market demand is known with a priori constant rate D per year. The production can start at any time at a constant rate P, which is greater than D; and there is a constant cost A associated with each production setup. The holding cost, unit purchase cost and unit salvage 32 cost are all time invariant. At any point of time, at least a pre-determined non-spoiled minimum inventory volume (I min ) must be maintained. The classical EMQ model can be directly applied to the perishable stockpile with zero min- imum inventory requirement as long as the inventory can be consumed before it expires. The shelf-life sets an upper bound on the regular EMQ cycle to guarantee the freshness as demon- strated in Fig. 3.2 part (a). As in Fig. 3.2, T is the EMQ cycle and T 1 is the production time T T’ T’’ T T’ Imin T T’ Tmin =Imin/D T’’ Tmin =Imin/D T’’ (a). Imin = 0 (b). “small” Imin > 0 Figure 3.2: Perishable Inventory System with I min = 0 or small I min . per cycle with inventory increasing at rate P− D; T 2 is the non-production period with inven- tory decreasing at rate D. (T 1 = D P T , T 2 = P−D P T , T = Q D where Q is the production batch size.) During T 1 period, only T 1 D amount can be consumed and the rest will be stocked and depleted gradually over the T 2 period. If a FIFO policy is exercised, the amount consumed in T 1 is pro- duced over the period T 00 = D P T 1 . Hence, the upper bound of the oldest item at the end of T is 33 T 0 = T 2 + P−D P T 1 = P−D P P+D P T , which must be less than the shelf-life T s ; that is, P−D P P+D P Q D ≤ T s . Thus we can obtain the optimal production batch size for the perishable inventory with zero I min as Q ∗ = min(EMQ ∗ ,T s D P 2 P 2 −D 2 ), where EMQ ∗ denotes the optimal solution of a regular EMQ model. When I min is non-zero and in a small amount, as in Fig. 3.2 part (b), with a FIFO policy, an I min amount is first consumed to meet the market demand over the period T min = I min D . Therefore, the upper bound of the oldest item should be modified to T 0 + T min = P−D P P+D P T + I min D , and the optimal production batch size for the perishable items with I min minimum inventory requirement is Q ∗ = min(EMQ ∗ ,(T s − I min D )D P 2 P 2 −D 2 ). At the same time, I min is upper bounded by T min ≤ T s and T min ≤ T 1 , which is I min ≤ min(T s D, D P Q ∗ ). This quantifies the situation where a simple extension of the classical EMQ model is applicable to a perishable inventory system. However, when such condition cannot be met, a more sophisticated model is required, which is the focus for the rest of this chapter. 3.3 Modified EMQ Model In this section, we introduce the modified EMQ model which incorporates both regular market demand and emergency demand (I min ) for a perishable stockpile. Fig. 3.3 gives an illustration of the inventory plot for the modified EMQ model we propose. We define the inventory cycle (T inv ) as the minimum length of time (in years) that an inventory pattern repeats. We only allow disposing once at the end of each T inv , at which we dispose all the inventory above the I min level so that the exact I min is reached at the beginning of the next T inv . We 34 Perishable Idisposal Ip3 T Tinv Imin Huge Volume (I_min) Last Production Cycle Underlying regular EMQ cycle Figure 3.3: Illustration for the Modified EMQ Model. use a FIFO issuing policy such that we always use the oldest items to satisfy the market demand as well as to dispose. As we discussed earlier in this chapter, the two key concerns in this model: the items are per- ishable while the huge volume of these perishable items must be maintained in the system, present challenges to traditional inventory models and need to be addressed explicitly. We propose the following ordering policy as demonstrated in Fig. 3.3 to accommodate these two requirements. For any given Q, we initially run a regular EMQ cycle (we call them the “underlying regular EMQ cycles” with cycle length T , where T = Q D ) and make some adjustment near the end of the inventory cycle. In order to maintain the freshness of the whole inventory stock in the system, within an inventory cycle, an I min amount of inventory must be produced and the same amount must be consumed, either by means of regular market demand or by disposing at the end of T inv . 35 This is called the stability condition of the system (3.1 and 3.2) and will be discussed in more detail in the next section. From a consumption perspective (3.1), the I min level must be greater than the demand rate times T inv to guarantee the necessity of this model. I min > D· T inv . (3.1) Otherwise, a traditional EMQ model is sufficient to solve the problem by depleting/refreshing the I min amount within T inv with the regular market demand rate as presented in section 3.2. To satisfy the stability condition from a production perspective (3.2), we require that the I min level must be less than the production rate times T inv to guarantee production feasibility. I min ≤ P· T inv . (3.2) In order to produce an I min amount within T inv while satisfying constraint 3.1, a last production cycle may be required to produce the amount to be disposed at the end of an inventory cycle. There are four quantities to be decided in this model: [a] how much to produce in each regular underlying EMQ cycle (the production batch size Q); [b] how long is each regular underlying EMQ cycle (T ); [c] when to initiate the last production cycle; and [d] how long is the last production cycle. 36 [b], [c] and [d] can be uniquely determined by a given Q. Therefore, we have a system with a single independent decision variable Q and an unique inventory plot corresponding to the given Q. A total cost relevant to maintaining such a system can be computed with the aid of the inventory plot at a given Q. Fig. 3.3 provides an illustrative example on the inventory plot. It demonstrates one case out of the 5 different possible cases. Next, we discuss the 5 cases (see Fig. 3.4) and the classification criteria. The classification is based on three criteria: 1. if a last production cycle is needed to replace the extra disposing part of inventory; 2. where the inventory cycle ends relative to a regular underlying EMQ cycle (in the produc- tion period – the uphill region in the inventory graph, or in the non-production period – the downhill region); 3. where the last production cycle starts (if it starts at the same downhill region of a regular underlying EMQ cycle as it ends or not). The disposing policy of the model requires that we dispose all the inventory above the I min level at the the end of each inventory cycle to restart a new cycle with exactly I min . The amount to be disposed is either the part of the minimum inventory (I min ) which cannot be consumed by the regular market demand within an inventory cycle (which is I disposal = I min − DT inv ) or the part of the inventory above the minimum inventory volume requirement generated by a regular EMQ cycle at the end of an inventory cycle (defined as I p3 ), whichever is bigger. Cases 1 and 2 in Fig. 3.4 depict the situation that I p3 is greater than I disposal and we dispose I p3 without a last production cycle. Cases 3, 4 and 5 fit the situation where I p3 is less than I disposal and we 37 Ip3 T Tinv T Tinv Ip3 Idisposal Ip3 T Tinv Idisposal Ip3 T Tinv Idisposal Ip3 T Tinv CASE 1 CASE 5 CASE 4 CASE 3 CASE 2 Figure 3.4: All 5 Cases of Possible Scenarios. 38 initiate a last production cycle to produce the amount I disposal − I p3 and dispose exactly I disposal at the end of each T inv . We further classify the cases based on where an inventory cycle ends relative to the underlying EMQ cycle. In cases 2 and 4, the inventory cycle ends at an uphill region (the production period); otherwise, in cases 1, 3 and 5, the cycle ends in a downhill region (the non-production period). Furthermore, the 3 cases which have a last production cycle, are further classified depending on where the last production cycle starts. If it starts at the same non- production period as it ends, it is case 5; otherwise, if it starts at some earlier regular EMQ cycle, we have cases 3 and 4. With the above three criteria, we can distinguish these 5 cases uniquely. We illustrate the classification process in Fig. 3.5. These enumerate all possible scenarios. Case 1, 2, 3, 4 and 5 If last production cycle is needed? If Tinv ends at a uphill region? If Tinv ends at a uphill region? If last production starts at the same downhill region? Case 1 Case 2 Case 4 Case 3 Case 5 NO NO NO NO YES YES YES YES Figure 3.5: Classification Criteria for 5 cases. 39 3.4 Constraint on the Inventory Cycle Length Parameter In this section, we present a crucial constraint on one important parameter (T inv ) in this perishable inventory system, to maintain the whole stockpile within the shelf life throughout the time hori- zon. Recall that the stability condition of the system states that to guarantee the freshness of the stockpile, the system requires refreshing the entire minimum inventory volume within each T inv . From the consumption’s perspective, within an inventory cycle, the minimum inventory volume must be consumed, either by the regular market demand or by disposing the remaining part at the end of the inventory cycle. From the production’s perspective, the minimum inventory volume must be produced within an inventory cycle. The amount of the minimum inventory produced in one cycle is completely consumed in the next cycle. This is equivalent to claim that, in any T inv , the system must produce I min to be used in the next cycle and consume I min which is produced in the previous cycle. At the beginning of each cycle (which is also the end of its previous cycle), the system reaches its minimum inventory volume requirement level by disposing any amount above it. The system runs a regular EMQ model with extra production in the last production cycle to produce the remaining amount of stock that cannot be consumed by the market demand and needs to be disposed at the end. Fig. 3.6 demonstrates the dynamics of the stability condition. I_min Consumption I_min Production (i-1)th Tinv I_min Consumption I_min Production i-th Tinv I_min Consumption I_min Production (i+1)th Tinv Figure 3.6: Illustration for the Stability Condition. 40 Since we are adopting a FIFO issuing policy, the I min amount produced in one inventory cycle will not be used until the next cycle starts. We assume a continuous and infinite time horizon, and we have a continuous age-distribution for the I min amount produced in an inventory cycle whose age is between 0 and T inv at the beginning of a cycle. The following constraint on T inv , stated in Lemma I, is able to guarantee the freshness of the stockpile. Lemma I: (Maximum Inventory Cycle Length) The maximum inventory cycle length is at most a half shelf-life T inv ≤ 1 2 T s ; this guarantees that • the complete stockpile is always within its shelf-life; • I min is always younger than half the shelf-life and the salvaged items are aged between the half shelf-life and full shelf-life; • the age distribution of the stockpile repeats itself every T inv . Idisposal Ip3 Tinv T1 Figure 3.7: Graph Illustration for an Extreme Boundary Case. Proof. We prove this by an extreme boundary case. Let us consider the extreme boundary case that the I min is exactly equal to the production rate times the inventory cycle, which means that 41 we have to keep producing all the time to guarantee the amount I min can be refreshed every T inv (see Fig. 3.7). That is I min = P· T inv . (3.3) Without loss of generality and for the convenience of the proof, we discretize the time horizon. That is, we divide an inventory cycle into T inv discrete periods from 0 to T inv − 1. Then we can express the age distribution of the I min at the end of an inventory cycle (we call it cycle 1) as follows, which represents that P units of inventory have been produced at every discrete period in cycle 1: Age= 0 P units; 1 P units; ... ...; T inv − 1 P units. (3.4) We observe the age distribution when it moves to the next inventory cycle (we call it cycle 2). At time 0 in cycle 2, we have P units of age T inv which were produced at time 0 in cycle 1 and we can use D units of them to fulfil the market demand. For this proof, we will deviate from our FIFO service (issuing) policy and leave the P− D unused units in stock. At time 1 in cycle 2, we again get P units of age T inv which were produced at time 1 in cycle 1; same as time 0, we can use D units for the market demand and leave P− D units unused in stock. We can continue this process along the discrete time horizon until T inv −1 in cycle 2. At time T inv −1 in cycle 2, before 42 we dispose the unconsumed inventory produced in cycle 1, we have a total of T inv ∗(P− D) units of stock to be disposed, with the following age distribution: Age= T inv for P-D units; T inv + 1 for P-D units; ... ...; 2T inv − 1 for P-D units. (3.5) We use the upper bound on cycle length 2T inv ≤ T s to ensure that all stocks left after disposing have the age distributed in [0, T s ]. We claim that this is the upper bound since we break the FIFO issuing policy in the process described above. Using a FIFO policy, we would have a younger inventory to dispose but the first I min will still have the same age distribution. If we abide with the FIFO issuing policy, at time 1 in cycle 2 we would first use the oldest left-over P−D units of age T inv + 1, which is produced at time 0 in cycle 1, to satisfy the market demand before we use the age T inv batch, which is produced at time 1 in cycle 1. Hence, we can claim that the oldest age of the disposal units at the end of cycle 2 will be at most 2T inv −1. Thus, a half shelf life is the upper bound for the inventory cycle to guarantee the freshness of the stock in the minimum inventory volume (I min is always younger than half shelf-life) as well as the potential salvage value (the age of salvaged items is always between half shelf-life and full shelf-life). 43 3.5 Total Cost Evaluation and Boundary Conditions In this section, we will discuss the detailed calculation on the total cost as well as the boundary point values for all 5 cases. For each case, the total cost depends continuously only with respect to Q. 3.5.1 Notation We first introduce the notation used in the calculation. Please note that we repeat some notation defined before to make this list complete. A: fixed cost associated with each production setup h : holding cost per unit per year D: demand rate on the regular market per year P: production rate per year T s : shelf-life of the inventory (in years) T inv : inventory cycle, the minimum time-span that an inventory pattern repeats (in years) I min : minimum inventory volume requirement v: purchasing cost per unit w: disposing (salvage) cost per unit. Across all 5 cases, the followings can be expressed by the above defined given parameters. Q: production batch size for a regular EMQ cycle T : the length of a regular EMQ cycle with production batch size as Q T 1 : production time in a regular EMQ cycle T 2 : non-production time in a regular EMQ cycle 44 I max : the maximum inventory level in a regular EMQ cycle Above are the parameters which represent the same notations as in a regular EMQ model. The variables below are notations specially used in our proposed model. I disposal : the minimum inventory amount to be disposed every T inv N: number of complete regular EMQ cycles in a T inv T p3 : the remainder of T inv divided by T I p3 : the inventory level of a regular underlying EMQ at the end of a T inv T disposal : production time for the extra inventory to be disposed, which is max(I disposal − I p3 ,0) δ: the non-production time in T p3 T p1 : from the start of the last production cycle to the end of the current underlying EMQ cycle I p1 : the inventory level at the beginning of the last production cycle N 1 : number of complete regular EMQ cycles within the last production cycle M: number of regular EMQ orders in a T inv Below are notations for different components in the total cost evaluation. TC: total cost in a T inv TC Inv : total inventory carrying cost in a T inv TC N : total ordering cost (or production setup cost) in a T inv TC Purchase : total purchase cost in a T inv TC Salvage : total salvage cost in a T inv T RC: total relevant (w.r.t. Q) cost in a T inv . 45 Now we are ready to calculate the inventory level for the different cases to prepare the total cost computation. Below are the basic quantities which share the same formula across all 5 cases. T = Q D , T 1 = Q P , T 2 = Q( 1 D − 1 P ); I max = Q(1− D P ), I disposal = I min − T inv · D; N =b T inv T c, T P3 = T inv %T = T inv − N· T ; I p3 = D·(T− T p3 ) for cases 1, 3, 5; where T inv ends at non-production time; (P− D)· T p3 for cases 2, 4; where T inv ends at production time. For cases 3, 4 or 5, we have the formula for the non-production time in the last incomplete EMQ cycle (T p3 ): δ= max(T p3 − T 1 ,0). Since we assume that at the end of each inventory cycle, the inventory level is set back to I min . At T inv , the inventory level I p3 is less than the required disposal amount I disposal for cases 3, 4 and 5, the time that is needed to produce the extra disposal amount (I disposal − I p3 ) is: T disposal = I disposal −I p3 P . Note that (I disposal − I p3 ) is just part of the amount produced in the last production cycle which cannot be covered by the production periods in a regular underlying EMQ cycle within an inventory cycle. Hence T disposal only occupies part of the last production cycle with the remaining time used to produce items to meet the regular market demand. Another way to express the classification of cases 3, 4 and 5 by δ and T disposal is: if δ= 0 ⇒ case 4; if δ> T disposal > 0 ⇒ case 5; if T disposal > δ> 0 ⇒ case 3. Next, based on the classification we described above, we present how to calculate the total cost for the 5 cases. 46 First, we can decompose the total cost of maintaining this perishable inventory system within a single inventory cycle into 4 parts: inventory holding cost (TC inv ), fixed ordering cost (TC N ), unit purchase cost (TC purchase ) and salvage cost on the disposal part (TC Salvage ). That is: TC = TC Inv + TC N + TC Purchase + TC Salvage . We first look at the computation on the purchase cost and salvage cost. The total purchase cost is the amount to produce within an inventory cycle times the unit price. The total salvage cost is the amount to dispose at the end of each inventory cycle times the unit salvage value. Hence we have the following formulas: For cases 1 and 2: TC Purchase =(I min + I p3 − I disposal )· v, TC Salvage = I p3 · w. For cases 3, 4 and 5: TC Purchase = I min · v, TC Salvage = I disposal · w. In general, these two parts are straightforward and easy to compute for all 5 cases. For cases 3, 4 and 5, TC Purchase and TC Salvage are fixed and independent of the production batch size (Q); hence they can be removed from the total relevant cost calculation. Next, we look at the computation on the total ordering cost and inventory holding cost, which are more complicated and will be our focus from now on in this subsection. We first give the general calculation formula here and then expand them with respect to Q for each case later. For cases 1 and 2, since there is no last production cycle, the number of orders is the number of complete regular EMQ cycles in an inventory cycle plus 1. And the total inventory carrying cost can be calculated by the area under the inventory plot. For case 1 where the T inv ends in the downhill region, the area under the inventory plot would be N+ 1 regular EMQ triangles minus the cut-off small triangle in the shadow, in Fig. 3.8. For case 2 where the T inv ends in the uphill 47 region, the area would be N regular EMQ triangles plus the small extra triangle in the shadow, in Fig. 3.9. The formula is as follows: TC N =(N+ 1)· A TC Inv = 1 2 (N+ 1)· T· I max − 1 2 (T− T p3 )· I p3 for case 1; 1 2 N· T· I max + 1 2 T p3 · I p3 for case 2. For cases 3 and 4, we can only use the non-production time (T 2 ) of regular EMQ cycles to produce the I disposal − I p3 amount (production time of the EMQ cycle is already in use to satisfy the regular demand D). The number of complete EMQ cycles that would be covered by the last production cycle is: N 1 =b T disposal −δ T 2 c. If M is the number of complete regular EMQ cycles before the last production period starts, then there are M+ 1 orders per inventory cycle. M = N− N 1 ; TC N =(M+ 1)· A. For cases 3 and 4 (see Fig. 3.10 and Fig. 3.12), it is obvious that the last production cycle must start on a non-production period (T 2 ) and the time of this last production cycle during the current T 2 is: T p1 =(T disposal − δ)%T 2 =(T disposal − δ)− N 1 · T 2 . Hence we have the height of the short parallel lateral (E 1 E 2 ) of the trapezoid (E 1 E 2 E 3 E 4 ), which is: I p1 = T p1 · D. The total area under the inventory plot for both cases 3 and 4 is M regular EMQ triangles minus the cut-off small triangle in the shadow plus the area of the trapezoid. The long parallel lateral is exactly the amount to be disposed (which is the length of E 3 E 4 which equals to I disposal ). Hence the total inventory carrying cost can be expressed as: TC Inv = 1 2 M· T· I max − 1 2 T p1 · I p1 + 1 2 (I p1 + I disposal )·(T inv −(M· T− T p1 )). 48 For case 5 (see Fig. 3.14), the last production cycle starts at the same downhill slope as the T inv ends and the area under the inventory plot would be similar to the calculation of cases 3 and 4 except for the width of the cut-off small triangle in the shadow and the width of the trapezoid: N 1 = 0, I p1 =(T− T p3 + T disposal )· D; TC N =(N+ 2)· A=(M+ 1)· A; TC Inv = 1 2 (N+ 1)· T· I max − 1 2 (T− T p3 + T disposal )· I p1 + 1 2 (I p1 + I disposal )· T disposal . Next, we present the detailed calculation on the total cost as a function with respect to Q and the boundary conditions for all 5 cases. Case 1: Idisposal Ip3 T1 T2 Tp3 T Tinv Figure 3.8: Graph Illustration for Case 1. T RC(Q) = I p3 ·(v+ w)+ 1 2 h[(N+ 1)· T· I max −(T− T p3 )· D(T− T p3 )] = ( Q D − T inv + N Q D )· D·(v+ w)+ 1 2 h[(N+ 1) Q 2 D (1− D P ) −D·( Q D − T inv + N Q D ) 2 ] 49 We let T RC(Q)= a· Q 2 + b· Q+ c, then a = 1 2 h·[(N+ 1)·( 1 D − 1 P )− D· (N+ 1) 2 D 2 ] = − h 2 (N+ 1) P ·(1+ N· P D ) b = (N+ 1)·(v+ w)+ h· T inv ·(N+ 1) = (N+ 1)·[(v+ w)+ h· T inv ] The T RC(Q) is a quadratic function of Q and the stationary point of the quadratic curve is ( where a< 0 ): Q ∗ case1 =− b 2a = P[(v+ w)+ h· T inv ] h(1+ N P D ) The boundary points for case 1 on one side is that T inv ends at the downhill region when I p3 is exactly the amount of I disposal , which is I min − T inv D. That is: I min − T inv D= D[(N+ 1) Q D − T inv ]; hence: Q lb = I min N+1 . The other boundary point lies when T inv ends at the point where T p3 = T 1 . That is: T inv = N Q D + Q P ; hence: Q ub = T inv N D + 1 P . Case 2: T RC(Q) = I p3 ·(v+ w)+ 1 2 h(N· T· I max + T p3 · I p3 ) = (T inv − N Q D )·(P− D)·(v+ w)+ 1 2 h[N Q 2 D (1− D P ) +(P− D)·(T inv − N Q D ) 2 ] 50 T1 T2 Tp3 T Tinv Idisposal Ip3 Figure 3.9: Graph Illustration for Case 2. We let T RC(Q)= a· Q 2 + b· Q+ c, then a = 1 2 h·[N( 1 D − 1 P )+(P− D)· N 2 D 2 ] = h 2 N·( 1 D − 1 P )(1+ N P D )] b = − N D ·(P− D)·(v+ w)− h(P− D)· T inv · N 1 D = −N 1 D· P (P− D)· P[(v+ w)+ h· T inv ] The T RC(Q) is a quadratic function of Q and the stationary point of the quadratic curve is ( where a> 0 ): Q ∗ case2 =− b 2a = P[(v+ w)+ h· T inv ] h(1+ N P D ) 51 The calculation of boundary points for case 2 is very similar to case 1. The lower bound of case 5 is the upper bound of case 4 and it is the case when T inv ends at the point where T p3 = T 1 . That is: T inv = N Q D + Q P ; hence: Q lb = T inv N D + 1 P . The other boundary point is that T inv ends at the uphill region when I p3 is exactly the amount of I disposal , which is I min − T inv D. That is: I min − T inv D=(P− D)(T inv − N Q D ); hence: Q ub = (PT inv −I min )D N(P−D) . Case 3: T RC(Q) = 1 2 h[M· T· I max − T p1 · DT p1 )+(I disposal + T p1 · D)·(T inv − M· T + T p1 )] We first compute T p1 , which can be expressed as a linear function of Q: T p1 = T disposal −(T p3 − T 1 )− N 1 · T 2 = I disposal + T inv ·(D− P) P +(N+ 1− N 1 )·( 1 D − 1 P )· Q = β 1 + α 1 · Q Where we define: α 1 = ( 1 D − 1 P )·(N+ 1− N 1 )=(M+ 1)·( 1 D − 1 P ) β 1 = I disposal − T inv ·(P− D) P = I disposal + T inv D P − T inv = I min P − T inv Please note that I min P is the total amount of time for producing within an inventory cycle. Hence, 52 Idisposal Tp1 T1 T2 Ip3 Ip1 Tp3 T Tinv E1 E2 E3 E4 Figure 3.10: Graph Illustration for Case 3. β 1 is the negative of the idle time (non-production time) in an inventory cycle. Hence we have: T RC(Q) = 1 2 h[M· Q 2 D ·(1− D P )−(α· Q+ β) 2 · D+ (I disposal +(α· Q+ β)· D)·(T inv − M· Q D + α· Q+ β)] We let T RC(Q)= a· Q 2 + b· Q+ c, then a = 1 2 h·[M·( 1 D − 1 P )− D· α 2 + α· D·(α− M D )] = − h 2 · M 2 ·( 1 D − 1 P )< 0 b = h·{−α· β· D+ 1 2 [(β· D+ I disposal )·(α− M D )+ αD·(T inv + β)]} = h 2 ·[(2M+ 1)·(T inv − I disposal P − D P · T inv )+ I disposal D ] 53 The T RC(Q) is a quadratic function of Q and the stationary point of the quadratic curve is ( where a< 0 ): Q ∗ case3 =− b 2a = (2M+ 1)·(T inv − I disposal P − D P · T inv )+ I disposal D 2M 2 ·( 1 D − 1 P ) Next, we show the calculation of the boundary point values for case 3. Without loss of generality, Range 1 Range 2 A B C D (M) (N) SD2 SC1 SD1 EA2 Range 3 EB3 SC3 Figure 3.11: Graph Illustration for the Boundary Value of Case 3. we assume that the N and M values are both fixed and as shown in Fig. 3.11, the T inv period ends at somewhere in the downhill slope between C and D; the last production cycle starts somewhere in the downhill region between A and B. From Fig. 3.11, we can infer that only if I disposal falls in region 1 or region 2 or region 3, it is feasible for case 3 with the parameter M and N valid (that is, T inv ends between C and D and the last production cycle starts between A and B). Otherwise, if I disposal is too short (below region 1) or too high (above region 2) with fixed N, the value of parameter M would increase or decrease accordingly. 54 When I disposal falls into region 1, the two boundary points of case 3 with M and N are: (1) T inv ends at D, and the last production cycle starts at SD1; (2) T inv ends at C, and the last production cycle starts at SC1. Or when I disposal falls into region 2, the two boundary points are: (1) T inv ends at D, and the last production cycle starts at SD2; (2) T inv ends at EA2, and the last production cycle starts at A. Or when I disposal falls into region 3, the two boundary points are: (1) T inv ends at EB3, and the last production cycle starts at B; (2) T inv ends at C, and the last production cycle starts at SC3. Next, we will calculate the Q values corresponding to these boundary points. For fixed N and M values, when T inv ends at point D, we have (N+ 1)· T = T inv . Since T = Q D D , we have (N+ 1)· Q D D = T inv . Hence Q D = T inv ·D N+1 , where Q D is the production batch size when T inv ends at point D. When T inv ends at point C, we have N· T + T 1 = T inv , so N· Q C D + Q P = T inv . Hence, Q C = T inv N D + 1 P , where Q C is the production batch size when T inv ends at point C. If T inv ends somewhere between C and D and the last production cycle starts at A, we have I disposal − I p3 P =(N− M+ 1)· T 2 +[T inv −(N· T + T 1 )], If we use Q A to represent the production batch size when the last production cycle starts from A, then: I disposal −[(N+ 1)· Q A D − T inv ]· D P = (N− M+ 1)· Q A ·( 1 D − 1 P ) +[T inv −(N· Q A D + Q A P )] 55 Hence we have: I disposal + D· T inv P −(N+ 1) Q A P = (N− M+ 1) Q A D −(N− M+ 1) Q A P − N Q A D − Q A P + T inv I disposal +(D− P)· T inv P = −(M− 1) Q A D +(M− 1) Q A P Q A = (P− D)· T inv − I disposal P·(M− 1)·( 1 D − 1 P ) If T inv ends somewhere between C and D and the last production cycle starts at B, we have I disposal − I p3 P =(N− M)· T 2 +[T inv −(N· T + T 1 )], If we use Q B to represent the production batch size when the last production cycle starts from B, with the similar calculation for Q A , we have: Q B = (P− D)· T inv − I disposal P· M·( 1 D − 1 P ) Case 4: T RC(Q) = 1 2 h[M· T· I max − T p1 · DT p1 )+(I disposal + T p1 · D)·(T inv − M· T + T p1 )] Again, T p1 can be expressed as a linear function of Q: T p1 = T disposal − N 1 · T 2 = I disposal − T inv ·(P− D) P +( 1 D − 1 P )·(N− N 1 )· Q = β 2 + α 2 · Q 56 Idisposal Tp1 T1 T2 Ip3 Ip1 Tp3 T Tinv E1 E2 E3 E4 Figure 3.12: Graph Illustration for Case 4. Where we define: α 2 = ( 1 D − 1 P )·(N− N 1 )= M·( 1 D − 1 P ) β 2 = I disposal − T inv ·(P− D) P Hence we have: T RC(Q) = 1 2 h[M· Q 2 D ·(1− D P )−(α· Q+ β) 2 · D+ (I disposal +(α· Q+ β)· D)·(T inv − M· Q D + α· Q+ β)] We let T RC(Q)= a· Q 2 + b· Q+ c, then a = 1 2 h·[M·( 1 D − 1 P )− D· α 2 + α· D·(α− M D )] = h 2 · M·(1− M)·( 1 D − 1 P )< 0 57 b = h·{−α· β· D+ 1 2 [(β· D+ I disposal )·(α− M D )+ αD·(T inv + β)]} = h· M·[T inv − I disposal P − D P · T inv ] The T RC(Q) is a quadratic function of Q and the stationary point of the quadratic curve is ( where a< 0 ): Q ∗ case4 =− b 2a = T inv − I disposal P − D P · T inv (M− 1)·( 1 D − 1 P ) Next, we calculate the boundary point values for case 4. Similar to the discussion for case 3, Idisposal I’disposal Range 2 Range 1 A B C E (M) (N) SC2 SE1 SE2 EA2 Range 3 I”disposal SE3 EB3 Figure 3.13: Graph Illustration for the Boundary Value Case 4. when I disposal falls into region 1, the two boundary points are: (1) T inv ends at C, and the last production cycle starts at SC1; (2) T inv ends at E, and the last production cycle starts at SE1. When I disposal falls into region 2, the two boundary points are: 58 (1) T inv ends at C, and the last production cycle starts at SC2; (2) T inv ends at EA2, and the last production cycle starts at A. When I disposal falls into region 3, the two boundary points are: (1) T inv ends at EB3, and the last production cycle starts at B; (2) T inv ends at E, and the last production cycle starts at SE3. The calculation of the Q values corresponding to these boundary points is as follow: for fixed N and M, when T inv ends at point E, we have N· T = T inv . Since T = Q E D , we have N· Q E D = T inv . Hence Q E = T inv ·D N , where Q E is the production batch size when T inv ends at point E. When T inv ends at point C, we have N· T + T 1 = T inv , so N· Q C D + Q P = T inv . Hence, Q C = T inv N D + 1 P , where Q C is the production batch size when T inv ends at point C. When T inv ends somewhere between C and D and the last production cycle starts from A, we have I disposal − I p3 P =(N− M+ 1)· T 2 If we use Q A to represent the production batch size when the last production cycle starts from A, then: I disposal −[T inv − N· Q A D ]·(P− D) P =(N− M+ 1)· Q A ·( 1 D − 1 P ) 59 Hence we have: I disposal −(P− D)· T inv P + N Q A D (P− D) P = (N− M+ 1)Q A ( 1 D − 1 P ) I disposal −(P− D)· T inv P = −(M− 1)Q A ( 1 D − 1 P ) Q A = (P− D)· T inv − I disposal P·(M− 1)·( 1 D − 1 P ) When T inv ends somewhere between C and D and the last production cycle starts from B, we have I disposal − I p3 P =(N− M)· T 2 If we use Q B to represent the production batch size when the last production cycle starts from B, with the similar calculation for Q A , then: Q A = (P− D)· T inv − I disposal P· M·( 1 D − 1 P ) Case 5: T RC(Q) = 1 2 h{(N+ 1)· T· I max −(T− T p3 + T disposal ) 2 · D +[I disposal +(T− T p3 + T disposal )· D]· T disposal } 60 Both T disposal and T− T p3 + T disposal can be expressed as linear functions of Q: T disposal = I disposal − I p3 P = I disposal + D· T inv −(N+ 1)· Q P = α 2 · Q+ β 2 T− T p3 + T disposal = T− T p3 + I disposal − D·(T− T p3 ) P = (N+ 1)( 1 D − 1 P )· Q− T inv ·(1− D P )+ I disposal P = α 1 · Q+ β 1 Where we define: α 1 = (N+ 1)( 1 D − 1 P ) β 1 = −T inv ·(1− D P )+ I disposal P = −T inv + T inv D+ I disposal P = I min P − T inv α 2 = − (N+ 1) P β 2 = I disposal + D· T inv P = I min P We note that β 2 represents the total production time within an inventory cycle and β 1 denotes the negative of the idling (non-production) time in an inventory cycle. 61 Idisposal T1 T2 Ip3 Ip1 Tp3 T Tinv Figure 3.14: Graph Illustration for Case 5. Hence we have: T RC(Q) = h 2 [(N+ 1)· Q 2 D ·(1− D P )−(α 1 · Q+ β 1 ) 2 · D+ I disposal ·(α 2 · Q+ β 2 )+(α 2 · Q+ β 2 )· D(α 1 · Q+ β 1 )] We let T RC(Q)= a· Q 2 + b· Q+ c, then a = h 2 ·[(N+ 1)·( 1 D − 1 P )− D· α 2 1 + D·(α 1 + α 2 )] = − h 2 N·(N+ 1)·( 1 D − 1 P )< 0 b = h 2 [−2α 1 · β 1 · D+ α 2 · I disposal + D·(α 1 · β 2 + α 2 · β 1 )] = h·(N+ 1)·[(1− D P )· T inv − I disposal P ] 62 The T RC(Q) is a quadratic function of Q and the stationary point of the quadratic curve is ( where a< 0 ): Q ∗ case5 =− b 2a = (1− D P )· T inv − I disposal P N·( 1 D − 1 P ) Here we calculate the boundary points for case 5. Again, only when I disposal falls into either Idisposal Range 2 Range 1 N I’disposal A1 A2 C Figure 3.15: Graph Illustration for Boundary Value of Case 5. region 1 or region 2, it is valid with the fixed parameter N. For both region 1 or 2, the lower bound of Q would be the situation that T inv ends at C, then: T inv =(N+ 1) Q D , hence Q lb = T inv D N+1 . The upper bound of Q depends on which region the I disposal falls into. When I disposal is in region 1, the upper bound is when T inv ends at A1; while when I disposal is in region 2, the upper bound is when T inv ends at A2. At A1, we have: [(N+ 1) Q D − T inv ]· D= I disposal ; so Q ub = I disposal +T inv D N+1 = I min N+1 . At A2, we have: (T inv −N Q D )(P−D)= I disposal ; which is: T inv (P−D)−I disposal = N Q D (P−D). So: Q= T inv D N − I disposal D N(P−D) = D N PT inv −(DT inv +I disposal ) P−D ; hence, Q ub = D N PT inv −I min P−D . 63 3.6 Local and Global Optimality Q Total Cost (1) Non-continuous (2) Non-differentiable Figure 3.16: Graph Illustration for the total cost w.r.t. Q With the total cost calculation formula in section 3.5 different cases, it is straightforward to plot the total relevant cost with respect to Q (see Fig. 3.16). The x-axis is the Q value and the y- axis is the total relevant cost. The below smooth plot represents the total relevant cost with respect to Q for a regular EMQ model without the I min constraint . It is a quadratic function over Q and the optimal Q value can be readily obtained at the stationary point (where the first order derivative is equal to zero). The top irregular plot denotes the total cost with respect to Q for our proposed model. We can see from the plot that the total cost is non-continuous and non-differentiable at those boundary points between two different cases. In this section, we propose an exact solution method to obtain the optimal Q that reaches the minimum total cost with polynomial complexity. We first discuss the local optimality property for each segment in the total cost plot. Then we 64 prove the property which guarantees the global optimality. Finally, we present the complete algorithm and demonstrate its complexity. 3.6.1 Local Optimality D2 A1 B1 A2 C1 D1 B2 O1 C2 O2 Figure 3.17: Graph Illustration for the Local Optimum Point Theorem I: (Local Optimality) Within each case with fixed N and M, the total cost is a quadratic function of Q and the local minimum is either at a boundary point or at the zero first order derivative point. Proof. The quadratic property of the total cost with fixed N and M for each case (corresponding to a segment in the plot) is evident from the calculation in section 3.5. The local minimum lies at either the boundary point with smaller value (left plot in Fig. 3.17) or the zero first order derivative point (stationary point, right plot in Fig. 3.17). 3.6.2 Global Optimality After we obtain the local optimality of each segment, we can compare the local minimum values from different segments and select the lowest one as the global optimum. However, since the 65 number of segments goes to infinity as Q decreases (or N increases), we need to show that only a limited number of segments are required to calculate the local optimal values as potential global optimality candidates. As in the regular EMQ model, the optimal Q value occurs at a point such that the total ordering cost is at the same level as the total inventory carrying cost. These two components are balanced. In our proposed model, we can use the analogy of the optimality condition for a regular EMQ model to explore the global optimum near the region where the ordering cost and the inventory carrying cost are most balanced. It is straightforward to demonstrate that only when the N value is low, cases 1, 2 and 5 may occur since low N leads to large Q, which also implies large regular underlying EMQ cycles that makes I p3 in a scale that is comparable with I disposal . When N gets larger, cases 3 and 4 segments alternate with each other. If we prove that the total cost will monotonously increase as N increases after a threshold value ¯ N, then it implies that only a limited number of local optimum need to be computed as global minimum candidates. Next we will prove this property. Theorem II: (Global Optimality) When N is greater than a threshold value ¯ N, the total cost will monotonously increase as N increases. Proof. We first prove that the total cost of a fixed N value for both cases 3 and 4 is bounded; then we prove that the bound monotonously increases as N increases after a certain threshold value. First, we use the graph to prove the bounds on the total cost with fixed N. Fig. 3.18 illustrates the general situation of the inventory plot with 4 regular EMQ orderings before the last production cycle. If we consider the area under the complete triangle for the last production cycle as fixed, then the rest area for the inventory carrying cost would be the area in shadow in Fig. 3.18. The shadow area is bounded by the shadow area of the three small triangles in Fig. 3.19 from below; 66 and is bounded by the shadow area of the four large triangles in Fig. 3.20 from above. This is evident since the last production cycle triangle is the same for all three situations (Fig. 3.18, 3.19 and 3.20) and the time T inv − T 0 disposal is evenly divided by 3 and 4 respectively in Fig. 3.20 and Fig. 3.19. (Please note that T 0 disposal here represents the length of the last production cycle, which is different from the T disposal defined and used in section 3.3 and 3.5.) Thus the area of one triangle in Fig. 3.18 is smaller than the one triangle area in Fig. 3.20 and is larger than the one triangle area in Fig. 3.19. Since the shadow area in Fig. 3.18 is between the area for 3 and 4 triangles; the upper and lower bounds on it is proved. Idisposal Tinv Figure 3.18: Graph Illustration for the Global Optimum Proof - General Case Next, we prove the monotonous increasing property of the bounds after a certain threshold. For both cases 3 and 4, the purchase cost and salvage cost are independent to Q; only ordering cost and inventory carrying cost are sensitive to Q. We let T triangle = T inv − T 0 disposal denote the time that is occupied by the regular EMQ cycle triangles. The total area for exactly N triangles within T triangle is: Area(N)= 1 2 DT(1− D P )T N where T = T triangle N So: Area(N)= D 2 (1− D P ) T 2 triangle N . 67 Idisposal Tinv T’disposal Figure 3.19: Graph Illustration for the Global Optimal Proof - Lower Bound. Idisposal Tinv T’disposal Figure 3.20: Graph Illustration for the Global Optimal Proof - Upper Bound. As N increases to N+ 1, the total inventory carrying cost decreases the amount of Area(N)− Area(N+ 1) and the total ordering cost increases by A (where A is the fixed cost associated with each production setup). Hence it is proved that the total cost at the boundary points (where T triangle can be exactly divided into integer number of EMQ cycles) is monotonous increasing beyond ¯ N where ¯ N is the smallest integer value satisfied by Area(N)− Area(N+ 1)≤ A. 68 3.7 Exact Solution Algorithm and Complexity Analysis Guaranteed by Theorem I and Theorem II, below we have the complete algorithm to reach the global optimum for the modified EMQ model we proposed in section 3.3. Initialization: 1. Calculate N U 2. Set global minimum g_min = BIG number Main Loop: for each N from 0 to N U 1. Calculate possible cases (at most 5 cases) and their boundary points 2. Record local optimum for each case (l_min N [case i]) 3. Get the local optimum for current N which is the minimum among all possible cases ( l_min N = min (l_min N [case i]) ) 4. If ( l_min N < g_min ) then : g_min = l_min N Figure 3.21: Complete Exact Algorithm for Solving the Modified EMQ Model Theorem III: (Complexity) The exact solution algorithm at most needs to explore 5 ¯ N local minimum points to reach the global minimum, where ¯ N =b q D 2 (1− D P ) T 2 triangle A c. Proof. From the proof in Theorem II, we have Area(N)= D 2 (1− D P ) T 2 triangle N , and the monotonous increasing threshold is given by ¯ N where ¯ N is the smallest integer value satisfied by Area(N)− Area(N+ 1)≤ A. So we have : D 2 (1− D P )T 2 triangle ( 1 N − 1 N+1 )≤ A D 2 (1− D P ) T 2 triangle A ≤ N(N+ 1)≤(N+ 1) 2 N+ 1≥ q D 2 (1− D P ) T 2 triangle A 69 Hence the smallest integer N which satisfies the above inequality is: d q D 2 (1− D P ) T 2 triangle A e− 1 which is: ¯ N =b q D 2 (1− D P ) T 2 triangle A c. For each N, there are at most 5 different cases, thus the maximum number of local minimum points we need to explore before reaching the global optimality is 5 ¯ N. 3.8 Computational Experiments In this section, we conduct numerical experiments for a potential anthrax attack on our proposed model. It serves two purposes. First, we compare our proposed model which combines two types of demand (an emergency demand – the I min minimum inventory requirement, and a regular market demand) into a single system, with a standard model which runs the two parts separately; second, we are interested in investigating how the government controlled parameters affect the firm’s profit and thus it provides the government a perspective to negotiate parameters with the firms producing the drugs. This section is organized as follows. We first discuss the parameter estimation in the ex- periments, then compare two different models (the proposed model and a standard model) and conduct a sensitivity analysis for our proposed model. 3.8.1 Parameter Estimation In a potential anthrax attack, the federal government is prepared to treat 10 million exposed per- sons. This represents a stockpile of 1.2 billion Cipro pills (the treatment regimen is two pills a day for 60 days), as I min used in our experiments. According to the Cipro Pharmarcy.com 70 website (htt p ://www.cipro f loxacinpharmacy.com/active.html): “Cipro has a shelf life of ap- proximately 36 months. However, materiel has been and is currently being tested through the DOD/FDA Shelf Life Extension Program (SLEP) and has received extensions up to 7 1/2 years from original expiration date and some lots have received up to 9 years from original expira- tion date. Materiel shows no signs of deteriorating based on yearly test.” Based on the above statement, we use 9 years (108 months) as the shelf-life of the drug, which indicates that the gov- ernment would pay for the production and storage of I min every 9 years. We define a parameter Y to represent the flexibility that the government gives to firms, which is defined as the number of months before the expiration that the government allows I min to be sold. With the constraint that the maximum cycle length (T inv in our model) must be at most half shelf-life (we use shelf-life minus the flexibility here to represent the actual allowable on-shelf period). In this section, we use the upper limit that T inv = 1 2 (T s −Y)= 108−Y 2 , which dictates the relationship between T inv and Y . We use 36 months for Y in the base case which represents that the government allows firms to resell the required I min 3 years prior to its expiration. Then T inv = 108− 36× 2= 36 months = 3 years. Please note that the government still pays the firm for producing and keeping the inventory every 9 years (108 months) regardless the flexibility (Y ) it allows. We also define p gov as the price that the government pays to the firms per pill for production and storage; p market as the price the firms can sell to the regular US domestic market per pill. We estimate p gov , p market , v and w according to Socolar and Sager (2002). Based on the same source, we assume that the government uses a tier pricing strategy to pay firms, 95 cents listed in table 3.1 represents the price per pill for the first million pills, the second million pills is 10 cent cheaper per pill and from the third million onwards, there is another 10 cent per pill discount. Hence as long as the first million pills’ price is given, the total price that the federal government pays to 71 Table 3.1: Estimated numerical value of the parameters used in the experiment. Parameter Unit Estimation I min million pills 1200 (= 10× 60× 2) T inv month 36 Y month 36 (=108 - 2× T inv ) p gov mil $ /mil pill 0.95 p market mil $ /mil pill 4.67 v mil $ /mil pill 0.2 w mil $ /mil pill -0.3(3 years)/-0.075(6 years) D mil pill / year 300 P mil pill / year 600 A mil $ /time 2 (= 100× h) h mil $ /mil pill/year 0.02 firms is fixed every year. We use 30 cents – the price of generic at Ranbaxy in India (Socolar and Sager, 2002) as the salvage value after 36 months to the secondary overseas market, and 7.5 cents if after 72 months (this will be used in the standard model’s calculation only). The parameters related to the firms’ manufacturing and inventory keeping are approximated according to Singh (2001). “In the US alone, Bayer sold $1.04 billion worth of Cipro in 1999”. This is equivalent to about 220 (≈ 1040 4.5 ) million pills in year 1999; with a 4% per year growth rate, the demand reaches 300 million in year 2007/2008 and we use this as the demand rate D in our model. “At best, Bayer offered to produce 200 million pills within 60 days”. This indicates a maximum production capacity as 1200 million per year, we use half of it – 600 million per year as the production rate allocated for our system. The inventory holding cost (h) mainly comes from the capital cost, with an expected 10% annual investment return rate, h≈ 0.2× 10%= 0.02 mil$/mil pill/year. We assume the ordering cost A is 100 times the h value. Thus, A= 2 for our experiments. In Table 3.1, we summarize the estimated parameter values as we discussed above and we use them as the base case in our experiments. 72 As we combine the two parts of demand into our proposed model: the regular market de- mand with constant rate D and the emergency demand I min required by the federal government, which are two distinct resources that generate revenues. According to our assumption, both rev- enues are independent on the production batch size (Q); hence the total revenue per year is not a function of Q. As the profit equals to revenue minus cost and the revenue is independent on Q, hence minimizing the cost is the same as maximizing the profit. We emphasize this equiva- lence relationship because we used minimizing the cost in the analysis in section 3.5 and we will use the total profit in the following experimental section (3.8.3) to illustrate how the government controlled parameters (I min , p gov and Y ) affect the total profit based on the analysis in 3.5. 3.8.2 Model Comparison Isalvage Imin Ts Cost = $78 Mil/Year Imin Shelflife P = 400 Regular EMQ For Market Demand P = 600 P = 200 Cost = $111.45 Mil/Year Our Model Standard Model System Part (I) System Part (II) Constant Production For Emergency Demand Figure 3.22: Total cost comparison of a standard model and our proposed model We use the base case parameter setting to run our proposed model and the resulting inventory plot is shown in the left of Fig. 3.22. We use the same parameters to run a standard model, which 73 uses two separate systems to meet the two parts of the demand respectively. At the same time, we split the production capacity 600 million pills per year into two parts for these two individual systems. We use 400 million pills per year for the regular market demand with the regular EMQ model and another 200 million pills per year constantly running for the 1.2 billion pills as I min , which is completely refreshed every 6 years and we can obtain 7.5 cents per pill at the secondary market when they are salvaged after 6 years. As the computation is demonstrated below, our model cost 78 million per year to run the combined system which saves 33.45 million a year compared with the standard model. This illustrates the significant benefit (30% cost-saving in a year) by using our proposed model instead of running two separate systems. For our proposed model, there are four parts in the total cost: • Inventory carrying cost: – Emergency part: 1200× 0.02= $24 mil/year; – Regular market part: ( 1 2 × 150× 2+ 1 2 × 300)× 0.02 = 6 million for 3 years, which is $2 million per year; • Ordering cost: 3 times in each T inv (3 years), which is $2 million per year; • Salvage cost: 300×(−0.3)=−90 million dollars in 3 years, which is -30 million dollars per year; • Variable cost: I min ×v= 1200×0.2= 240 million dollars for 3 years, which is $80 million per year. Hence the total annual cost of the proposed system is: 24+ 2+ 2− 30+ 80= $78 million. For the standard model, we calculate the cost for the two systems separately: 74 • Regular market EMQ model: – Inventory carrying cost: q 2hAD(1− D P ))= $2.45 mil/year; – Variable cost: 0.2× 300= $60 mil/year; • Emergency demand: – Inventory holding cost: 1200× 0.02= $24 mil/year; – Variable cost: 1200 6 × 0.2= $40 mil/year; – Salvage cost: −0.075× 1200 6 =−$15 mil/year Hence the total annual cost of the standard model is: 2.45+60+24+40−15= $111.45 million. 3.8.3 Sensitivity Analysis In previous subsection, we demonstrated the advantage of our proposed model over a standard model by using the base case parameter setting. Now we are interested in the sensitivity analysis on the government controlled parameters: I min , Y and p gov . We use different salvage strategies to investigate how the profit changes for different I min , Y and p gov values. All figures in this subsec- tion are plotted as the profit from the manufacturer’s perspective over the parameters controlled by the government. The government can use these plots to negotiate the parameter settings with the firms. Fig. 3.23 and Fig. 3.24 plot the relationship between the total profit and the flexibility (Y ) given by the government at different I min level. The x-axis represents Y which ranges from 0 to 69 months and it corresponds to T inv from 54 to 19.5 months. The y-axis represents the total profit and different plot lines denote different I min levels. 75 Profit (Mil $/year) Y (months) Figure 3.23: Profit v.s. Y with Q max = 1600 at different I min level In Fig. 3.23, we limit the maximum production batch size (Q max ) as 1600 million pills. Based on this figure, we have the following observations. • For any fixed flexibility (Y ), the higher I min level, the more profit the firm can gain; this is due to the extra revenue obtained from the government is higher than the additional cost required for maintaining the extra I min . • For any fixed I min level, the more flexibility allowed (larger Y , which means firms can salvage the I min amount earlier), the higher profit firms can gain; this is due to the extra flexibility given to firms that they can refresh their inventory in a shorter period hence reduce the average running cost of the proposed inventory system. However, the profit increases at different slopes with different Y values. This phenomenon will be explained in a later discussion. 76 • The two dots, which on two different I min level plot lines and at the same profit level, demonstrate the tradeoff between the high flexibility for the firms and the low minimum inventory requirement for the government. If the government would pay less to the firms by requiring a smaller amount of minimum inventory, it must allow more flexibility (longer time before the expiration to salvage the pills) for firms to obtain the same level of prof- itability. Profit (Mil $/year) Y (months) Figure 3.24: Profit v.s. Y which allows salvage at most 300 million pills In Fig. 3.24, we use a different salvage strategy by assuming that we can earn salvage value for the first 300 million disposed pills only. Additional pills will be disposed at a zero salvage value. Compared with Fig. 3.23, this figure shares the same plot shape when Y is low. But the plots start to have a negative slope at the points when the salvage amount reaches 300 million pills. The turning points are at different Y values for different I min plots; the lower the I min , the 77 larger the turning point’s Y value. This is explained by the fact that a less minimum inventory requirement will reach a given salvage amount with higher flexibility (or, a larger Y ); or from the mathematical formula: the salvage amount = I min − DT inv = I min − D 9−Y 2 = I min − 4.5D+ DY 2 . For a fixed salvage amount, a smaller I min comes with a larger Y . When I min is fixed, as the flexibility Y increases, the inventory cycle (T inv ) is reduced; this results in more salvage amount at the end of each T inv . If the salvage amount exceeds the 300 million limit, only the first 300 million is valuable; however, we still need to pay the unit cost for those salvaged stocks above the 300 million. Thus, after the turning point, the more we are forced to salvage, the less profit we gain. Profit (Mil $/year) Y (months) Figure 3.25: Profit v.s. Y with Q max = 1600 at different unit government-paid price Fig. 3.25 illustrates the relationship between the profit and the flexibility (Y ) at different prices the government would pay firms. Again we limit the maximum production batch size (Q max ) as 78 1600 million pills and use the base case I min as 1.2 billion pills. Similar to Fig. 3.23, we have the following observations. • The higher unit price the government pays, the higher profit the firm gains. • The two dots, which on two different p gov level plot lines and at the same profit level, demonstrate the tradeoff between the high flexibility for the firms and the low price for the government. If the government would pay less to the firms by reducing the unit price of a pill, it must allow more flexibility (longer time before the expiration to salvage the pills) for firms to obtain the same level of profitability. Y Profit Regular EMQ Imin Shelflife Ts Imin Idisposal Idisposal Imin Ts Ts Imin (a) (b) (c) (d) Figure 3.26: Illustration of the different slopes in the profit v.s. Y plot Fig. 3.26 provides a graph explanation on the different slopes at different segments on the profit v.s. flexibility (Y ) plot. It uses the base case parameters and gives the inventory plots at 79 the optimal production batch size Q ∗ for different Y values. Through the sensitivity analysis, the marginal profit gained over a fixed small4Y in (a), (b) and (d) segments all come from the cost saving on the regular EMQ part, which is trivial, compared with the marginal profit gain over a fixed small 4Y at the (c) segment, which is from the additional gain due to the extra salvage amount. In summary, in this experimental section, we estimate the model parameters for a potential anthrax attack scenario and use them to quantitatively demonstrate the advantage of our proposed model over a standard model. We also perform a sensitivity analysis on the government controlled parameters to provide the aid for the negotiation on the parameter settings between the federal government and the firms. 80 Chapter 4 Stochastic Routing Problem to Minimize Unmet Demand 4.1 Problem Statement In this chapter, we address a stochastic routing problem motivated by the problem of distribut- ing medical supplies in large-scale emergency response. The distinguishing characteristics of large-scale emergencies are high demand for supplies, low frequency of occurrence, and high un- certainty in many aspects (e.g. when and where they happen). This makes it virtually impossible to predict the occurrence and impact of large-scale emergencies. In this setting, a vehicle rout- ing/dispatching model provides a method to plan for the first responder’s reaction to the emer- gency, and in the event an emergency occurs, a methodology to quickly adjust the dispatching policy to provide an efficient response. In a large scale emergency, we consider that the objective of the routing problem is to min- imize the unsatisfied demand that can occur because resources and first responders are over- whelmed. In fact, unmet demand in an emergency situation can result in loss of life, an impact that outweighs other commonly used VRP objectives such as the travel time or number of vehicles used. The initial supplies in push packages will typically not be sufficient to cover all demands 81 in the recommended response times. Priorities can be given to these demands to maximize the survivability ratio in the affected region. An additional important aspect in meeting the demand in a large-scale emergency is the re- sponse time. Since the administration of the medication to the population within a time-frame makes an appreciable health difference in a large-scale emergency, it is important to associate each demand with a time window, to address situations in which a late delivery directly leads to loss of life, such as in an anthrax attack. The unpredictable nature of large-scale emergencies leads to many stochastic elements in the problem, including stochastic demand and travel times. First, at a given demand point (e.g. a neighborhood block), the quantity of required emergency supplies (antidotes, protective equip- ments, medication, etc.) is often proportional to the unknown size of the population and/or num- bers of casualties. Second, the casualties exposures, or demands among “worried well” are hard to accurately predict. In addition, in an emergency situation, the traffic condition can become highly uncertain, due to unpredictability in people’s behavior. At the same time, the emergency event itself may directly affect road conditions. For example, an earthquake may destroy a seg- ment of the freeway. When we consider the problem of routing inventory in response to large-scale emergencies, both routes and delivery quantities are important in creating an effective dispatching plan. These two aspects of the solution are related to the following uses of efficient inventory routing plans for a large-scale emergency: First, from a planning perspective, efficient routes can be used by agencies in mock trial runs of an emergency event to pre-design the routing policies and provide training opportunities. Second, from an operational perspective, the models can be used to update the pre-designed routes and promptly determine exact delivery quantities to dispensing 82 or treatment sites to generate real-time near-optimal dispatching solutions after an emergency situation has happened. 4.2 Two Stage Solution Process We decompose our problem into a two-stage problem. In the first planning stage, we formulate the problem into a mixed integer programming (MIP) model well in advance of any possible emergency to generate pre-planned routes, which can be used for mock trial runs and training purposes. In the second operational stage, at the time of the emergency when more information is revealed, we need to quickly respond to the event and generate the delivery requirements with the planned routes as well as making adjustments on the routes if necessary. We refer to the second stage model as the recourse strategy. In this work, we present three recourse operational models that differ in how much re-optimization we can do in the operational stage. We consider an LP recourse strategy in which we can only adjust the quantity sent on each vehicle, a knapsack re- course strategy which, in addition to setting the quantity, allows to modify the routes by skipping customers in each route, and for comparison purposes we also consider a re-planning strategy in which we look for the optimal routing solution for the realized demand. The operational stage problem should balance the need for familiarity with a planned strategy, which can be used by responders for training purposes, with efficiency in the solution and a quick solution time for a rapid response. 83 4.3 A Priori Routing Problem In this section, we address how we mathematically model the problem at the planning stage to accommodate the special requirements of routing for large-scale emergencies. As mentioned earlier, we aim to minimize the unmet demand to maximize the life-saving. Such a problem is related to the routing problem with profits, since minimizing the unmet demand maximizes the collected profits when the profit of each node is its demand level. The route generation is also constrained by deadlines associated with every demand node. However, in our problem, the unmet demand is not only caused by the deadlines, but also potentially by the vehicle capacity and the insufficient supply at the depot. In our problem setting, the large scale emergency makes both the demand level and the traffic condition highly unpredictable; hence uncertainty exists both in demand and in travel time. The overwhelming demand with limited resources and the urgency of the timely delivery make the problem constrained by both the vehicle capacity (in demand dimension) and the service deadline (in temporal dimension). We formulate the route planning stage as a stochastic programming problem where we quantify all the unmet demand, due to either insufficient capacity, supply, or late delivery. We use chance constraints to represent the uncertainty in demand and travel times. In addition to minimizing the unmet demand, the planning problem also considers a sec- ondary objective to guide the construction of a complete route. Instead of ignoring the nodes not selected by the routing problem with profits due to tight deadline, insufficient supply, or limited vehicle capacity, we include them at the end of each route according to this secondary objective. These complete routes that visit all demand nodes even after the deadline or with an empty truck, provide a blueprint for recourse actions in the operational stage. The secondary objective is the 84 arrival time at each node. We use a significantly small coefficient κ to weigh this secondary objective so it will not interfere with the route generation for nodes that can be satisfied within constraints, which is guided by the minimum unmet demand primary objective. We illustrate this model for the planning stage with the following figures. Figure (a) in Fig. 4.1 shows the routes generated by a classical VRP, which visit all the demand nodes and aims to minimize the total travel time. Figure (b) in Fig. 4.1 demonstrates the routes obtained from a routing problem with profits on the same topological graph, which visit the nearby high-profit nodes within the given deadline while ignoring those far-away low-profit ones. Figure (c) in Fig. 4.1 illustrates the result by the routing problem to minimize the unmet demand and then complete the route as we just proposed. It generates exactly the same result as the routing problem with profits before the deadline. After the deadline, it visits the remaining far-away nodes guided by the secondary objective. 4.3.1 Notation We consider a set K of vehicles and a set D of demand nodes. We identify an additional node, node 0, as the supply node (depot) and let C= D∪{node 0} represent the set of all nodes. Indexed on sets K and C, we define the following deterministic parameters: n: initial number of vehicles at the supply node (depot) s : amount of supplies at the supply node (depot) c k : load capacity of vehicle k dl i : service deadline at demand node i. 85 25 10 30 15 DEPOT DEPOT 15 30 10 25 15 5 20 15 5 20 (b) : Example of a Routing Problem with Profit (a) : Example of a Vehicle Routing Problem DEPOT 15 30 10 25 15 5 20 (c): Example of a Routing Problem to Minimize Unmet demand Figure 4.1: A priori route model at the planning stage. We use M as a large constant used to express nonlinear relationships through linear con- straints. We also consider the following two parameters to represent the uncertain travel time and demand, respectively τ i, j,k : time required to traverse arc(i, j) for vehicle k ζ i : amount of demand needed at node i. Finally, we define the binary and non-negative decision variables as follows, indexed on sets K, C: Binary: 86 X i, j,k : flow variables, equal to 1 if (i, j) is traversed by vehicle k and 0 otherwise S i,k : service variables, equals to 1 if node i can be serviced by vehicle k Non-negative: Y i, j,k : amount of commodity traversing arc(i, j) using vehicle k U i : amount of unsatisfied demand of commodity at node i T i,k : visit time at node i of vehicle k δ i,k : delay incurred by vehicle k in servicing i. 4.3.2 Deterministic Model The deterministic, minimize unmet demand problem can be expressed as follows DP : minimize ∑ i∈D U i + κ ∑ i∈D,k∈K T i,k subject to constraints (4.1)−(4.17), where the constraints are explained in detail below. The objective of model DP is to minimize the weighted sum of the total unmet demands over all demand nodes and the total visit time at demand nodes of all vehicles. The κ value usually is set to be very small to make the total travel time a secondary objective compared with the unmet demand quantity. However, the travel time is a necessary term in the objective function to guide the route generation after the deadline. Since we model the routing problem in response to a large-scale emergency, the service start times (arrival times) directly associate with when the supply will be shipped and used at the dispensing sites. We would like to serve the dispensing sites as early as possible for life-saving purposes, so the arrival time is a much more important 87 indicator of the service quality than the conventional objectives such as travel times or operational time. We group the constraints into four parts: route feasibility constraints, time constraints, de- mand flow constraints and node service constraints. The following constraints (4.1) - (4.6) char- acterize the vehicle flows on the path and enforce the route feasibility. ∑ i∈D ∑ k∈K X 0,i,k ≤ n (4.1) ∑ i∈D ∑ k∈K X i,0,k ≤ n (4.2) ∑ j∈D X 0, j,k = ∑ j∈D X j,0,k = 1 (∀k∈ K) (4.3) ∑ j∈C ∑ k∈K X i, j,k = 1 (∀i∈ D) (4.4) ∑ j∈C ∑ k∈K X j,i,k = 1 (∀i∈ D) (4.5) ∑ j∈C X i, j,k = ∑ j∈C X j,i,k (∀i∈ D k∈ K) (4.6) Constraints (4.1) and (4.2) specify that the number of vehicles to service must not exceed the available quantity ready at the supply node at the beginning of the planning horizon. The number of vehicles to service is stated by the total number of vehicles flowing from and back to the depot. Constraint (4.3) represents each vehicle flow from and back to the depot only once. Constraints (4.4) and (4.5) state that each demand node must be visited only once. Constraint (4.6) requires that all vehicles who flow into a demand point must flow out of it. We here highlight that the two constraints (4.4) and (4.5) can be relaxed to address split delivery (“=” is replaced by “≥” ). These two constraints state that each demand node must be 88 visited at least once, which allows split delivery, e.g. multiple visits to a demand node. This is an extension of the classical VRPs and has been investigated by different groups of researches (Dror and Trudeau, 1989, 1990; Dror, 1994; Archetti et al., 2006). Dror and Trudeau (1989, 1990) have analyzed the saving generated by allowing split deliveries in a VRP. Archetti et al. (2006) demonstrated a tabu search heuristic to solve problems up to 199 nodes; most instances were solved less than 10 minutes and a few cases used around 3 hours. In a large-scale emergency it is likely that the demand at points which are seriously affected will not be able to be satisfied by a single truckload. Hence we consider split delivery, which also provides greater solution flexibility to potentially service more demand nodes before the deadline. Constraints (4.7)-(4.10) guarantee schedule feasibility with respect to time considerations. T 0,k = 0 (∀k∈ K) (4.7) (T i,k + τ i, j,k − T j,k )≤(1− X i, j,k )M (∀i, j∈ C k∈ K) (4.8) 0≤ T i,k ≤ ∑ j∈C X i, j,k M (∀i∈ D k∈ K) (4.9) 0≤ T i,k − δ i,k ≤ dl i ∑ j∈C X i, j,k (∀i∈ D k∈ K) (4.10) The fact that all vehicles leave the depot at time 0 is specified by constraint (4.7). Constraint (4.8) enforces the time continuity based on the node visiting sequence of a route. Constraint (4.9) sets the visit time to be zero if the vehicle does not pass a node. The variable δ i,k represents the delay of the visit time if a vehicle reaches the node after the deadline and is set to zero if it arrives before the deadline in constraint (4.10). 89 This model primarily accommodates the emergency situation where late deliveries could lead to fatalities. To maximize the likelihood of saving lives, medication should be received by the affected population within the specified hours of the onset of symptoms to impact the patient survival. This is the rationale behind the preference of using a hard deadline constraint instead of the soft deadline. However, for problems where late deliveries are possible we can translate the proposed model to soft deadlines, having the penalty on the violation represent the worsening in patient condition due to late arrival. Constraints (4.11)-(4.13) state node service constraints. δ i,k ≤(1− S i,k )M (∀i∈ D k∈ K) (4.11) S i,k ≤ ∑ j∈C X i, j,k (∀i∈ D, k∈ K) (4.12) S i,k M≥ ∑ j∈C Y j,i,k − ∑ j∈C Y i, j,k ! (∀i∈ D, k∈ K) (4.13) Binary decision variables S i,k are used to indicate whether a node i can be serviced by vehicle k (when it equals to 1). That is, if the vehicle k visits node i before the deadline, then the vehicle can drop off some commodities at this node. However, the vehicle does not necessarily do it when S i,k equals to 1 since there might not be enough supply at the depot so the vehicle may not carry any commodities when it visits a later node in the route. We use these binary variables to keep the feasible region of this problem non-empty all the time. Constraints (4.4) and (4.5) will still enforce each node to be visited once and only once no matter before or after the deadline; however, those visits after the deadline cannot service the node any more. Constraint (4.11) states the deadline constraint and it can only be violated when S i,k equals to zero. Constraint (4.12) 90 illustrates the relationship between the binary flow variables and the binary service variables. It implies the service variable can only be true when a vehicle physically passes a node. Constraint (4.13) requires that no commodity flows in a node after the deadline. On the other hand, there is no compulsory dropping-off commodities at nodes visited before the deadline since there may not be enough supplies to meet the demand. It establishes the connection between the commodity flow and the vehicle flow. Constraints (4.14)-(4.16) state the construction on the demand flows. s− ∑ k∈K " ∑ j∈C Y 0, j,k − ∑ j∈C Y j,0,k # ≥ 0 (4.14) ∑ k∈K " ∑ j∈C Y j,i,k − ∑ j∈C Y i, j,k # +U i − ζ i ≥ 0 (∀i∈ D) (4.15) X i, j,k c k ≥ Y i, j,k (∀i, j∈ C, k∈ K) (4.16) Constraint (4.14) requires the total shipment of commodity from the depot not exceeding its current supply inventory level. Constraint (4.15) enforces the balanced material flow requirement for the demand nodes. Constraint (4.16) allows the flow of commodities as long as there is sufficient vehicle capacity. It also connects the commodity flow and the vehicle flow. X i, j,k ,S i, j binary; Y i, j,k ≥ 0; U i ≥ 0; T i,k ≥ 0; δ i,k ≥ 0. (4.17) Constraint (4.17) states the binary and non-negativity properties of the decision variables. 91 4.3.3 Stochastic Model There are many ways to address parameter uncertainty in optimization problems, leading to dif- ferent models and requiring different information on the uncertainty. In particular, we can obtain an answer to a routing problem under uncertainty through the solution of a representative deter- ministic problem, by using stochastic programming and chance constrained models, by a robust optimization approach, or a markov-decision process model. The type of information available on the uncertain parameters is key in determining which model is most appropriate for a given application. In this work we ignore this application specific concern and assume that any infor- mation that is needed is available. The chance constrained programming relies on probabilistic information, supposed to be available, and tries to find a solution which is “optimal” in a proba- bilistic sense. By adopting this approach, we believe that it is not necessary to look for a solution that is always feasible, since the worst case (if there is one) is very unlikely to occur. However, Robust optimization aims generally to provide the best solution feasible for all the uncertainty considered. The parameters τ i, j,k in constraint (4.8) and ζ i in constraint (4.15) represent the uncertain travel time and demand parameters of our problem, respectively. If we ignore the uncertainty and replace these random quantities by representative values, such as their mean μ τ i, j,k and μ ζ i or mode values, we can solve a deterministic problem DP to obtain a simple solution for this problem. This deterministic solution will be helpful as a benchmark to compare the quality of routes and demonstrate the merits of other more sophisticated methods we discuss next. There are two other ways to handle uncertainty that for this problem lead to the solution of a single deterministic prob- lem DP: chance constrained programming and robust optimization. The solution of this routing 92 problem through other methods of representing uncertainty, such as stochastic programming and markov-decision processes require more involved solution procedures and will not be explored in this paper. In chance constrained programming (CCP) we assume that the parameters τ i, j,k and ζ i are unknown at the time of planning but follow some known probability distributions. We assume they are uniformly and independently distributed. We let α D and α T represent the confidence level of the chance constraints defining the unmet demand at each node and the arrival time of each vehicle at each node respectively. Thus, the constraints with stochastic parameters must hold with these given probabilities. For a given distribution on τ i, j,k and ζ i , we can rewrite constraint (4.8) and constraint (4.15) in the chance constrained fashion with levels α T and α D as follows: P τ|(T i,k + τ i, j,k − T j,k )≤(1− X i, j,k )M ≥ 1− α T (∀i, j∈ C k∈ K) (4.18) P ( ζ| ∑ k∈K " ∑ j∈C Y j,i,k − ∑ j∈C Y i, j,k # +U i − ζ i ≥ 0 ) ≥ 1− α D (∀i∈ D) (4.19) We call this chance constrained model (CCP model), which is modified based on the DP model in section 3.2, by replacing constraints (4.8) and (4.15) with constraints (4.18) and (4.19). Under some assumption of their distribution, constraint (4.18) and constraint (4.19) can be transformed to their deterministic counterpart. From this point onward in this paragraph, we use short no- tation τ and ζ to substitute τ i, j,k and ζ i for simplicity. For example, we assume τ and ζ follow a lognormal distribution with mean μ τ and standard deviation σ τ and mean μ ζ and standard de- viation σ ζ respectively. The logarithm log(τ), log(ζ) are normally distributed as normal(μ 0 τ ,σ 0 τ ) and normal(μ 0 ζ ,σ 0 ζ ). The relationship between the parameters of lognormal distribution and nor- mal distribution is stated as: μ 0 = logμ− 1 2 σ 02 , σ 02 = log( μ 2 +σ 2 μ 2 ). We let κ T and κ D represent 93 the Z value for the normal distribution corresponding to the confidence level α T and α D and we call them “safety factors” in the later experimental results section. Therefore, the deterministic counterpart of constraint (4.18) and constraint (4.19) can be expressed as: (T i,k + e μ 0 τ i, j,k +κ T σ 0 τ i, j,k − T j,k )≤(1− X i, j,k )M (∀i, j∈ C k∈ K) (4.20) ∑ k∈K " ∑ j∈C Y j,i,k − ∑ j∈C Y i, j,k # +U i ≥ e μ 0 ζ i +κ D σ 0 ζ i (∀i∈ D) (4.21) The robust optimization model assumes that the uncertain parameters τ i, j,k and ζ i are only known to belong to a given uncertainty setU. The robust optimization approach, introduced by Ben-Tal and Nemirovski (1998), and more recently extended to integer programming by Bert- simas and Sim (2003) and VRP by Sungur et. al (2006), requires that the solution satisfy the constraints with uncertain parameters for all possible values in the uncertainty setU. That is, we rewrite constraint (4.8) and constraint (4.15) as follows: T i,k + τ i, j,k − T j,k ≤(1− X i, j,k )M ∀τ i, j,k ∈U, (∀i, j∈ C k∈ K) (4.22) ∑ k∈K " ∑ j∈C Y j,i,k − ∑ j∈C Y i, j,k # +U i − ζ i ≥ 0 ∀ζ i ∈U, (∀i∈ D) (4.23) These infinitely many constraints, indexed over every τ i, j,k ∈U and ζ i ∈U can be expressed by single constraints by substituting the maximum possible uncertainty parameters τ ∗ i, j,k = max τ i, j,k ∈U τ i, j,k and ζ ∗ i = max ζ i ∈U ζ i . With these maximum uncertainty parameters, these robust constraints can be written as: T i,k + τ ∗ i, j,k − T j,k ≤(1− X i, j,k )M (∀i, j∈ C k∈ K) (4.24) 94 ∑ k∈K " ∑ j∈C Y j,i,k − ∑ j∈C Y i, j,k # +U i − ζ ∗ i ≥ 0 (∀i∈ D) (4.25) As we mentioned earlier, both the chance constrained model and the robust optimization model can be expressed as a deterministic equivalent problem, with changing demand and travel time values. Therefore these problems are just as difficult to solve as a single deterministic routing problem. The only difference between these two models of addressing uncertainty is the value of the uncertainty used, for instance e μ 0 τ i, j,k +κ T σ 0 τ i, j,k versus τ ∗ i, j,k for the time continuity constraints. Note that varying κ T and κ D we can represent τ ∗ i, j,k and ζ ∗ i values as well as the mean values μ τ i, j,k and μ ζ i , thus in the remainder of the paper we only consider the chance constrained model for different confidence levels. 4.4 Recourse Strategies After an emergency occurs we can use the observed outcome of the uncertain parameters to make a quick adjustment to the preplanned routing solution to respond to the event. This response needs to provide the delivery requirements (how much to load on each vehicle) based on the planning result as well as adjustments to the routes to make the delivery more effective. This second stage model must be able to be solved quickly so that the adjustments do not delay the emergency response. We present three different recourse strategies and their corresponding mathematical models in the following subsections. 4.4.1 Notation In the second operational stage, the uncertain parameters in the planning stage (τ i, j,k and ζ i ) have been revealed and become deterministic as defined below: 95 t i, j,k : actual time required to traverse arc (i, j) for vehicle k d i : actual amount of demand needed at node i. Since the planned route is assumed as given at the operational stage, the vehicle flows X i, j,k become parameters. In our first recourse strategy, the value of S i,k is uniquely determined once X i, j,k are given. Hence, they also are no longer decision variables. However, in our second recourse strategy, we introduce a new binary variable Z i, j,k : flow variables, equal to 1 if (i, j) is traversed by vehicle k and 0 otherwise and establish its relationship with S i,k . Hence, S i,k remain as variables. 4.4.2 LP Recourse In the first recourse strategy, we strictly follow the pre-planned route obtained in the first stage but we can adjust the supply carried by each vehicle and the amount to be dropped at different demand nodes according to the revealed demand level and travel time. In this approach, the commodity flow variables become the only amount to be determined. This is the easiest strategy with least flexibility among the three we propose in this paper. Since it is a linear programming problem, it can be solved very efficiently by a CPLEX solver. Model DLP: minimize ∑ i∈D U i + κ ∑ i∈D,k∈K T i,k 96 subject to: constraints (4.13)-(4.14), (4.16), (4.26)-(4.27). ∑ k∈K " ∑ j∈C Y j,i,k − ∑ j∈C Y i, j,k # +U i − d i ≥ 0 (∀i∈ D) (4.26) Y i, j,k ≥ 0; U i ≥ 0. (4.27) 4.4.3 Knapsack Recourse The alternative recourse strategy is inspired by the classical recourse strategy B in Bertsimas (1992), which assumes the demand level will be revealed before the fleet departs from the depot hence the zero demand customers will be skipped. In our strategy, we share the same assumption that the actual demand and travel time needed will become known before the vehicles leave the depot. Besides deciding the amount carried by each vehicle, the recourse allows for skipping low demand nodes when a direct visit to the next high demand node could result in a less total unmet demand level. This recourse strategy can be modeled into an MIP as below. We use the same set of parameters and decision variables as in the first stage planning model except that X i, j,k becomes a parameter (its value will be given from the pre-planning solution) and we add a new set of vehicle flow variables Z i, j,k . We keep the same objective function and replace X i, j,k with Z i, j,k in the time constraints, node service constraints and demand flow constraints (constraint (4.8) - (4.17)). In the route feasibility constraints group, we will establish the relationship between the new vehicle flow variable Z i, j,k and the pre-planned routes X i, j,k . 97 DP : minimize ∑ i∈D U i + κ ∑ i∈D,k∈K T i,k subject to constraints (4.28)−(4.43), The following constraints (4.28)-(4.32) enforce the route feasibility. ∑ i∈D ∑ k∈K Z 0,i,k = ∑ i∈D ∑ k∈K X 0,i,k (4.28) ∑ i∈D ∑ k∈K Z i,0,k = ∑ i∈D ∑ k∈K X i,0,k (4.29) ∑ j∈D Z 0, j,k = ∑ j∈D Z j,0,k = 1 (∀k∈ K) (4.30) ∑ i, j∈R k ,ibefore j Z i, j,k ≤ ∑ i, j∈R k ,ibefore j X i, j,k (∀k∈ K) (4.31) ∑ j∈C Z i, j,k = ∑ j∈C Z j,i,k (∀i∈ D k∈ K) (4.32) Constraints (4.28)-(4.30) require the number of vehicles dispatched in the second-stage to be the same as the solution obtained from the first-stage pre-planned solution. R k in constraint (4.31) represents the subset of nodes visited by the pre-planned route k. It is defined as R k = {i : i∈ D|X i, j,k = 1}. The expression i before j states the sequence of visiting. Constraint (4.31) allows skipping of some nodes in a given route. Constraint (4.32) guarantees the connectivity of the vehicle flow. Constraints (4.33)-(4.36) enforce schedule feasibility with respect to time considerations. T 0,k = 0 (∀k∈ K) (4.33) (T i,k +t i, j,k − T j,k )≤(1− Z i, j,k )M (∀i, j∈ C k∈ K) (4.34) 98 ∑ j∈C Z i, j,k M≥ T i,k ≥ 0 (∀i∈ D k∈ K) (4.35) dl i ∑ j∈C Z i, j,k ≥ T i,k − δ i,k ≥ 0 (∀i∈ D k∈ K) (4.36) Constraints (4.37)-(4.39) address node service constraints. (1− S i,k )M≥ δ i,k (∀i∈ D k∈ K) (4.37) ∑ j∈C Z i, j,k ≥ S i,k (∀i∈ D, k∈ K) (4.38) ∑ j∈C Y j,i,k − ∑ j∈C Y i, j,k ≤ S i,k M (∀i∈ D, k∈ K) (4.39) Constraints (4.40)-(4.42) state the construction on the demand flows. s− ∑ k∈K " ∑ j∈C Y 0, j,k − ∑ j∈C Y j,0,k # ≥ 0 (4.40) ∑ k∈K " ∑ j∈C Y j,i,k − ∑ j∈C Y i, j,k # +U i − d i ≥ 0 (∀i∈ D) (4.41) Z i, j,k c k ≥ Y i, j,k (∀i, j⊆ C, k∈ K) (4.42) Z i, j,k ,S i, j binary; Y i, j,k ≥ 0; U i ≥ 0; T i,k ≥ 0; δ i,k ≥ 0. (4.43) Constraint (4.43) states the binary and non-negativity properties of the decision variables. The second recourse strategy will lead to a result at least as good as the first one by providing extra freedom on executing the route which allows some skipping on the low demand level nodes while still preserving the advantages of customer-vehicle familiarity and the readiness for training purposes as well. The trade-off is that it requires additional computational effort since the vehicle 99 flow variables are introduced into the model again. However, the flexibility on the vehicle flow is restricted within the fixed route, which dramatically decreases the admissible search space compared with the first stage model. Hence, it can be solved efficiently. We name it as knapsack recourse strategy because we propose an approximation solution approach by solving a knapsack problem. The problem can be solved efficiently by this approximation algorithm and the details of the approach will be discussed in the next section. 4.4.4 Reroute Recourse The last proposed recourse strategy is to solve a deterministic routing problem to minimize the unmet demand used for the first planning stage by plugging in the actual travel time and demand value. The advantage of this strategy lies in that it provides the largest search space and most flex- ibility regarding the current emergency scenario. But it loses the pre-knowledge of the planned routes which makes it difficult for training and driver learning on the routes. Since familiarity of the routes and training are important aspects of reliable solutions, it is impractical to apply this re-plan strategy in real life. However, this strategy can provide us with a bound to compare and evaluate the previous two strategies in the experimental section 4.6.4. 4.5 Solution Approaches 4.5.1 Tabu Search for A Priori Routing Problem Tabu search is a local search procedure which iteratively moves from a solution to its best neigh- bor until some stopping criteria are satisfied. The search keeps a tabu list which prohibits revisit- ing a recently explored node unless some aspiration criteria are met to avoid cyclic movement. It 100 allows a solution to temporarily move to a worse position to escape a local optimum. Tabu search has been successfully used in solving hard optimization problems in many fields. A comprehen- sive review on this technique and its applications can be found in Glover and Laguna (1997). Tabu search was first introduced for solving VRP by Willard (1989). Later on, different groups of re- searchers designed a variety of neighborhood/moves and adopted some problem-specific mech- anism to significantly improve the performance (Osman, 1993; Gendreau et al., 1994; Xu and Kelly, 1996; Toth and Vigo, 1998; etc.). Tabu search has also been applied to major variants of VRP, e.g. VRP with time windows (Taillard et al., 1997), VRP with split delivery (Archetti et al., 2006), the pick up and delivery problem (Bianchessi and Righini, 2006), as well as the stochas- tic VRP (Gendreau et al., 1996). It has been shown that tabu search generally yields very good results on a set of benchmark problems and some larger instances (Gendreau et al., 2002). Applying tabu search to a particular problem requires a fair amount of problem-specific knowledge. Given the success of tabu heuristic in the classical VRP and its variants, and the similar structure of our model, we believe this approach holds much promise in solving our prob- lem. The algorithm we propose here uses some ideas from the standard VRP, and incorporates new features taking into account the unmet demand objective of our problem. 4.5.1.1 Object-Oriented Tabu Search Framework We design an object-oriented Tabu search framework and apply several design patterns (Gamma et al., 1995; Alexandrescu, 2001). Object-oriented (OO) programming is a fast growing conven- tion in the software development and creates a more flexible, elegant, and ultimately reusable design. It uses abstraction to create models based on the real world and utilizes three established 101 paradigms: modularity, polymorphism and encapsulation. In the OR (Operations Research) re- search community, the easily-understood and reusable OO models are inheritably harder than in industry, not only because of the high requirement on the computing performance, but also be- cause of the diversity of the interested topics. However, there are some research communities that aim to adopt the OO design into the implementation, and promote the reuse of basic functional models and libraries (www.coin− or.org). We have noticed that although many research groups have been interested in and implemented the Tabu search for the VRP, the design of an efficient and flexible framework has been rarely dis- cussed in the literature. We would like to fill this gap and provide a detailed OO design for the tabu search framework in solving the routing problem. Fig. 4.2 gives the highest level of the CommandFactory Command InitVehicleCommand InitDepotCommand InitDemandPtCommand InitDeadlineCommand InitEngineCommand RoutePropertyCommand RouteManagerCommand UnmetNodeCommand DepotRemainSupplies Parser VRPTabuParser Connector DataCenter <Singleton> BizObjData Engine (Algorithm) VehicleContainer DemandPtContainer DepotContainer Figure 4.2: The architecture of the Tabu search implementation. architecture of the Tabu search. It separates the system into two parts: the Connector, which pro- vides the interface for outside access; and the DataCenter, which encapsulates the internal control 102 and computation process. The Connector is composed of a CommandFactory and a VRPTabu- Parser. It takes the input information wrapped in a Command, parses the information and then triggers the proper actions of the DataCenter. There are mainly two types of Commands: the first type is for the initial object data information, e.g., the property of the depot, the demand point and the vehicle; the other is for a given initial solution, e.g., the existing routes, the unassigned nodes, etc. The initial solution information is not compulsory. However, it allows for the saving of the initial solution. We define the DataCenter as an Singleton class, which “ensures a class only has one instance, and provides a global point of access to it” (Gamma et. al., 1995). It contains a BizObjData and a Engine. The BizObjData records all the information of the Business Object Data, which only takes the value through the trigger of Commands and cannot be altered during the search process. The Engine is the core optimization component. It completely isolates the data and algorithms needed for the Tabu search procedure within itself. The detailed decomposition of the Engine part is depicted in Fig. 4.3. There are four main components in the Engine: BizObjData, HelperData, InitSolGenerator, and TabuSearch (which is the inheritance of the LocalSearch). The ‘*’ next to the BizObjData represents that this is a ‘smart pointer’ (Alexandrescu, 2001; www.boost.org) to the BizObjData object in the upper level (DataCenter). The same representation also applies to the BizObjData and the HelperData under the SolDataSet. The HelperData contains a DistanceManager, which will be calculated once at the time the Engine is initiated and will not be updated during the later optimization process. The InitSolGenerator produces an initial feasible solution if it is not given. The TabuSearch proceeds with the search process. The details will be discussed in next section. 103 Engine BizObjData * DistanceManager HelperData InitSolGenerator LocalSearch TabuSearch SolDataSet InitSolAlgo BizObjData * HelperData * VRPSolution RouteManager UnassignNodeManager DepotRemainSupplyMger BestSol (VRPSolution) CurSol (SolDataSet) MoveManager TabuTenure Terminator Aspirator Figure 4.3: The architecture of the optimization engine for tabu search. By closing the discussion of the OO design for the tabu search of the routing problem, we would like to briefly summarize the pros and cons of it. First, it separates the data and the operations, and the algorithms and the controls. Such a separation makes the program easy to understand, modify and extend. Second, the highly modular design also facilitates the unit test to identify potential problems and limitations. At the same time, the overhead on the data struc- ture management is inevitable and hence the computational performance (in terms of computing times) may deteriorate. However, with the fast advancement of the hardware platform and the fairly large memory requirement and computational efforts of Tabu heuristic, the overhead of the extra memory consumption and the data structure management are not the bottleneck compared with the tabu algorithm itself. Hence, we believe the OO design will benefit the implementation. 104 4.5.1.2 Tabu Search Procedure Because the problem has insufficient supplies at the depot in most scenarios, we use an Unas- signNodeManager list of nodes to keep track of all the nodes with unmet demand. We also use a RouteManager list of solutions to keep track of all the incomplete (may not include all the nodes) but feasible (meet both capacity and time constraints) routes. We initiate a solution by using a visit-nearest-node heuristic and put all the unvisited nodes into the UnassignNodeManager list. The detailed tabu search procedure is shown in Fig. 4.4. 4.5.1.3 Neighborhood Structure and Tabu Moves A key element in designing a tabu heuristic is the definition of the neighborhood of a solution or equivalently the possible moves from a given solution. Beside adopting the standard 2-opt ex- change move and the λ-interchange move, we implement a new DEM-move to accommodate the special needs of the potentially unmet demands. The DEM-move inserts an unassigned demand node into a current route or exchanges an unassigned node with some (one to three) consecutive node(s) in a route by abiding both deadline and capacity constraints. Fig. 4.5 illustrates some examples of the DEM move. The upper left graph represents an unassigned node and an existing route before the DEM move. The lower left graph gives the result after an insertion DEM move. The unassigned node has been inserted between the first and second nodes of the existing route. The lower left graph is one of the upper left graph’s neighbors. This unassigned node could also be inserted between the second and third nodes or between the third and fourth nodes as long as the deadline and capacity constraints can be met after the insertion. The two graphs at the right hand side show an example of the exchange DEM move. In a λ-interchange move, with λ = 3, 105 Initialize BizObjData (Command driven) Initialize Engine (Command driven) helperData.init() Æ distanceMger.init() initSolGenerator.init() Æ get an initial feasible solution Optimize (Command Driven) tabuSearch.cloneSol (initSolGenerator.getSol()) moveMger.genInitAllMoves() tabuTenure.init() (tabu list is empty) curMove = moveMger.selectBestMove Y N N Y N Y Y N isBestMoveEmpty tabuTenure.isTabu(curMove) tabuSearch.isAspired(curMove) curSol.move(curMove) tabuTenure.updateTenure(curMove) curSol.isBetter(bestSol) bestSol = curSol tabuSearch.isTerminiate() moveManager.updateMoves() END of Tabu, post processing, output result N Y Figure 4.4: Tabu Search Procedure. 106 it includes a combination of vertex reassignments to different routes and vertex interchanges be- tween two routes; the reassigning/interchanging segment is up to three consecutive vertices in a route. (See Fig. 4.6 for the illustration with λ= 1 insertion at the left and λ= 2 exchange at the right.) We also maintain a MoveManager list to record all possible moves from the current solu- tion. After a solution moves to its neighbor, instead of reconstructing the whole neighborhood of the new solution, we only update the non-overlapping neighbors. That is, we eliminate those moves that are relevant (sharing the same route or sharing the same unassigned node, etc.) to the move that just has been executed, and generate new feasible moves that are relevant to the updated route(s) only. This will significantly reduce the computational effort of exploring the neighborhood. Exchange Insertion Insertion DEM Move Exchange DEM Move Figure 4.5: Illustration of the DEM Moves. 107 4.5.1.4 Tabu Tenure, Aspiration and Stopping criteria The Tabu tenure is the number of iterations a recent move keeps prohibited. It means a move cannot be re-executed (or be moved back) within the tabu tenure unless it is “aspired”. In this implementation, we apply a random tabu duration (uniformly random generated between 5 to 10 steps) and a “best-ever” aspiration scheme (if it leads to a new best solution so far). We stop the heuristic when any of the following conditions has been met: (1) after a given number of non-improving steps, and (2) a prescribed number of iterations has been executed. Exchange Insertion Insertion Ȝ -Interchange Move Exchange Ȝ -Interchange Move Figure 4.6: Illustration of the λ-interchange Moves. 4.5.1.5 Post-processing for a Complete Route Construction After we obtain the tabu search solution, we complete the route by adopting a variation of a “next- earliest-node” heuristic to insert all unassigned nodes after the deadline of each route (keeping the before-deadline part intact). In this post-processing heuristic, each unassigned node picks a 108 route to locate itself after the deadline, where the summation of its arrival time and the arrival time to the next node is minimized at the time of assignment. The quality of this tabu heuristic is evaluated in section 4.6.2. 4.5.2 Approximation Heuristic for Knapsack Recourse In the operational stage, the first LP recourse strategy can be solved by CPLEX, while the third re-plan strategy shares the same model as the planning stage, hence can be solved by the tabu heuristic we presented in the previous subsection. Therefore, in this subsection, we focus on presenting an approximation algorithm for the second knapsack recourse strategy. After the actual demand and travel time have been revealed at the beginning of the second stage, the knapsack recourse strategy follows the first-stage pre-planned routes and allows skip- ping low demand level nodes. Hence, there are two decisions to make: (1) which nodes to be skipped and, (2) how much supply to load on each vehicle under the given total supply and dead- line, to minimize the total unmet demand. We propose the following solution procedure: For each pre-planned route, we enumerate all the feasible routes by skipping some nodes and sticking with the original sequence. A feasible route must be finished within the pre-specified deadline. This is equivalent to assigning a binary digit to each node in the pre-planned route, 1 if selected, 0 otherwise. Therefore, there are at most 2 n K feasible routes for each vehicle (where n is the number of demand nodes and K is the number of vehicles). We use L k to denote the set of feasible routes for vehicle k. For every feasible route l for vehicle k in L k , an associated demand level a k,l is the summation of the demand of nodes in the route. With a known total supply quantity at the depot, we can solve an integer programming problem to decide the route to be used for each vehicle. 109 Parameters: s : amount of supplies at the supply node (depot) a k,l : the total demand on route l of vehicle k. Binary Decision variables: W k,l : equals to 1 if route l is used by vehicle k, 0 otherwise. IP : maximize ∑ k∈K,l∈L k a k,l W k,l subject to constraints (4.44)−(4.45), ∑ l∈L k W k,l = 1 (∀l∈ K) (4.44) ∑ k∈K,l∈L k a k,l W k,l ≤ s (4.45) The IP aims to maximize the demand that can be serviced, which is the same as minimizing the total unmet demand. Constraint (4.44) states that only one route can be selected for each vehicle. Constraint (4.45) enforces the total demand to be serviced will not exceed the supply available at the depot. This problem is called a multiple-choice knapsack problem, which is a 0-1 knapsack problem with the addition of the disjoined multiple-choice constraints. Pisinger (1995) presented a simple partitioning algorithm and incorporated it into a dynamic programming framework to solve the problem efficiently. Problems with 1000 different classes of items and each class with 100 items can be solved within 200 seconds. The solution of this IP will specify which route calculated in the first enumeration process should be used (hence we know which node(s) would be skipped); and the associated parameter 110 a k,l states the quantity of supply loaded on vehicle k. The constraints from the temporal perspec- tive are considered in the first enumeration procedure (the definition of a feasible route); and the ones from the demand/supply perspective are stated in the second IP model. 4.6 Computational Experiments In this section, we demonstrate how the proposed models and two stage solution approaches presented in this work are used. We define a problem scenario as a randomly generated network with a single depot and a set of dispersed demand nodes whose mean demand quantity is also randomly generated. The size of a problem scenario is defined as the number of the demand nodes. We conduct three sets of experiments. The first is to evaluate the performance of the tabu heuristic for the first stage model on different size problems by comparing the heuristic results with CPLEX bounds. The second and third are both to perform simulations on 10 different scenarios with size 50 and av- erage the results. Our interests in the second experiment are primarily focused on comparing the alternatives between short-and-risky routes and long-and-safe routes as well as the investigation on how the value of the safety factors influences the quality of the routes. The purpose of the third experiment is to compare the effect of different recourse strategies through simulations. We first describe how to generate the input parameters and the simulation process, then demonstrate and discuss the experimental results. 111 4.6.1 Numerical Setup and Simulation Processes We test our model and heuristic on 10 different problem scenarios. For each scenario we uni- formly generate: 50 demand nodes and 1 depot node in a 200 by 200 square domain; and a mean demand quantity for each demand node ranging from 5 to 15. We service this demand with a fleet of 10 uncapacitated vehicles. To further evaluate the performance of our tabu heuristic, we test it on two smaller problem sizes: with 10 and 20 demand nodes with fleet size as 3 and 4 respec- tively. We use the Euclidean distances between any two of the nodes and assume a symmetric complete graph topology. The mean travel time between any two nodes is proportional to the distance. In the DP model, we use the mean value of the demand quantities and travel times as its parameters. In the CCP model, we use a lognormal distribution with the same mean value as was used in the DP model; the standard deviation is set to be proportional to the mean value for demand (20% of the mean value) and inversely proportional to the mean value of the travel time (σ = UB−μ 100 μ, UB is the upper bound for the inversely proportional transformation, whose value is dependent on the graph topology, which must be greater than the longest arc in the graph). We restrict our analysis to this type of travel time distribution because we aim to compare the tradeoffs between short-and-risky routes and long-and-safe routes. Note that the DP routes in our experiments, tend to be short-and-risky since they do not take into account the variability, while the CCP routes tend to favor longer routes that have smaller variance. For the chance constrained model, we set the confidence level as 95% for a priori route comparison and 85% for the recourse strategies comparison, thus setting the values κ D and κ T in constraints (4.20)-(4.21) to 1.65 and 1.04 respectively, which we refer to as the “safety factor”. Since the routes generated 112 from the model are sensitive to both the deadline and the total supply at the depot, we vary these two parameters and observe the results. We use 70%, 80%, 90%, 100% and 120% of the base quantity as the available supply quantity. The base quantity of the total supply at the depot is defined to be the summation of the demand quantity at all the demand nodes, which is 500 on average. The deadline is set to 40%, 50%, 60%, 80%, 100% and 120% of the base route length. The base route length is defined as the average length of all the edges in the graph times 5 (50 demand nodes are served by 10 vehicles; so on average, each vehicle serves 5 demand nodes). We call one combination of the deadline and the total supply parameters a test case. So for each problem scenario, we have 30 test cases, which are identified by the different combinations of the deadline on the demand nodes (6 types: 40%, 50%, 60%, 80%, 100%, 120%) and the total supply at the depot (5 types: 70%, 80%, 90%, 100%, 120%). For each test case, both the deterministic and the CCP version of the problem are solved by the tabu heuristic. In this simulation experiment, we evaluate and compare the effectiveness of deterministic and chance-constrained a priori routing model as well as three recourse strategies we proposed in section 4.3 and 4.4. We present the process of the simulation as a flowchart in Fig. 4.7. First, with different given deadline and total supply combination and other input parameters (include the stochastic travel time and demand parameters), we use the tabu search heuristic to solve the first stage deterministic or chance-constrained model to obtain the set of routes and commodity flow values. However, at the pre-emergency stage, only the preplanned routes are of interest since they provide the guide for training and preparation purposes; the commodity flow solution which determines how much to load on each vehicle is not very meaningful at the planning stage since the actual demand may vary significantly. Hence we discard the commodity flow solution and only keep the preplanned routes. Then we randomly generate the realization of those stochastic 113 parameters (travel time and demand) following the distribution we used to model the first stage CCP formulation. For each test case, 100 parameter realizations have been generated to simulate different scenarios. The simulated data and the preplanned routes are sent to the three different recourse strategy models to obtain the new commodity flow solution as well as the actual total unmet demand. In summary, we average each test case (a supply and deadline combination) over 10 different problem scenarios with each scenario over 100 replications. Hence the final result is the average value of 1000 raw data points. CCP Para MIP Solver (Tabu Search) X* ccp, Y* ccp Para realization, X* ccp LP Recourse (CPLEX) Discard Y* Y* LP, Total Unmet Y* ks, Total Unmet Average & Compare (Repeat 10) (Repeat 100) (Over 1000) (DL & Sup) Preplan Routes Realization Response Knapsack Recourse Replan (Tabu) Y* replan, Total Unmet Figure 4.7: Flowchart of Simulation Process. 4.6.2 Evaluation of Tabu Search Quality To evaluate the quality of our tabu heuristic, we first run this search process over 10 random problem scenarios with problem size 50. The average of the percentage of unmet demand as well as the remaining supply quantity level at the depot for the deterministic formulation is recorded 114 for each test case. Note that each entry in the table is an average of 10 routes. The results are presented in Table 4.1 and Table 4.2. Table 4.1: Percentage of unmet demand over total demand for the deterministic model. DL200 DL250 DL300 DL400 DL500 DL600 (40%) (50%) (60%) (80%) (100%) (120%) sup350(70%) 30.78% 30.63% 30.63% 30.61% 30.59% 30.59% sup400(80%) 21.19% 20.75% 20.73% 20.71% 20.64% ∗ 20.64% ∗ sup450(90%) 12.76% 11.12% 10.82% 10.82% 10.82% 10.78% sup500(100%) 5.31% 2.20% 1.87% 1.77% 1.75% 1.75% sup600(120%) 4.98% 0.20% 0.00% ∗ 0.00% ∗ 0.00% ∗ 0.00% ∗ Table 4.2: Remaining supply at depot for the deterministic model. DL200 DL250 DL300 DL400 DL500 DL600 (40%) (50%) (60%) (80%) (100%) (120%) sup350(70%) 0.9 0.2 0.2 0.2 0.1 0.1 sup400(80%) 2.4 0.2 0.2 0.1 0 ∗ 0 ∗ sup450(90%) 9.7 1.5 0.1 0.1 0.1 0.1 sup500(100%) 22.1 6.5 4.7 4.4 4.3 4.2 sup600(120%) 120.3 95.9 94.8 ∗ 94.8 ∗ 94.8 ∗ 94.8 ∗ These results show that, in terms of minimizing the unmet demand, the 6 cases which are indicated by a “*”, reach the minimum. They either meet all the demands (0% unmet demand) or deliver all the supplies (zero remaining supply at the depot) over the 10 problem scenarios we tested. From the first three rows of the above tables, the results are very close to the best we can do regarding to the unmet demand, since the remaining supplies are almost empty and the percentage of the unmet demand is near the simple lower bound given by the shortage of supplies to meet the total demand. Only the DL200-sup500 and DL200-sup600 grids give a higher percentage of unmet demand than what is due to the shortage in supplies. Another experiment to evaluate the quality of tabu is to run this search over 5 random problem scenarios over each problem size on 3 sizes (10, 20 and 50). We also run these problems using 115 CPLEX 9.0 for one hour and record the lower bound and upper bound obtained from the CPLEX solver to compare with the tabu result. The results for size 10, 20 and 50 are presented in Tables 4.3, 4.4 and 4.5 respectively. From these tables, we can see that for smaller size problems (size 10), the CPLEX solver can obtain a feasible solution within 1 hour for most of the instances (40 out of the 45 instances). The Tabu heuristic results obtained within 200 seconds generate better solutions than the CPLEX results for about half of the instances (17 out of the 40 instances). The CPLEX solver cannot obtain a feasible solution within 1 hour for all larger size problems (size 20, 50). However, Tabu search can generate a solution within 200 seconds close to the lower bound obtained from the CPLEX. The lower bound provided by the CPLEX solver primarily reflects the shortage of the supply at the depot compared with the total demand. These tables also show that there are bigger gaps between the CPLEX lower bound and the tabu result for the tight deadline and low total supply level combination. From both sets of evaluation experiments, we suspect the big gaps between the tabu results and trivial lower bound from CPLEX at test cases with sufficient supply and short deadline is due to the tight deadline which might prevent the timely delivery even with enough supplies. A tighter lower bound with respect to a short deadline is a subject for future research. However, in general, we conclude that the tabu search method we proposed is very effective in minimizing the unmet demand, which is our primary concern in this model. 116 Table 4.3: Tabu Heuristic Result with CPLEX Bounding for 10 Customers and 3 Vehicles. 10 Nodes 3 Vehs CPLEX = 3600 sec Deadline 150 - (40%) 250 - (80%) 400 - (120%) CPLEX LB CPLEX UB Tabu CPLEX LBCPLEX UB Tabu CPLEX LBCPLEX UB Tabu Supply 70 - (70%) case 1 14633 N/A 14792 14633 N/A 14792 14633 14664 14792 case 2 30523 30537 33002 30522 30538 30913 30522 30537 30913 case 3 40162 40176 41774 40162 40176 41774 40162 40180 41774 case 4 22043 22059 24996 22044 22059 22275 22043 22059 22275 case 5 34762 41121 37764 34762 34784 37764 34763 34781 35964 90 - (90%) case 1 4.11 N/A 14792 3.08 N/A 11.86 3.39 32.01 11.86 case 2 10522 18801 18789 10522 10538 11008 10523 10539 11008 case 3 20162 22984 23133 20162 20179 21912 20162 20181 21912 case 4 2044 29233 17952 2043 2060 6582 2044 2063 6582 case 5 14763 30712 21961 14763 14782 16301 14763 14781 15811 120 - (120%) case 1 3.11 34074 14792 3.30 5964 11.86 3.25 27.60 11.86 case 2 2.30 10536 18789 2.14 31.45 8.16 2.67 16.20 8.16 case 3 1.70 8980 16539 1.85 21.94 8.24 2.10 20.24 8.24 case 4 3.29 N/A 17952 3.15 21.35 11.37 3.65 20.67 11.37 case 5 2.47 37776 21961 2.57 27.70 10.70 3.02 26.46 10.70 Table 4.4: Tabu Heuristic Result with CPLEX Bounding for 20 Customers and 4 Vehicles. 20 Nodes 4 Vehs CPLEX = 3600 sec Deadline 200 - (40%) 400 - (80%) 600 - (120%) CPLEX LB CPLEX UB Tabu CPLEX LBCPLEX UB Tabu CPLEX LBCPLEX UB Tabu Supply 140 - (70%) case 1 36461 N/A 36499 36461 N/A 36498 36461 N/A 36498 case 2 70551 N/A 70851 70552 N/A 70753 70552 N/A 70753 case 3 75171 N/A 75400 75171 N/A 75344 75171 N/A 75344 case 4 56362 N/A 60026 56362 N/A 58375 56363 N/A 58375 case 5 57491 N/A 58152 57492 N/A 59671 57492 N/A 59671 180 - (90%) case 1 1.15 N/A 14822 1.08 N/A 20.78 1.21 N/A 20.78 case 2 30558 N/A 30978 30552 N/A 30846 30552 N/A 30846 case 3 35171 N/A 36337 35171 N/A 35360 35171 N/A 35360 case 4 16362 N/A 53166 16362 N/A 20168 16362 N/A 20168 case 5 17492 N/A 28996 17491 N/A 19945 17491 N/A 19945 240 - (120%) case 1 1.15 N/A 14822 1.15 N/A 20.78 1.21 N/A 20.78 case 2 1.77 N/A 7800 1.83 N/A 20.39 1.69 N/A 20.39 case 3 0.98 N/A 15332 0.96 N/A 21.89 0.99 N/A 21.89 case 4 2.53 N/A 53166 2.53 N/A 22.97 2.51 N/A 22.97 case 5 1.58 N/A 28996 1.42 N/A 23.19 1.59 N/A 23.19 117 Table 4.5: Tabu Heuristic Result with CPLEX Bounding for 50 Customers and 10 Vehicles. 50 Nodes 10 Vehs CPLEX = 3600 sec Deadline 200 - (40%) 400 - (80%) 600 - (120%) CPLEX LB CPLEX UB Tabu CPLEX LBCPLEX UB Tabu CPLEX LBCPLEX UB Tabu Supply 350 - (70%) case 1 139752 N/A 140100 139752 N/A 139810 139752 N/A 139810 case 2 154132 N/A 154190 154132 N/A 154190 154132 N/A 154190 case 3 137861 N/A 137960 137861 N/A 137960 137861 N/A 137960 case 4 147863 N/A 147930 147863 N/A 147940 147863 N/A 147940 case 5 150453 N/A 150710 150453 N/A 150550 150453 N/A 150550 450 - (90%) case 1 39751 N/A 40569 39751 N/A 40125 39751 N/A 40125 case 2 54132 N/A 59333 54132 N/A 54195 54132 N/A 54195 case 3 37861 N/A 38950 37861 N/A 37949 37861 N/A 37949 case 4 47863 N/A 48825 47863 N/A 47924 47863 N/A 47924 case 5 50453 N/A 50766 50453 N/A 50536 50453 N/A 50536 600 - (120%) case 1 1.49 N/A 40569 1.49 N/A 47.23 1.49 N/A 47.23 case 2 1.63 N/A 40470 1.63 N/A 47.77 1.63 N/A 47.77 case 3 1.22 N/A 38950 1.22 N/A 47.55 1.22 N/A 47.55 case 4 3.36 N/A 30781 3.36 N/A 48.34 3.36 N/A 48.34 case 5 2.70 N/A 29562 2.70 N/A 58.99 2.70 N/A 58.99 4.6.3 Evaluation of A Priori Route 4.6.3.1 The Deterministic Routes vs the CCP Routes In this section, we first evaluate the efficiency and robustness of our chance constrained model in determining pre-planned routes. We compare solutions from the CCP model against routes generated by the DP model that uses the mean values of the demand quantity and travel times as problem parameters. As the simulation process we described in section 4.6.1, we compute the DP and CCP routes (which do not allow split delivery) with the tabu heuristic for 30 test cases over each one of 10 problem scenarios and then generate 100 instances of the realization. We evaluate how well the DP and CCP routes can meet the demand in each instance by using the LP recourse strategy. We 118 solved the linear program DLP with CPLEX 9.0 with default settings on a 3.2 GHz CPU with 2GB RAM. The results are summarized in table 4.6 and Fig. 4.8. Table 4.6: Average percentage of the unmet demand over the total demand for DP and CCP models (lognormal distribution) (%) DL200 DL200 DL250 DL250 DL300 DL300 DL400 DL400 DL500 DL500 DL600 DL600 (40%) (40%) (50%) (50%) (60%) (60%) (80%) (80%) (100%) (100%) (120%) (120%) DP CCP DP CCP DP CCP DP CCP DP CCP DP CCP sup350 (70%) 31.79 35.07 31.32 31.94 31.03 30.91 30.57 30.68 30.78 30.63 31.04 30.63 sup400 (80%) 24.28 30.53 23.20 24.85 22.69 22.09 21.83 21.01 22.13 20.73 21.87 20.70 sup450 (90%) 19.67 28.75 17.48 21.00 16.35 15.63 14.59 12.16 13.98 11.21 14.45 11.00 sup500 (100%) 17.85 28.07 14.78 19.70 13.23 13.28 11.11 6.51 10.15 4.17 9.36 3.34 sup600 (120%) 17.22 28.38 14.37 19.69 11.57 13.00 10.54 5.53 8.55 2.55 7.78 1.44 Figure 4.8: Simulation comparison for the DP and CCP models. The results suggest that, when the deadline is very tight (between 40% to 60% of the base route length), the deterministic routes outperform the CCP routes. There is no benefit for utilizing a conservative routing strategy under tight deadlines since the longer less risky trips (arcs) in most generated instances are longer than the deadlines. In this case, the short-and-risky routes generated by the DP model at least have a chance in some of the generated instances of arriving 119 on-time especially for the demand points in the beginning of the route. In the situation where there is limited supply (e.g., when the total supply is only 70% of the total demand), the DP and CCP routes yield roughly the same amount of unmet demand even under a more relaxed deadline. Because of the insufficient supplies at the depot, the extended deadline cannot provide a better coverage. Under moderate deadlines and supply quantities (when the percentage of the total supply over the total demand and the percentage of the deadline over the base route length both are between 80% and 120%), we can conclude that CCP routes outperform DP routes 2% to 6% in terms of the percentage of the unmet demand over the total demand, according to our simulation experiments. As the deadline and the total supply are relaxed, the percentage of the unmet demand is decreasing for both DP and CCP routes; at the same time, the advantage of CCP routes over DP routes gradually increases. The simulation results on the percentage of the unmet demand from the CCP routes are approaching to its lower bound due to the shortage of the supply, with an increasing deadline. We believe that it is because of, under the relaxed enough deadline and supply quantities, the conservative nature of the CCP model shows its merits. It leads to similar number of nodes in different routes in the CCP model. In contrast, the DP routes are more prone to have an uneven number of nodes between different paths leading to a higher chance of having unmet demand when they are tested over different realization instances. 4.6.3.2 The Role of the Safety Factor Another experiment we conduct is to fix the deadline and the total supply at the depot and vary the confidence level, which has a one-to-one correspondence to the κ T and κ D values in constraints (4.20) and (4.21). We use the same value for both κ T and κ D . For one specific network topology 120 and demand distribution (a problem scenario), we fix the deadline at 80% of the base route length and the total supply at 90% of the total demand. We generate different sets of CCP routes with different safety factors, and run the same simulation as described above. We plot the percentage of the unmet demand over the total demand in Fig. 4.9. Since we are using a lognormal distribution, the deterministic case (mean value) lies at e μ 0 + σ 02 2 . The corresponding safety factor is bigger than 0.5 as compared with exactly 0.5 for the normal distribution. The value of the safety factor for the deterministic case in the lognormal distribution is heavily dependent on both the mean and the variance of the distribution. For the specific problem scenario we are testing, the estimated safety factor is 0.75 by averaging the variance value of all the edges in the graph in a statistical sense. It is shown as a dot on the plot. The point pinned by a square gives the best average unmet demand ratio, which is also the value we used for the CCP model for the comparison purpose in section 4.6.3.1. It corresponds to the confidence level at 95% and the safety factor as 1.65. The plot in Fig. 4.9 infers that with an increasing “contingency stuffing” (a bigger safety factor) starting from the deterministic point, it will provide better coverage. However, if we are too conservative as further increasing the safety factor after the “best CCP point”, the quality of the result deteriorates. Under a given deadline, if the safety factor is too big such that all the short and risky arcs have been “stretched” to overpass the deadline, then the route planning based on those “over-stretched” arcs is not very meaningful in providing any guidance since every node would miss its deadline. Hence it is crucial for planning with the “optimal” safety factor, which is neither too opportunistic nor too conservative. As we discussed earlier, the robust optimization, which aims to optimize the worst-case sce- nario, also shares the same problem structure and computational complexity as the deterministic and the CCP formulations under a bounded uncertainty set or a given set of scenarios. Recently, 121 different groups of researchers are interested in develop more elaborated uncertainty sets in order to address the issue of over-conservative in worst case models while maintaining the computa- tional tractability. Also several works have established the link between the chance-constrained technique and robust optimization (Klopfenstein and Nace, 2006; Calafiore and Campi, 2005; Erdogan and Iyengar, 2006; Chen et.al., 2006). In our work, the only difference between the deterministic counterpart of the robust model and the CCP model are the different safety factors used by each model. In chance-constrained optimization, the constraints are enforced up to a pre- specified level of probability. Robust optimization seeks a solution which satisfies all possible constraint instances. It achieves this by enforcing that the strictest version of every constraint is satisfied by substituting the largest uncertainty parameter value in each constraint. This maximal value for the uncertainty parameter has a one-to-one correspondence with a safety factor value in the chance-constrained formulation and lies somewhere on the plot in Fig. 4.9. Figure 4.9: Simulation comparison for different safety factors. 122 4.6.4 Comparison of Recourse Strategies In this experiment, we follow the simulation process in section 4.6.1 to evaluate and compare three recourse strategies discussed in section 4.4 by allowing split-delivery in the a priori route generation. The result is presented in table 4. The routes obtained in the planning stage are from the chance-constrained model. The number in each grid is the percentage of the unmet demand quantity over the total demand level. For every test case, the re-plan column gives the best result (lowest total unmet demand), since it regenerates both the routes and the commodity flows at the same time, which provides the largest searching space. We note that the solution for the re- planned routes may not necessarily be optimal since we used the Tabu heuristic to identify them. Nevertheless, the re-planned routes outperforms the other two strategies. However, with different realization parameters, this approach may generate significantly different routes. It loses the advantage of the familiarity of preplanned routes. On the other hand, both the LP and knapsack recourse strategies serve for training and preparation purposes. Since the knapsack recourse strategy provides extra flexibility by allowing to skip low-demand and/or far-away nodes in a route, it has a larger search space than the LP strategy. The LP can be solved to optimality using CPLEX and the knapsack strategy can be solved to a near optimal solution by the approach we presented in section 4.5.2. In our numerical simulation, the knapsack recourse strategy always outperforms the LP recourse. When the problem has a tight total supply level these 3 strategies perform similarly, regardless of the deadline. All are very close to the trivial lower bound – the shortage of the supply compared with the total demand. This shows that these models did their best to send out all the available supply within the given deadline in such conditions. On the other hand, when there are tight 123 Table 4.7: Unmet Demand Percentage Comparison Between 3 Recourse Strategies for 50 Cus- tomers and 10 Vehicles. DL 200 DL 250 DL 300 replan LP Knapsack replan LP Knapsack replan LP Knapsack SUP 350 30.3% 30.7% 30.3% 30.3% 30.5% 30.3% 30.3% 30.3% 30.3% SUP 400 20.4% 22.1% 20.4% 20.4% 21.2% 20.4% 20.4% 20.8% 20.4% SUP 450 10.5% 16.7% 11.3% 10.5% 13.7% 10.8% 10.5% 12.5% 10.6% SUP 500 1.9% 15.2% 8.1% 1.9% 10.5% 5.7% 1.9% 8.0% 4.2% SUP 600 0.0% 15.2% 6.8% 0.0% 10.3% 5.6% 0.0% 7.5% 3.2% DL 400 DL 500 DL 600 replan LP Knapsack replan LP Knapsack replan LP Knapsack SUP 350 30.3% 30.3% 30.3% 30.3% 30.3% 30.3% 30.3% 30.3% 30.3% SUP 400 20.4% 20.4% 20.4% 20.4% 20.4% 20.4% 20.4% 20.4% 20.4% SUP 450 10.5% 11.2% 10.5% 10.5% 10.9% 10.5% 10.5% 10.7% 10.5% SUP 500 1.9% 5.4% 3.2% 1.9% 4.1% 2.7% 1.9% 3.5% 2.5% SUP 600 0.0% 4.5% 1.5% 0.0% 2.9% 1.0% 0.0% 2.1% 0.8% deadlines, as the total supply level goes up, the gaps between the re-plan strategy and the knapsack strategy increase. This increase, but more pronounced, is also observed when comparing the re- plan and the LP strategy. This extra solution search space gives the re-plan and knapsack model more chances to reach a better result. In summary, the knapsack recourse strategy provides a nice trade-off between maintaining the familiarity of the preplanned routes and an efficient solution and quick solution times. Indeed, it is the most efficient recourse strategy that uses the preplanned routes and the knapsack approximation algorithm proposed allows to obtain solutions quickly. 124 Chapter 5 Conclusion and Future Work 5.1 Conclusion In this research, we studied the mathematical models on the distribution process of medical sup- plies to mitigate the impacts and consequences of large-scale emergencies. We analyzed the characteristics of the large-scale emergency and emphasized the special requirements addressed to the inventory management problem on the massive perishable supplies and the routing prob- lems in response to large-scale emergencies. For the perishable inventory management system, trivial extension of a regular EMQ model is adequate when the required minimum inventory is not significant compared with the amount consumed by the regular market demand rate within the shelf-life. Since it can be timely and completely depleted and refreshed by the regular market demand. However, when we consider the VMIs in the SNS for the large-scale emergency setting, the minimum inventory requirement is in a scale which is comparable with the total market consumption within the shelf-life. A more sophisticated inventory management strategy is required to provide the fresh and massive 125 stockpile throughout the time horizon in the system. Hence in this work, we modeled the perish- able inventory management problem with a minimum inventory volume constraint as a modified Economic Manufacturing Quantity (EMQ) model. We discussed the policies and assumptions adopted in this model from both the regular perishable inventory management context and the special constraints on the minimum stock size and maximum inventory cycle length enforced by the large-scale emergency response context. Different possible scenarios are discussed and the calculation on the total cost and boundary conditions for each scenario is presented. The total cost is decomposed into four components: inventory holding costs, fixed ordering costs, purchasing costs and salvage costs. With the aid of inventory plots, we formulated the problem to mini- mize the total relevant cost w.r.t. the production batch size as an unconstrained non-continuous non-differentiable optimization problem. We proved the existence of the local as well as global minimum of the total cost with respect to the production batch size. Hence, an exact solution procedure is proposed and its complexity is proved to be polynomial. We estimated the param- eters in the modified EMQ model for a potential anthrax attack scenario from various sources and used them to compare the proposed model with a standard model to show a significant cost saving of running our system as around 33 million US dollars per year, which saved about 30% of the cost to a standard model. We performed sensitivity analysis on some government controlled parameters in the system and observed that at a given profitability level of the firm, there are tradeoffs between the less amount paid by the government to firms (either by reducing the I min requirement or by reducing the unit price p gov ) with the higher flexibility the government allows to firms (longer time before the expiration to salvage the pills). We decomposed the routing problem into two stages: planning and operational stages. We proposed a chance constrained model to minimize the unmet demand for the planning stage and 126 three different recourse strategies for the operational stage. Unique features and requirement from large-scale emergencies on this problem have been analyzed and embedded into our model. We demonstrated that in the planning stage, a combined routing with profit and regular stochastic routing model is adequate to address both the high uncertainty associated with the emergency scenarios and the minimum unmet demand as the primary objective. Solution approaches for both stages were presented. A tabu search heuristic is proposed to solve the planning stage model. A linear programming solver and a multi-choice knapsack solver is applied to solve the models in the operational stage. We illustrated the effectiveness of the models and solution approaches by computational experiments and concluded that the first stage chance-constrained model combined with the second stage knapsack recourse strategy were most effective in solving our problem at hand. The influence of the safety factor for the planning stage model at different values is also illustrated by simulation results. The link between the deterministic model, the chance constrained model and the robust optimization model by varying the value of the safety factor were established. We demonstrated that it is crucial to plan with a “optimal” safety factor value, which is neither too conservative nor too opportunistic. 5.2 Future Work For the inventory management problem, our current model assumed an uniform unit price on the items sold to the regular market and uniform unit salvage value on the items disposed at the end of each inventory cycle regardless of their age. An interest in the future research is to combine the inventory model with a revenue model to address the potential economic impact with more sophisticated issuing, pricing and disposing strategies which incorporate the age distribution of 127 the stockpile. Another interesting direction in the future work is to extend the deterministic demand rate to a stochastic one. Since it is more realistic to assume that the demand is random and it follows a certain probabilistic distribution and address it with some stochastic analysis techniques. For the routing problem, from the modeling perspective, we note that in classical stochastic programming models with recourse the first stage solution is informed by the outcome of the recourse stage. We did not pursue this approach for two reasons: in a large scale emergency situation the uncertainty is so significant that it is difficult to represent it with scenarios accurately. The second reason is because such a stochastic model with recourse would be considerably more difficult to solve. Nevertheless, a comparison of our preplanned solution against the solution from a stochastic model with recourse is an interesting topic of future research. From the solution approach perspective, we are interested in developing a tighter lower bound with short deadline, which does not only help to evaluate the effectiveness of our heuristic result, but could also facilitate to eliminate those infeasible routes in the earlier stage to improve the efficiency of the heuristic. By further investigating the bounding techniques specially geared to our model, an exact solution approach, such as branch and cut, can be developed to solve the planning stage routing model. 128 Bibliography [1] P. Abad, “Optimal pricing and lot-sizing under conditions of perishability and partial back- ordering,” Management Science, vol. 42, no. 8, pp. 1093–1104, Aug. 1996. [2] ——, “Optimal pricing and lot-sizing under conditions of perishability, finite production and partial backordering and lost sale,” European Journal of Operational Research, vol. 144, no. 3, pp. 677–685, Feb. 2003. [3] T. Abdemaguid and M. M. Dessouky, “Genetic algorithm approach to the integrated inventory-distribution problem,” International Journal of Production Research, vol. 44, pp. 4445–4464, 2006. [4] A. Alexandrescu, Modern C++ Design – Genertic Programming and Design Pattern Ap- plied. Addison-Wesley, 2001. [5] C. Archetti, A. Hertz, and M. Speranza, “A tabu search algorithm for the split delivery vehicle routing problem,” Transportation Science, vol. 40, pp. 64–73, 2006. [6] C. Archetti, A. Hertz, and M. G. Speranza, “Metaheurisitics for the team orienteering problem,” Journal of Heuristics, vol. 13, pp. 49–76, 2007. [7] F. Baita, W. Ukovich, R. Pesenti, and D. Favaretto, “Dynamic routing-and-inventory prob- lems: a review,” Transportation Science, vol. 32, no. 8, pp. 585–598, 1998. [8] A. Ben-Tal and A. Nemirovski, “Robust convex optimization,” Mathematics of Operations Research, vol. 23, no. 4, pp. 769–805, 1998. [9] R. W. Bent and P. V . Hentenryck, “Scenario-based planning for partially dynamic vehicle routing with stochastic customers,” Operations Research, vol. 52, no. 6, pp. 977–987, 2004. [10] D. Bertsimas, “Probabilistic combinational optimization problems,” Ph.D. dissertation, Operation Research Center, Massachusetts Institute of Technology, Cambridge, MA, 1988. [11] ——, “A vehicle routing problem with stochastic demand,” Operations Research, vol. 40, no. 3, pp. 574–585, May 1992. [12] D. Bertsimas and G. V . Ryzin, “A stochastic and dynamic vehicle routing problem in the euclidean plane,” Operations Research, vol. 39, no. 4, pp. 601–615, 1991. [13] ——, “Stochastic and dynamic vehicle routing in the euclidean plane with multiple capac- itated vehicles,” Operations Research, vol. 41, no. 1, pp. 60–76, 1993. 129 [14] D. Bertsimas and M. Sim, “Robust discrete optimization and network flows,” Mathemati- cal Programming, vol. 98, pp. 49–71, 2003. [15] D. Bertsimas and D. Simchi-levi, “A new generation of vehicle routing research: robust algorithms, addressing uncertainty,” Operations Research, vol. 44, no. 2, pp. 286–304, Mar. 1996. [16] N. Bianchessi and G. Righini, “Heuristic algorithms for the vehicle routing problem with simultaneous pick-up and delivery,” Computers and Operations Research, vol. 34, pp. 578–594, 2006. [17] J. Branke, M. Middendorf, G. Noeth, and M. Dessouky, “Waiting strategies for dynamic vehicle routing,” Transportation Science, vol. 39, pp. 298–312, 2005. [18] A. Burnetas and C. Smith, “Adaptive ordering and pricing for perishable products,” Oper- ations Research, vol. 48, no. 3, pp. 436–443, May-Jun. 2004. [19] S. E. Butt and T. M. Cavalier, “A heuristic for the multiple tour maximum collection prob- lem,” Computers and Operations Research, vol. 21, no. 1, pp. 101–111, 1994. [20] S. E. Butt and D. M. Ryan, “A optimal solution procedure for the multiple tour maxi- mum collection problem using column generation,” Computers and Operations Research, vol. 26, no. 4, pp. 427–441, 1999. [21] G. Calafiore and M. Campi, “Uncertain convex programs: randomized solutions and con- fidence levels,” Mathematical Programming, vol. 102, pp. 25–46, 2005. [22] R. Carraway, T. Morin, and H. Moskowitz, “Generalized dynamic programming for stochastic combinatorial optimization,” Operations Research, vol. 37, no. 5, pp. 819–829, 1989. [23] I.-M. Chao, B. Golden, and E. Wasil, “A fast and effective heuristic for the orienteering problem,” European Journal of Operational Research, vol. 88, no. 3, pp. 475–489, 1996. [24] ——, “The team orienteering problem,” European Journal of Operational Research, vol. 88, no. 3, pp. 464–474, 1996. [25] D. Chazan and S. Gal, “A markovian model for a perishable product inventory,” Manage- ment Science, vol. 23, no. 5, pp. 512–521, Jan. 1977. [26] X. Chen, M. Sim, and P. Sun, “A robust optimization perspective on stochas- tic programming,” Optimization Online, 2005, http://www.optimization- online.org/DB HTML/2005/06/1152.html. [27] Y . H. Chun, “Optimal pricing and ordering policies for perishable commodities,” European Journal of Operational Research, vol. 144, no. 1, pp. 62–82, Jan. 2003. [28] G. Clarke and J. Wright, “Scheduling of vehicles from a central depot to a number of delivery points,” Operations Research, vol. 12, pp. 568–581, 1964. 130 [29] M. Dror, “Modeling vehicle routing with uncertain demands as a stochastic program: prop- erties of the corresponding solution,” European Journal of Operational Research, vol. 64, pp. 432–441, 1993. [30] M. Dror and M. Ball, “Inventory/routing: reduction from an annual to a short-period prob- lem,” Naval Research Logistics, vol. 34, no. 6, pp. 891–905, 1987. [31] M. Dror, G. Laporte, and P. Trudeau, “Vehicle routing with stochastic demands: properties and solution framework,” Transportation Science, vol. 23, no. 3, pp. 166–176, Aug. 1989. [32] ——, “Vehicle routing with split deliveries,” Discrete Applied Mathematics, vol. 50, pp. 239–254, 1994. [33] M. Dror and P. Trudeau, “Stochastic vehicle routing with modified savings algorithm,” European Journal of Operational Research, vol. 23, pp. 228–235, 1986. [34] ——, “Savings by split delivery routing,” Transportation Science, vol. 23, pp. 141–145, 1989. [35] ——, “Split delivery routing,” Naval Research Logistics, vol. 37, pp. 383–402, 1990. [36] J. Dunn, “Scheduling underway replenishment as a generalized orienteering problem,” Master’s thesis, Naval Postgraduate School, Monterey, California, USA, 1992. [37] E. Erdoˇ gan and G. Iyengar, “Ambiguous chance constrained problems and robust opti- mization,” Mathematical Programming, vol. 107, pp. 37–61, 2006. [38] D. Feillet, P. Dejax, and M. Gendreau, “Traveling salesman problems with profits,” Trans- portation Science, vol. 39, no. 2, pp. 188–205, 2005. [39] M. Ferguson, V . Jayaraman, and G. Souza, “Note: An application of the eoq model with nonlinear holding cost to inventory management of perishables,” European Journal of Op- erational Research, vol. 180, no. 1, pp. 485–490, Jul. 2007. [40] M. Fischetti, J. S. Gonzalez, and P. Toth, “A branch-and-cut algorithm for the generalized traveling salesman problem,” Operations Research, vol. 45, no. 3, pp. 378–394, 1997. [41] ——, “Solving the orienteering problem through branch-and-cut,” INFORMS Journal of Computing, vol. 10, no. 2, pp. 133–148, 1998. [42] B. Fries, “Optimal ordering policy for a perishable commodity with fixed lifetime,” Oper- ations Research, vol. 23, no. 1, pp. 46–61, Jan.-Feb. 1975. [43] F. Fumero and C. Vercellis, “Synchronized development of production, inventory, and distribution schedules,” Transportation Science, vol. 33, no. 3, pp. 330–340, 1999. [44] G. Gallego and G. van Ryzin, “Optimal dynamic pricing of inventories with stochastic demand over finite horizons,” Management Science, vol. 40, no. 8, pp. 999–1020, Aug. 1994. 131 [45] E. Gamma, R. Helm, R. Johnson, and J. Vlissides, Design Patterns – Elements of Reusable Object-Oriented Software. Addison-Wesley, 1995. [46] M. Gendreau, A. Hertz, and G. Laporte, “A tabu search heuristic for the vehicle routing problem,” Management Science, vol. 40, pp. 1276–1290, 1994. [47] M. Gendreau, G. Laporte, and J. Y . Potvin, The Vehicle Routing Problem. SIAM Mono- graphs on Discrete Mathematics and Applications, SIAM Publishing, 2002, ch. Meta- heuristics for the capacitated VRP, pp. 129–154. [48] M. Gendreau, G. Laporte, and R. Seguin, “An exact algorithm for the vehicle routing problem with stochastic demands and customers,” Transportation Science, vol. 29, no. 2, pp. 143–155, May 1995. [49] ——, “Stochastic vehicle routing,” European Journal of Operational Research, vol. 88, pp. 3–12, 1996. [50] ——, “A tabu search heuristic for the vehicle routing problem with stochastic demands and customers,” Operations Research, vol. 44, no. 3, pp. 469–477, May 1996. [51] M. Gendreau, G. Laporte, and F. Semet, “A branch-and-cut algorithm for the undirected selective traveling salesman problem,” Networks, vol. 32, pp. 263–273, 1998. [52] ——, “A tabu search heuristic for the undirected selective traveling salesman problem,” European Journal of Operational Research, vol. 106, pp. 539–545, 1998. [53] B. Giri and K. Chaudhuri, “Determinisitic models of perishable inventory with stock- dependent demand rate and nonlinear holding cost,” European Journal of Operational Research, vol. 105, no. 3, pp. 467–474, Mar. 1998. [54] G. Glover and M. Laguna, Tabu Search. Kluwer, Boston, MA, 1997. [55] A. Goel and V . Gruhn, “Large neighborhood search for rich vrp with multiple pickup and delivery locations,” in MEC-VNS: 18th Mini Euro Conference on VNS, 2005. [56] C.-H. Goh, B. Greenberg, and H. Matsuo, “Two-stage perishable inventory models,” Man- agement Science, vol. 39, no. 5, pp. 633–649, May 1993. [57] B. Golden, L. Levy, and R. V ohra, “The orienteering problem,” Naval research logistics, vol. 34, no. 3, pp. 307–318, 1987. [58] B. Golden, Q. Wang, and L. Liu, “A multifaceted heuristic for the orienteering problem,” Naval research logistics, vol. 35, pp. 359–366, 1988. [59] S. Graves, “The application of queuing theory to continuous perishable inventory systems,” Management Science, vol. 28, no. 4, pp. 400–406, Apr. 1982. [60] E. Hadjiconstantinou and D. Roberts, The Vehicle Routing Problem. SIAM Monographs on Discrete Mathematics and Applications, SIAM Publishing, 2002, ch. Routing under uncertainty: an application in the scheduling of field service engineers, pp. 331–352. 132 [61] M. Hayes and J. Norman, “Dynamic programming in orienteering: route choice and sitting of controls,” Journal of the Operational Research Society, vol. 35, pp. 791–796, 1984. [62] V . N. Hsu, “Dynamic economic lot size model with perishable inventory,” Management Science, vol. 46, no. 8, pp. 1159–1169, Aug. 2000. [63] V . N. Hsu and T. Lowe, “Dynamic economic lot size models with period-pair-dependent backorder and inventory costs,” Operations Research, vol. 49, no. 2, pp. 316–321, Mar.- Apr. 2001. [64] H. Hwang, “Food distribution model for famine relief,” Computers and Industrial Engi- neering, vol. 37, no. 1, pp. 335–338, 1999. [65] T. Ilhan, S. Iravani, and M. Daskin, “The orienteering problem with stochastic profits,” IIE Transactions, vol. 40, pp. 406–421, 2008. [66] P. Jaillet, “Stochastic routing problem,” in Stochastics in Combinatorial Optimization, G. Andreatta, F. Mason, and P. Serafini, Eds. World Scientific, New Jersey, 1987, pp. 178–186. [67] P. Jaillet and A. Odoni, Vehicle Routing: Methods and Studies. North-Holland, Amster- dam, 1988, ch. The Probabilistic Vehicle Routing Problem. [68] A. J´ ez´ equel, “Probabilistic vehicle routing problems,” Master’s thesis, Department of Civil Engineering, Massachusetts Institute of Technology, 1985. [69] H. Jia, “Models and solution approaches for facility location of medical supplies for large- scale emergencies,” Ph.D. dissertation, University of Southern California, 2006. [70] H. Jula, M. M. Dessouky, and P. Ioannou, “Truck route planning in non-stationary stochas- tic networks with time-windows at customer locations,” IEEE Transactions on Intelligent Transportation Systems, 2005, to appear. [71] E. Kao, “A preference order dynamic program for a stochastic travelling salesman prob- lem,” Operations Research, vol. 26, pp. 1033–1045, 1978. [72] H. Katagiri and H. Ishii, “Fuzzy inventory problems for perishable commodities,” Euro- pean Journal of Operational Research, vol. 138, no. 3, pp. 545–553, May 2002. [73] O. Klopfenstein and D. Nace, “A robust approach to the chance-constrained knapsack problem,” Optimization Online, 2006, http://www.optimization- online.org/DB HTML/2006/03/1341.html. [74] R. Kopach, B. Balcioglu, and M. Carter, “Tutorial on constructing a red blood cell in- ventory management system with two demand rates,” European Journal of Operational Research, vol. 185, no. 3, pp. 1051–1059, Mar. 2006. [75] V . Lambert, G. Laporte, and F. Louveaux, “Designing collection routes through bank branches,” Computers and Operations Research, vol. 20, pp. 783–791, 1993. 133 [76] G. Laporte, “The vehicle routing problem: an overview of exact and approximate algo- rithms,” European Journal of Operational Research, vol. 59, pp. 345–358, 1992. [77] G. Laporte, F. Laporte, and H. Mercure, “Models and exact solutions for a class of stochas- tic location-routing problems,” European Journal of Operational Research, vol. 39, pp. 71–78, 1989. [78] G. Laporte, F. Louveaux, and H. Mercure, “The vehicle routing problem with stochastic travel times,” Transportation Science, vol. 26, no. 3, pp. 161–170, Aug. 1992. [79] G. Laporte and S. Martello, “The selective traveling salesman problem,” Discrete applied mathematics, vol. 26, pp. 193–207, 1990. [80] R. C. Larson, The McGraw-Hill Handbook of Homeland Security. The McGraw-Hill Companies, 2005, ch. Decision Models for Emergency Response Planning. [81] R. C. Larson, M. Metzger, and M. Cahn, “Responding to emergencies: Lessons learned and the need for analysis,” Interfaces, vol. 36, no. 6, pp. 486–501, 2006. [82] A. Leifer and M. Rosenwein, “Strong linear programming relaxations for the orienteering problem,” European Journal of Operational Research, vol. 73, pp. 517–523, 1994. [83] Y .-C. Liang and A. Smith, “An ant colony approach to the orienteering problem,” Depart- ment of industrial and systems engineering, Auburn University, Auburn, AL, Tech. Rep., 2001. [84] L. Liu and Z. Lian, “(s, s) continuous review models for products with fixed lifetimes,” Operations Research, vol. 47, no. 1, pp. 150–158, Jan.-Feb. 1999. [85] N. Mladenovic and P. Hansen, “Variable neighborhood search,” Computers and Operations Research, vol. 24, pp. 1097–1100, 1997. [86] S. Nahmias, “Perishable inventory theory: A review,” Operations Research, vol. 30, no. 4, pp. 680–708, Jul.-Aug. 1982. [87] P. Nandakumar and T. Morton, “Near myopic heuristics for the fixed-life perishability problem,” Management Science, vol. 39, no. 12, pp. 1490–1498, Dec. 1993. [88] I. H. Osman, “Metastrategy simulated annealing and tabu search algorithms for the vehicle routing problem,” Annals of Operations Research, vol. 41, pp. 421–451, 1993. [89] J. D. Papastavrou, “A stochastic and dynamic routing policy using branching processes with state dependent immigration,” European Journal of Operational Research, vol. 95, pp. 167–177, 1996. [90] W. Pierskalla and C. Roach, “Optimal issuing policies for perishable inventory,” Manage- ment Science, vol. 18, no. 11, pp. 603–614, Jul. 1972. [91] D. Pisinger, “A minimal algorithm for the multiple-choice knapsack problem,” European Journal of Operational Research, vol. 83, pp. 394–410, 1995. 134 [92] F. Raafat, “Survey of literature on continuously deteriorating inventory models,” The Jour- nal of the Operational Research Society, vol. 42, no. 1, pp. 27–37, Jan. 1991. [93] R. Ramesh, Y .-S. Yoon, and M. Karwan, “An efficient four-phase heuristic for the gen- eralized orienteering problem,” Computers and Operations Research, vol. 18, no. 2, pp. 151–165, 1991. [94] ——, “An optimal algorithm for the orienteering tour problem,” ORSA Journal of Com- puting, vol. 4, no. 2, pp. 155–165, 1992. [95] N. Ravichandran, “Stochastic analysis of a continuous review perishable inventory system with positive lead time and poisson demand,” European Journal of Operational Research, vol. 84, no. 2, pp. 444–457, Jul. 1995. [96] D. Rosenfield, “Disposal of excess inventory,” Operations Research, vol. 37, no. 3, pp. 404–409, May-Jun. 1989. [97] ——, “Optimality of myopic policies in disposing excess inventory,” Operations Research, vol. 40, no. 4, pp. 800–803, Jul.-Aug. 1992. [98] N. Secomandi, “A rollout policy for the vehicle routing problem with stochastic demands,” Operations Research, vol. 49, no. 5, pp. 796–802, 2001. [99] M. Sobel and R. Zhang, “Inventory policies for systems with stochastic and deterministic demand,” Operations Research, vol. 49, no. 1, pp. 157–162, Jan.-Feb. 2001. [100] D. Socolar and A. Sager, “Cipro prices expose windfall profits,” Health Services Depart- ment, Boston University School of Public Health, June 2002. [101] W. Stewart and B. Golden, “Stochastic vehicle routing: a comprehensive approach,” Euro- pean Journal of Operational Research, vol. 14, pp. 371–385, 1983. [102] I. Sungur, F. Ord´ o˜ nez, and M. M. Dessouky, “A robust optimization approach for the ca- pacitated vehicle routing problem with demand uncertainty,” University of Southern Cali- fornia, Working paper 2006-06, 2006. [103] M. R. Swihart and J. D. Papastavrou, “A stochastic and dynamic model for the single- vehicle pick-up and delivery problem,” European Journal of Operational Research, vol. 114, pp. 447–464, 1999. [104] E. D. Taillard, P. Badeau, M. Gendreau, F. Guertin, and J. Y . Potvin, “A tabu search heuristic for the vehicle routing problem with soft time windows,” Transportation Science, vol. 31, pp. 170–186, 1997. [105] H. Tang and E. Miller-Hooks, “Algorithms for a stochastic selective travelling salesperson problem,” Journal of the Operational Research Society, vol. 56, pp. 439–452, 2005. [106] ——, “A tabu search heuristic for the team orienteering problem,” Computers and Opera- tions Research, vol. 32, pp. 1379–1407, 2005. 135 [107] F. M. Tasgetiren and A. E. Smith, “A genetic algorithm for the orienteering problem,” in Proceeding of 2000 Congress Evolutionary Computation, 2000, pp. 1190–1195. [108] E. Tekin, U. Gurler, and E. Berk, “Age-based vs. stock level control policies for a perish- able inventory system,” European Journal of Operational Research, vol. 134, no. 2, pp. 309–329, Oct. 2001. [109] F. Tillman, “The multiple terminal delivery problem with probabilistic demands,” Trans- portation Science, vol. 3, pp. 192–204, 1969. [110] P. Toth and D. Vigo, “The granular tabu search (and its application to the vehicle routing problem),” DEIS, Universit` a di Bologna, Italy, Tech. Rep. OR/98/9, 1998. [111] ——, The Vehicle Routing Problem. SIAM Monographs on Discrete Mathematics and Applications, SIAM Publishing, 2002. [112] P. Trudeau and M. Dror, “Stochastic inventory routing: route design with stockouts and route failures,” JTS, vol. 26, no. 3, pp. 171–184, 1992. [113] T. Tsiligirides, “Heuristic methods applied to orienteering,” Journal of Operational Re- search Society, vol. 35, no. 9, pp. 797–809, 1984. [114] T. Vaughan, “A model of the perishable inventory system with reference to consumer- realized product expiration,” The Journal of the Operational Research Society, vol. 45, no. 5, pp. 519–528, May 1994. [115] S. Viswanathan and K. Mathur, “Integrating routing and inventory decisions in one ware- house multicustomer multiproduct distribution systems,” JMS, vol. 43, no. 3, pp. 294–312, 1997. [116] H. M. Wagner and T. M. Whitin, “Dynamic version of the economic lot size model,” Management Science, vol. 5, no. 1, pp. 89–96, Oct. 1958. [117] Q. Wang, X. Sun, B. Golden, and J. Jia, “Using artificial neural networks to solve the orienteering problem,” Annuls of Operational Research, vol. 61, pp. 111–120, 1995. [118] C. Waters, “Vehicle-scheduling problems with uncertainty and omitted customers,” Jour- nal of the Operations Research Society, vol. 40, pp. 1099–1108, 1989. [119] H. Weiss, “Optimal ordering policies for continuous review perishable inventory models,” Operations Research, vol. 28, no. 2, pp. 365–374, Mar.-Apr. 1980. [120] J. A. G. Willard, “Vehicle routing using r-optimal tabu search,” Master’s thesis, The Man- agement School, Imperial College, London, 1989. [121] C. Williams and B. E. Patuwo, “A perishable inventory model with positive order lead times,” European Journal of Operational Research, vol. 116, no. 2, pp. 352–373, Jul. 1999. 136 [122] ——, “Analysis of the effect of various unit costs on the optimal incoming quantity in a perishable inventory model,” European Journal of Operational Research, vol. 156, no. 1, pp. 140–147, Jul. 2004. [123] T. Wu, “Optimization models for underway replenishment of a dispersed carrier battle group,” Master’s thesis, Naval Postgraduate School, Monterey, California, USA, 1992. [124] J. Xu and J. P. Kelly, “A network flow-based tabu search heuristic for the vehicle routing problem,” Transportation Science, vol. 30, pp. 379–393, 1996. [125] H. Yan and T. Cheng, “Optimal production stopping and restarting times for an eoq model with deteriorating items,” The Journal of the Operational Research Society, vol. 49, pp. 1288–1295, 1998. 137
Abstract (if available)
Abstract
Rapid and efficient wide-scale distribution of medical supplies plays a critical role in assuring the effectiveness in managing the risks of large-scale emergencies such as a bio-terrorism attack. Important issues in the design of such an efficient distribution network involve deciding how to route distribution vehicles and how to manage these inventories. The high uncertainty, limited supply and overwhelming demand associated with a large-scale emergency may result in significant unmet demand, which is a direct representation of the most undesirable consequence -- loss of life. Solving appropriate vehicle routing and inventory management problems in a coordinated manner can ensure the design of a logistic network capable of efficiently distributing medical supplies to decrease the potential fatalities in responding to a large-scale emergency. In this work, we develop models and solution approaches to solve a perishable inventory management problem and a vehicle routing problem in the context of emergency response to minimize unmet demand.
Linked assets
University of Southern California Dissertations and Theses
Conceptually similar
PDF
New approaches for routing courier delivery services
PDF
Vehicle routing and resource allocation for health care under uncertainty
PDF
Models and solution approaches for facility location of medical supplies for large-scale emergencies
PDF
Dynamic programming-based algorithms and heuristics for routing problems
PDF
Routing for ridesharing
PDF
Models and algorithms for energy efficient wireless sensor networks
PDF
Models and algorithms for pricing and routing in ride-sharing
PDF
Models and algorithms for the freight movement problem in drayage operations
PDF
Evaluating city logistics using two-level location routing modeling and SCPM simulation
PDF
Train routing and timetabling algorithms for general networks
PDF
Personalized Pareto-improving pricing-and-routing schemes with preference learning for optimum freight routing
PDF
Computational validation of stochastic programming models and applications
PDF
High-performance distributed computing techniques for wireless IoT and connected vehicle systems
Asset Metadata
Creator
Shen, Zhihong
(author)
Core Title
Routing and inventory model for emergency response to minimize unmet demand
School
Viterbi School of Engineering
Degree
Doctor of Philosophy
Degree Program
Industrial and Systems Engineering
Publication Date
08/05/2008
Defense Date
03/28/2008
Publisher
University of Southern California
(original),
University of Southern California. Libraries
(digital)
Tag
emergency response,mixed integer programming,nonlinear programming,OAI-PMH Harvest,perishable inventory management,stochastic vehicle routing with profit
Language
English
Contributor
Electronically uploaded by the author
(provenance)
Advisor
Dessouky, Maged M. (
committee chair
), Ordonez, Fernando (
committee chair
), Koenig, Sven (
committee member
)
Creator Email
shenz@usc.edu,zhihong.shen@gmail.com
Permanent Link (DOI)
https://doi.org/10.25549/usctheses-m1543
Unique identifier
UC191252
Identifier
etd-Shen-2319 (filename),usctheses-m40 (legacy collection record id),usctheses-c127-106172 (legacy record id),usctheses-m1543 (legacy record id)
Legacy Identifier
etd-Shen-2319.pdf
Dmrecord
106172
Document Type
Dissertation
Rights
Shen, Zhihong
Type
texts
Source
University of Southern California
(contributing entity),
University of Southern California Dissertations and Theses
(collection)
Repository Name
Libraries, University of Southern California
Repository Location
Los Angeles, California
Repository Email
cisadmin@lib.usc.edu
Tags
emergency response
mixed integer programming
nonlinear programming
perishable inventory management
stochastic vehicle routing with profit