Inventory and scheduling problems in supply chain management

c 2013 Selim Bora ALL RIGHTS RESERVED INVENTORY AND SCHEDULING PROBLEMS IN SUPPLY CHAIN MANAGEMENT by SELIM BORA A dissertation submitted to the Graduate School-New Brunswick Rutgers, The State University of New Jersey in partial fulfillment of the requirements for the degree of Doctor of Philosophy Graduate Program in Operations Research written under the direction of Endre Boros and approved by New Brunswick, New Jersey January, 2013 ABSTRACT OF THE DISSERTATION Inventory and Scheduling Problems in Supply Chain Management by SELIM BORA Dissertation Director: Endre Boros This dissertation deals with real life inventory and supply chain management problems. Three chapters are dedicated to different problems, two of them dealing with vessel scheduling and container management, where as the last one deal with inventory management of critical supplies at hospitals, during time of surging demand. We propose heuristics and exact approaches for these problems that are both efficient in terms of time, and accurate enough to our needs. Even though inventory and supply chain management field is a popular field of research, these problems have not been addressed yet, and known results cannot be directly applied to any of them. ii Acknowledgements This dissertation is the fruit of years of research that I have done with many people, whose contribution in numerous ways to the research and the making of the thesis deserves special mention. In the first place, I would like to state my gratitude to Endre Boros for his supervision, advice, and guidance through out the work. I gratefully acknowledge Lei Lei for her advice, supervision, and crucial contribution to the problems dealt in Chapter 3 and Chapter 4 of this dissertation. She was the one to get me started on my dissertation, and allowed me to gain experience in various aspects of Supply Chain Management field. Many thanks to Wei Xiong and Tayfur Altiok for introducing me to the problem dealt in Chapter 2 of this thesis, and to Melike Baykal-Gursoy for her supervision after Tayfur Altiok’s unfortunate passing away. In addition, I would like to thank many people who have been with me along this road. My mother and my father, for supporting me in every way possible... Tayfur Altiok, for showing me the way... Melike Baykal-Gursoy, for being there when needed... Endre Boros, for being kind, helpful and patient to me, and not giving up on me... iii Contents Abstract ii Acknowledgements iii List of Figures vi List of Tables vii 1 Introduction 1 2 Outbreak Detection Using Hidden Markov Models and Inventory Management 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 A Known Disease . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.1 The Algorithm for the Optimal Target Levels . . . . . . . . . . . . 2.2.2 Scenarios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Joint Disease Identification and Inventory Management . . . . . . . . . . 2.3.1 The EM Algorithm for Disease Identification . . . . . . . . . . . . 2.3.2 Influence of Lead Time . . . . . . . . . . . . . . . . . . . . . . . . 2.3.3 Deterministic Model for Forecasting . . . . . . . . . . . . . . . . . 2.3.4 Comparison of EM Algorithm and the Deterministic Model Based on Revision Frequency . . . . . . . . . . . . . . . . . . . . . . . . . 3 A Case of Container-Vessel Scheduling Problem 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . 3.2 Problem Definition . . . . . . . . . . . . . . . . . . 3.2.1 Backward Heuristic . . . . . . . . . . . . . . 3.2.2 Greedy Heuristic . . . . . . . . . . . . . . . 3.2.3 Improved Greedy Heuristic . . . . . . . . . 3.2.4 Bender’s Type Decomposition Approach . . 3.3 Computational Results . . . . . . . . . . . . . . . . 3.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . iv . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 5 9 13 15 18 18 20 22 23 27 27 31 38 39 40 40 41 43 4 Container Vessel Scheduling with Time Windows 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . 4.2 Literature . . . . . . . . . . . . . . . . . . . . . . . 4.3 Problem Definition . . . . . . . . . . . . . . . . . . 4.4 Minimum Cost Flow Network Formulation . . . . . 4.5 Minimum Cost Flow Heuristic . . . . . . . . . . . . 4.5.1 Stage 1 . . . . . . . . . . . . . . . . . . . . 4.5.2 Stage 2 . . . . . . . . . . . . . . . . . . . . 4.5.3 Stage 3 . . . . . . . . . . . . . . . . . . . . 4.6 Decomposition Approach . . . . . . . . . . . . . . 4.7 Results . . . . . . . . . . . . . . . . . . . . . . . . . 5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 45 47 51 57 58 59 59 60 60 64 65 v List of Figures 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 2.10 2.11 2.12 2.13 Schematic relationship between the five subgroups in the model . . . . . Inventory control model used in the hospital model . . . . . . . . . . . . Case 1, Daily demand, total cost per day and the optimized target inventory levels . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Case 2, Daily demand, total cost per day and the optimized target inventory levels . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Case 3, Daily demand, total cost per day and the optimized target inventory levels . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . No Lead Time, Daily demand and Target Inventory Level based on hidden Markov model and EM algorithm . . . . . . . . . . . . . . . . . . . . . . 30 Day Lead Time, Daily demand and Target Inventory Level based on hidden Markov model and EM algorithm . . . . . . . . . . . . . . . . . 60 Day Lead Time, Daily demand and Target Inventory Level based on hidden Markov model and EM algorithm . . . . . . . . . . . . . . . . . Daily Renewals, Daily demand and Target Inventory Level based on Hidden Markov and Deterministic Approach, 60 Day Lead Time . . . . . . Weekly Renewals, Daily demand and Target Inventory Level based on Hidden Markov and Deterministic Approach, 60 Day Lead Time . . . . Bi-Weekly Renewals, Daily demand and Target Inventory Level based on Hidden Markov and Deterministic Approach, 60 Day Lead Time . . . . Monthly Renewals, Daily demand and Target Inventory Level based on Hidden Markov and Deterministic Approach, 60 Day Lead Time . . . . Hidden Markov and Deterministic Combined, Daily demand and Target Inventory Level based on Hidden Markov and Deterministic Approach . . 10 . 11 . 16 . 17 . 17 . 21 . 21 . 22 . 24 . 24 . 25 . 25 . 26 3.1 An example of the network for |V1 | = 2, |V2 | = 1, |V3 | = 2, |N | = 3, and |T | = 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 4.1 An example of the network for |T | = 7, |P | = 2, Tp,q = 1 ∀p, q ∈ P, and time windows for delivery are indicated by brackets around the nodes . . 59 vi List of Tables 2.1 2.2 3.1 3.2 3.3 3.4 3.5 4.1 4.2 For |T |=2, given that h < p, 2nd and 4th scenario need comparison, and we see that no intermediate value will yield a better cost . . . . . . . . . . 14 For |T |=3, given that h < p, This time 2nd and 8th scenarios need comparison, yielding us the same result . . . . . . . . . . . . . . . . . . . . . . 14 |N |=3, |T |=4, vessel type same, number of vessels in each period shown, as well as the optimal objective function value for each case, indicating that even a much simpler version of the original problem is not submodular |N |=10, |T |=10, 3 different vessel types, number of vessels each method solves for vary for XpressMP, Backward and Greedy as 10, 100, 1 in respect. Runs are terminated after 2 hours or when 0.1% gap from the best bound is reached . . . . . . . . . . . . . . . . . . . . . . . . . . . . . |N |=10, |T |=10, 3 different vessel types, number of vessels each method solves for vary for XpressMP and Improved Greedy as 10 and 1 in respect. All runs are terminated after 2 hours or when 0.1% gap from the best bound is reached . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . |N |=10, |T |=10, 3 different vessel types, number of vessels each method solves for vary for XpressMP and Improved Greedy as 10 and 1 in respect. All runs are terminated after 2 hours or when 0.1% percent gap from the best bound is reached . . . . . . . . . . . . . . . . . . . . . . . . . . . . . |N |=15, |T |=10, 3 different vessel types, number of vessels each method solves for vary for XpressMP and Improved Greedy as 10 and 1 in respect. All runs are terminated after 2 hours or when 0.1% percent gap from the best bound is reached . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 42 42 42 43 |N | = 20, |V | = 40, |T | = 30, 0.1% gap termination condition or 2 hrs for both XpressMP and Decomposition . . . . . . . . . . . . . . . . . . . . . . 64 |N | = 40, |V | = 80, |T | = 30, 0.1% gap termination condition or 2 hours for both XpressMP and Decomposition . . . . . . . . . . . . . . . . . . . . 64 vii DEDICATION PAGE Dedicated to Dr. Tayfur Altiok. . . viii 1 Chapter 1 Introduction The importance of supply chain management (SCM) cannot be denied. Global markets are expanding beyond borders and redefining the way demand and supplies are managed. Therefore, management of a network of all business processes and activities involving procurement of raw materials, manufacturing and distribution of finished goods, whether intercontinental or international is getting more and more complex every day. In addition to this, the cost margins are increasing even further. The goal of this work is to tackle some problems existing in SCM literature, which are all relevant to some real-life problems, and develop techniques to give insight for future decisions. The research carried out in the first problem addresses the issue of managing critical medical supply inventories in hospital under surge scenarios. It is based on an approach 2 that optimizes inventory control parameters in hospital settings under scenarios such as the pandemic flu with surging demand for medical supplies. In the decade before the 2009-10 influenza pandemic caused by the new H1N1 virus (pH1N1), the spread of the much more lethal H5N1 ”avian” influenza in Asia and parts of Africa raised concerns about the potentially devastating impact of a severe global influenza outbreak (Salomon and Webster [45] and Chan [8]). In response, most developed countries and many private corporations made considerable investments over the last decade in the purchase of antiviral medications (AVMs) to treat those infected with influenza during a pandemic. To date, most of the literature has addressed either pandemic mitigation or preservation of healthcare workforce capacity during the peak of an outbreak. However, to the best of our knowledge, there is no evidence-based research available in the medical literature that can guide healthcare facilities to establish sufficient medical supply in order to maintain adequate surge capacity for flu patients or other pandemic. In the inventory management and supply chain literature, most existing models make the assumption of independence and time-homogeneity of the demand for medical supply. However, during a pandemic, demand is uncertain and definitely non-stationary due to surge dynamics. It grows exponentially in the early part of the pandemic and will decline when there are less people that can be infected by the disease (herd immunity). The non-stationary demand pattern raises a unique problem for inventory management since it becomes difficult to decide on when to order and how much to order of each of the critical medical supplies. This is also a significant research issue in the inventory management literature. Non-stationarity of the demand gives rise to complexities that make the problem mathematically and practically intractable. 3 The next set of problems are related to container-vessel scheduling. In the existing literature, it is possible to find many publications on vessel scheduling. Even though, the basic assumptions are very similar, the problems that have been dealt with are all different either because there’s an additional assumption, or the solution approach is different. In this research, we try to come up with a general models which covers many of these assumptions, and work with special cases, which appear quite common as subproblems in the literature. We propose heuristics and exact formulations based on a Bender’s type decomposition approach. Issues with container-vessel scheduling are vast. One of the issues is related to empty containers piling up at the ports. Mallon and Magaddino [37] suggest that empty containers dwell in container terminals for the longest on average, and indicate that therefore they are the main source of congestion at and around marine ports. Unfortunately, there is not much research on reuse or redistribution of empty containers, and we hope to contribute to this area, more precisely to the problem of distributing empty containers as they become available at terminals. The other vessel scheduling problem we deal with originates from the demand-supply coordination problem encountered by a major U.S. corporation in petrochemical industry such as oil and gasoline. The cost to produce and deliver gasoline products to the market consists of three major components: the transportation cost of crude oil to refiners, the operation cost of refinery processing, and the cost of marketing and distribution (Riley [40]). An oil company typically operates many tens of refineries, with several million 4 barrels of crude oil per day and several billion dollars on crude transportation per year. As the U.S. retail gasoline prices continue to rapidly elevate, effectively coordinating the demand and supply of gasoline products has therefore become even more crucial to oil companies (Cheng and Duran [10]). In this particular study, the distribution scheduling problem encountered in this process is very complicated due to the involvement of heterogeneous vessels (e.g., in terms of their loading capacities, discharging and berthing times, and operating costs) and the fact that each vessel has multi-level of loading capacities such that a load beyond the normal/base capacity will result in an extra overload cost. In addition, the inventories of the refinery depots must be managed, as well as, satisfying demand somehow. These inventory management and distribution processes are often referred to an integrated scheduling problem. For both of the problems, we start by stating the full version of the problem, which is hard to tackle, then show that an easier formulation is possible, which allows us to apply Bender’s type decomposition approach, that has a subproblem similar to a network flow problem, assuring integrality, and therefore computational efficiency. This formulation allows us to get the optimal solution, and we show that this approach performs better than a state of the art solver, such as XpressMP. In addition, we propose greedy heuristics for the same type of problems, which yield solutions much faster than both XpressMP and the Bender’s type decomposition approach, with a relatively higher, but still acceptable relative error from the best optimal solution. 5 Chapter 2 Outbreak Detection Using Hidden Markov Models and Inventory Management 2.1 Introduction The studies that have been done related to this problem can be categorized in four groups. The classical problem of finding optimal inventory levels for base-stock or (s,S) policies with backlogging has been dealt by Veinott Jr [48], where the model has the following features: There is a general demand process with no stationarity or independence assumptions, partial or complete backlogging of unfilled demand, a fixed nonnegative delivery lag (which may be positive only under complete backlogging), a nonstationary linear ordering cost, a nonstationary holding and shortage cost function, discounting of 6 future costs, and nonstationary restrictions like budget and storage limitations. The objective is to choose an ordering policy that minimizes the expected discounted costs over an infinite time horizon. They give these conditions to ensure that the base stock ordering policy is optimal and that the base stock levels in each period are easy to calculate. Karlin [29] formulates a dynamic inventory model in which the demand distributions may change from period to period. He shows that the optimal policy at each stage is characterized by a single critical number which also could vary in successive periods. He develops the dependence of the critical numbers as a function of stochastic ordering amongst distributions is developed under various conditions. Even though most of his studies are conducted under the assumption of linear purchasing cost, the possibility of convex purchasing cost is allowed. Also, Karlin and Scarf [30] partially characterize the structure of optimal policy where the lead time is 1 period, which in summary is as follows: Optimal ordering quantity is a decreasing function of inventory on hand at the beginning of a period and zero, outside a specified interval. Moreover, the rate of decrease (as a function of the inventory on hand) is strictly less than 1. With the additional assumption that demands are exponentially distributed, they also present a steady-state analysis of the dynamics of lost-sales systems that use base-stock policies. Morton [38] has extended this study to lost-sales model with deterministic lead times. He has shown that optimal ordering policy is a function of entire pipeline with following characteristics: 1) There’s a compact region around origin (all components of the pipeline vector are 0) such that order quantity is strictly greater than 0 if and only if 7 pipeline vector is in this region. 2) Order quantity decreases at a rate strictly between zero and one with respect to each component of the pipeline vector. 3) The rate of decrease in the order quantity per component is higher for components in the pipeline that are scheduled to arrive later in time. Song and Zipkin [47] deals with a similar problem, where the demand forms Markovmodulated Poisson process. The lead time however can be deterministic or stochastic. Next, Graves [24] considers an adaptive base-stock policy for a single-item inventory system, where the demand process is nonstationary. He then shows how the single-item model extends to a multi-stage, or supply-chain context. Some work in the literature incorporate forecast information. In Schoenmeyr and Graves [46], they use a case study with real data to demonstrate that there are significant benefits from the inclusion of the forecast process when determining the optimal safety stocks. They also conduct a computational experiment to explore how the placement and size of the safety stocks depend on the nature of the forecast evolution process. Also, Graves [24] deal with the inventory management of single item for base-stock policy, where the demand process is nonstationary and they use an autoregressive integrated moving average process for it. Buzacott and Shanthikumar [7] study the robustness of the parameters safety lead time and safety stock, and conclude that safety time is usually only preferable to safety stock when it is possible to make accurate forecasts of future required shipments over the lead time, otherwise safety stock is more robust in coping with changes in customer requirements in the lead time or with fluctuations in forecasts 8 of lead time demand for material requirements planning controlled productions systems. Chang [9] deals with the interchangeability of safety stocks and safety time. Last but not the least, Levi et al. [34] extend the model from backlogging to lost-sales, and devise approximation algorithm with computational efficiency which admits a worst-case performance guarantee. In our research however, we also wish to be able to deal with any disease, whether an existing one, or a totally new disease. Le Strat and Carrat [32] use a hidden Markov model for monitoring diseases, and show that they can accurately define the states, the distribution parameters according to the states, and the transition probabilities in between states. We try to extend this idea for monitoring, and currently identifying states of a disease and when an outbreak is going on. The research carried out in this project addresses the issue of inventory management that basically deals with how much of each unit to maintain on hand using and approach that optimizes inventory control system parameters in hospital settings under scenarios such as the pandemic flu with surging demand for medical supplies. This is a critical issue since too much of medical supply inventory can cost a lot while short inventories will not satisfy the demand. Finding the right amount of stock is always a challenge in inventory management due to uncertainties involved in demand. We then extend this study for case of encountering a disease for the first time. Inventory control with unknown demand is a field not limited only to hospitals of course. 9 Zheng and Federgruen [51] develop an efficient algorithm for reorder point/reorder quantity policies, where they allow backlogging and deal with poisson distribution. However, the demand distribution has fixed parameters. Geary et al. [22] deal with demand amplification(i.e. bullwhip), and present 10 published causes of bullwhip that could be avoided by reconstructing the supply chain. This study, makes our point about accurate inventory management of hospital supplies stronger. The contribution of this research is in providing the proof of a concept to hospital management that tools such as the one developed in this research will be instrumental in managing inventories once deployed in full scale. This research aims to provide the proof of concept for developing an approach to handle inventory management under surge scenarios, whether new or old, in hospitals. The objective is to help hospitals in making decisions regarding maintaining stocks of medical supplies. This is achieved by developing a formal procedure to help effectively control inventories of critical medical supplies and minimize inventory management costs while maintaining an acceptable customer service level in pandemic-like scenarios. 2.2 A Known Disease In essence, using a system of ordinary differential equations, a compartmental SEIR model (Anderson et al. [1]) is developed describing the transmission dynamics of a pandemic in a large population. Individuals in a population are divided into standard modeling compartments such as susceptible (S), exposed (E), infectious (I), recovered 10 (R), and disease-induced dead (D). As shown in Figure 2.1, susceptible individuals may become infected, incubate the infection, progress to become fully contagious, and finally either recover and develop immunity from the disease or die. Figure 2.1: Schematic relationship between the five subgroups in the model The arrows that connect the boxed subgroups represent movement of individuals. Susceptible (S) will first become exposed (E). Exposed (E) will become infectious (I) after an incubation period. Infectious (I) can either recover (R) after a recovery period or die (D). Following the schematic representation in Figure 2.1, the disease spread and control dynamics can be described by the following differential equations for a population of size N (Anderson et al. [1]). s = −βIS/N (2.1) E = βIS/N − κE (2.2) I = κE − (γ + δ)I (2.3) R = γI (2.4) D = δI (2.5) Then, the SEIR model is incorporated into a virtual hospital simulator using ARENA 11 simulation tool to identify daily random demand over the course of the pandemic. A dynamic programming model is introduced to obtain the optimum target inventories of the three products mentioned earlier in the virtual hospital model. The optimization will be cost minimization based on inventory holding cost, shortage cost and the cost of changing the target inventory level. The model is to be incorporated into the simulation model to manage the inventories of the virtual hospital introduced earlier. Finally a dynamic programming algorithm is used to optimize the target inventory levels as well as reorder points for the duration of the pandemic period. Clearly, the values of these parameters vary over time due varying demand over the pandemic duration. We have employed a periodic review model where the inventory of each product is monitored daily and a replenishment order is placed every time the inventory is observed below its target level. The order is place at 13:00 in the afternoon and is assumed to arrive by the end of the day, presumably before the beginning of the next day. The inventory management concept used here is presented in Figure 2.2. Figure 2.2: Inventory control model used in the hospital model Next, we present the optimization model. We use the disease parameters and related differential equations (2.1) in Arena, replicating many times the flow of sick people to 12 the hospital, and estimate the parameters related with the demand distribution function of our choice, which in this study has been assumed to be Poisson. Karlin and Scarf [30] shows that the random variable, y, measuring stock available at the next period, has the following stationary distribution: . F (y) = φ(y)[1 − φ(x̄ − y)] , 1 − [1 − φ(y)][1 − φ(x̄ − y)] (2.6) where x̄ is the target inventory level and φ(y) is the integral of density function. Accordingly, the expected quantity of unsatisfied demand is given by: . E(penalty) = Z x̄ ϕ(ξ)dξ 0 Z x̄ (ξ − y)f (y)dy + 0 Z ∞ ϕ(ξ)dξ x̄ Z ∞ (ξ − y)f (y)dy (2.7) x̄ Also, the amount on hand at the end of a day is given by: . E(handling) = Z 0 x̄ Z y (y − ξ)ϕ(ξ)f (y)dξdy, (2.8) 0 where ϕ(ξ) is the daily demand density function. Now we have the expected on-hand inventory level and the expected shortage level, we are ready to evaluate the cost minimizing objective function for feasible levels of target inventory of the product. The following algorithm will be used to generate the optimal target inventory level of the medical supply we have under consideration. 13 2.2.1 The Algorithm for the Optimal Target Levels 1. Given the disease parameters, run Arena to obtain the demand distribution function and its parameters for every day 2. Compute x̄t for every day t, minimizing the expected total cost per day, that is E(penalty) ∗ p + E(handling) ∗ h 3. Use dynamic programming to determine when the target inventory level must be changed. The dynamic programming formulation is done in such a way that the inventory level is the state, and the days are the stages in the model. So, the possible states (target inventory values) are x̄1 , x̄2 , . . . , x̄T . Also, let us define T Ct (I¯t ) to be the total cost up to day t, where the current inventory level is I¯t . Then, the following dynamic programming recursion can be established: . T Ct (I¯t ) = minx̄t P (x̄t ) + T Ct+1 (I¯t+1 ), (2.9) . P (x̄t ) = (It − Dt )+ ∗ h + (Dt − It )+ ∗ p + CC(x̄t ), (2.10) . ¯ CC(x̄t ) = M, x̄t 6= It−1 (2.11) where and where We do not need to consider other values due to the following result: Assuming that daily demand will equal to precomputed inventory levels, i.e. Dt = x̄t , no other target 14 values needs consideration for dynamic programming. We show this result by induction, skipping the case where T = 1. For T = 2, we are facing the following scenarios, where x̄1 < x̄3 < x̄2 and costs on individual days are given in Table 2.1: Scenario x̄1 x̄2 x̄1 , x̄2 x̄3 Day 1 0 h(x̄2 − x̄1 ) 0 + CC(x̄2 ) h(x̄3 − x̄1 ) Day 2 p(x̄2 − x̄1 ) 0 0 p(x̄2 − x̄3 ) Table 2.1: For |T |=2, given that h < p, 2nd and 4th scenario need comparison, and we see that no intermediate value will yield a better cost For T = 3, where x̄1 < x̄3 < x̄2 . An x̄4 value that is intermediate will be considered yielding following scenarios given in Table 2.2: Scenario x̄1 x̄2 x̄3 x̄1 , x̄2 x̄1 , x̄3 x̄2 , x̄3 x̄1 , x̄2 , x̄3 x̄4 > x̄3 x̄4 < x̄3 Day 1 0 h(x̄2 − x̄1 ) h(x̄3 − x̄1 ) 0 + CC(x̄2 ) 0 + CC(x̄3 ) h(x̄2 − x̄1 ) 0 + CC(x̄2 ) h(x̄4 − x̄1 ) h(x̄4 − x̄1 ) Day 2 p(x̄2 − x̄1 ) 0 p(x̄2 − x̄3 ) 0 p(x̄2 − x̄3 ) 0 + CC(x̄3 ) 0 + CC(x̄3 ) p(x̄2 − x̄4 ) p(x̄2 − x̄4 ) Day 3 p(x̄3 − x̄1 ) h(x̄2 − x̄3 ) 0 h(x̄2 − x̄3 ) 0 0 0 h(x̄4 − x̄3 ) p(x̄4 − x̄3 ) Table 2.2: For |T |=3, given that h < p, This time 2nd and 8th scenarios need comparison, yielding us the same result The relationship between consecutive days can either be reduced to case |T | = 2 or |T | = 3, which proves the result. 15 Thus, we are able to obtain optimal target inventory values for the medical supply under consideration for days 1 through T of the surge period. This algorithm can easily be programmed in an operational environment and implemented for each product on a daily basis. Below we provide some scenario analysis. 2.2.2 Scenarios In this section, we present three case studies of the virtual hospital for a product, to show the impact of holding cost, shortage cost and target level change cost. In all of these cases, we assume the following epidemiological parameter settings for our differential equations: Incubation Period: 1.9 days Infectious Period: 4.1 days Population Size: 100,000 Case 1: Holding Cost: 8/unit/day Shortage Cost: 16/unit Cost of Changing Target Level: 1,000 The daily demand, and the optimized target inventory levels for the product and the total cost per day are shown in Figure 2.3. In this case, we have kept the values of the cost parameters at relatively low levels. This ensures frequent changes in target levels 16 to adapt the changing environment due to surge. Figure 2.3: Case 1, Daily demand, total cost per day and the optimized target inventory levels Case 2: Holding Cost: 8/unit/day Shortage Cost: 16/unit Cost of Changing Target Level: 6,000 In this case, the cost of changing target inventory levels is increased by six fold. The daily demand, and the optimized target inventory levels for the product and the total cost per day are shown in Figure 2.4. Observe the impact of large target changing costs in the form of much lesser changes of target values. Compare the blue line in Figure 2.4 to Figure 2.3. Unit shortage cost is relatively higher than the unit holding cost and therefore the optimization model tries to keep the target levels close to the demand and yet changes them much less infrequently as compared to Case 1. Case 3: In this case, we will look into the impact of high unit holding costs. The cost parameters 17 Figure 2.4: Case 2, Daily demand, total cost per day and the optimized target inventory levels are given below: Holding Cost: 16/unit/day Shortage Cost: 8/unit Cost of Changing Target Level: 6,000 Figure 2.5: Case 3, Daily demand, total cost per day and the optimized target inventory levels Observe the impact of large unit holding cost which discourage holding high inventories and push the target levels to low values as observed in Figure 2.5 for the same product. 18 2.3 Joint Disease Identification and Inventory Management We have shown that our algorithm is able to deal with inventory management of any disease during a surge, as long as, we have some information about the disease. What if we were dealing with a new disease, where no preexisting data exists? Next, we try to tackle the question of identifying a disease, different states of it(i.e. regular, surge), how frequent each occurs, and parameters related to the distributions of these states. For that, we build on the idea of Le Strat and Carrat [32], we approach the problem by using two different forecasting tools, and then applying an optimization model to decide on the inventory control policies. 2.3.1 The EM Algorithm for Disease Identification We use EM algorithm to estimate the parameters of the unknown disease. We assume the random variable associated with daily incoming patients is exponential. We assume our hidden Markov model to have 2 states, one for regular times, and one for during time of surge. The EM algorithm computes the steady-state and transition probabilities between the states, and distribution parameters associated with the states, with the following steps: 1. Let yt , t = 1, ..., n are a realization of the stochastic process Y = (Yt ; t = 1, ...n, ). The idea is to associate each Yt an unobserved random variable St that determines 19 the conditional distribution of Yt ; if St = j, then the conditional distribution of Yt has density fj (yt ; θj ), where fj belongs to a given parameterized family, and θj are parameters to be estimated. We assume that the unobserved sequence St follows an m-state homogeneous Markov chain of order 1 with stationary transition probabilities αij = P (St = j|St−1 = i), i, j = 1, ..., m. 2. Let ψ = (θ1 , ..., θm , α11 , ..., αmm ) denote the complete parameter vector to be estimated. The estimation of the parameters by the EM algorithm starts by initializing the parameters (0) (0) (0) (0) (α1 , ..., αm ), (θ1 , ..., θm ) as well as the transition probability matrix (0) [αjk ]1≤j,k≤m 3. Let vjk (t) = P (St−1 = j, St = k|y1 , ..., yn ) and uj (t) = P (St = j|y1 , ..., yn ) be the conditional expectations given the observations and current parameters estimates. The values of vjk and uj can be obtained by computing: aj (t) = f (y1 , ..., yt , St = j) and bj (t) = f (yt+1 , ..., Yn |St = j) with the recursive forward-backward formulae. 4. At each iteration ω, ω = 1, ..., I, the E-Step is followed by the M-Step. 5. E-Step: Specifically, for j = 1, .., m: (0) (ω) aωj (1) = αj fj1 (y1 ; θj ) aωj (t) = (ω) k=1 αk (t Pm (ω) (ω) − 1)αkj fjt (yt ; θj ), t = 2, ..., n bωj (n) = 1 bωj (t) = (ω) (ω) (ω) k=1 αjk fkt+1 (yt+1 ; θk )bk (t Pm + 1), t = n − 1, ..., 1 The likelihood is calculated by Ln (ψ) = (ω) (ω) j=1 aj bj (t), Pm t = 1, ..., n 20 Finally, for t = 1, ..., n, (ω) (ω) uωj (t) = aj (t)bj (t) Pm (ω) l=1 al (n) ω (t) = vjk αjk fk (yt ;θk aj (t−1)bk (t) Pm (ω) l=1 al (n) (ω) (ω) (ω) (ω) 6. M-Step: The transition probabilities, stationary probabilities and θj : (ω+1) αjk = (ω+1) = αj 2.3.2 Pn (ω) t=2 vjk (t) Pm (ω) t=2 l=1 vjl (t) Pn (ω) t=1 ujk (t) Pm Pn (ω) k=1 t=1 uk (t) Pn Influence of Lead Time Once the parameters are found, the future demand is estimated based on state parameters, and distribution parameters. Based on future predictions, inventory target levels that would minimize the total shortage and holding cost are computed, and then dynamic programming is applied to decide when to change the inventory target level to another quantity, just like in the previous section. However, there is a shortcoming with this approach, which is displayed below. Basically, the hidden Markov model approach requires a learning period, and until it has been achieved, the predictions towards future are not very reliable. The Figures 2.6, 2.7 and 2.8 show how this approach is affected, based on the lead time we provide the algorithm to learn the disease. As it can be seen from Figures 2.6, 2.7 and 2.8, the hidden Markov model for disease identification performs really good, when we allow up to 60 days for learning. However, in real life, especially during a time of surge, especially when dealing with an unknown 21 Figure 2.6: No Lead Time, Daily demand and Target Inventory Level based on hidden Markov model and EM algorithm Figure 2.7: 30 Day Lead Time, Daily demand and Target Inventory Level based on hidden Markov model and EM algorithm disease, there is no such time to sit idly by, so for quicker response, we also incorporated a deterministic approach similar to that of Schoenmeyr and Graves [46]. 22 Figure 2.8: 60 Day Lead Time, Daily demand and Target Inventory Level based on hidden Markov model and EM algorithm 2.3.3 Deterministic Model for Forecasting We’ll require the following variables: V ariable Group It T ype Explanation t∈T Z+ inventory on hand at the beginning of day t Pt t ∈ T Z+ number of supplies ordered on day t Lt t ∈ T Z+ base stock level for day t xt t∈T {0, 1} xt = 1 if base stock level is changed to another value on day t 23 Then, we get the following model: min X ht (It − Dt ) + t∈T X pt (Dt − It ) + t∈T X CCt xt (2.12a) t∈T s.t.I0 = D0 (2.12b) It = It−1 − Dt + Pt ∀ t ∈ T (2.12c) It ≥ Lt t ∈ T (2.12d) xt ≤ Lt − Lt−1 t ∈ T (2.12e) (2.12a) is the objective function, consisting of three parts, holding, shortage, and inventory target level change cost respectively. (2.12b) and (2.12c) are inventory balance constraints, where as (2.12d) is to make sure that our beginning of the day inventory is always at least as much as our base stock level, and (2.12e) is to control the binary variable associated with inventory target level change. We replace Dt for future periods with the most recent realized value of Dt , as revision takes place. 2.3.4 Comparison of EM Algorithm and the Deterministic Model Based on Revision Frequency The Figures 2.6, 2.7 and 2.8 show that deterministic approach requires much less learning time for predictions to become accurate, however it is possible to see that hidden Markov model predictions eventually start yielding a better result. Of course, one important parameter for deterministic approach, as well as the hidden Markov approach, is the policy revision frequency, as our knowledge extends, and thus our predictions could 24 be better. Next set of Figures 2.9, 2.10, 2.11 and 2.12 display how well both hidden Markov approach, and deterministic approach does, based on the revision frequency. Note that, the duration set for predictions to be recomputed also define when a possible inventory level change decision is made. Figure 2.9: Daily Renewals, Daily demand and Target Inventory Level based on Hidden Markov and Deterministic Approach, 60 Day Lead Time Figure 2.10: Weekly Renewals, Daily demand and Target Inventory Level based on Hidden Markov and Deterministic Approach, 60 Day Lead Time 25 Figure 2.11: Bi-Weekly Renewals, Daily demand and Target Inventory Level based on Hidden Markov and Deterministic Approach, 60 Day Lead Time Figure 2.12: Monthly Renewals, Daily demand and Target Inventory Level based on Hidden Markov and Deterministic Approach, 60 Day Lead Time As expected, going over the predictions and inventory level change decisions performed more frequently, yield better results. However, it must be kept in mind that, updating predictions is a computational task, which still requires considerable amount of time. In 26 addition to that, changing the target inventory level in daily life is not such a straightforward task and it brings along additional costs. Based on our results, we have decided that a weekly or bi-weekly are the best options, which is based on the cost of inventory change level, as well as the holding and shortage costs. Because of the fact that a deterministic approach initially performs better, and we don’t have time to learn, especially when facing a new disease, we propose an approach that brings together two ideas listed above. We initially start with the deterministic approach, giving our hidden Markov model enough time to learn the disease(i.e. 60 days in this case), and then switch, so that we can utilize the strengths of both approaches. Figure 2.13 shows the increase in performance. Figure 2.13: Hidden Markov and Deterministic Combined, Daily demand and Target Inventory Level based on Hidden Markov and Deterministic Approach 27 Chapter 3 A Case of Container-Vessel Scheduling Problem 3.1 Introduction Our goal is to minimize the operating costs related to shipping and handling of goods. The fleet size is not fixed, nor an initial amount is set, so one of the tasks we have at hand is to determine the number of vessels that will be used within the planning horizon. Shipping costs can be divided into two categories:1) The fixed cost related to either purchase or lease of a vessel,2) the overloading cost which is incurred if the vessels carry above a certain capacity. There are two more costs that we need to watch out for. Each shipment made to a port may incur a holding or penalty cost based on the demand. If the demand is not met on time, it cannot be satisfied at a later time period, and therefore we need to pay penalty for each unit. Also, if the port is forced to hold 28 some inventory, then a holding cost is charged. In addition to all these cost factors, we also need to consider the fact that each vessel is available for a certain amount of time within a period, and therefore even if a vessel has enough capacity, it may not have enough time to visit all the ports we desire. Optimally solving distribution operations scheduling problem is not an easy task. Previous work related to industrial shipping varies a lot. Here, we focus on the existing results that are closely related to our work. A large summary of works related to various types of vessel scheduling and routing problem can be found in the literature survey by Christiansen et al. [16]. Two more recent surveys can be found more specifically in the area of combined inventory management and routing (Andersson et al. [2]) and on fleet composition and routing (Hoff et al. [26]). Xinlian et al. [50] presents an algorithm which combines the linear programming technique with that of dynamic programming to improve the solution to linear model for fleet planning. Even though their approach is similar, the problem they are dealing with requires demand satisfaction and initial fleet is already given, and the decision is to whether add new vessels to the existing fleet or not. Cho and Perakis [13] presented a better formulation to the original fleet deployment problem proposed by Ronen [43]. In this formulation, just like we do, there is a single loading port, finite number of customer ports, and a finite planning horizon. However, they require the demand to be met, and the fleet size is constant. The costs incurred 29 are due to routes chosen, shipping cargoes, and unloading time. They show that this formulation is better for computational efficiency. Cho and Perakis [12] present a study regarding fleet size and design of optimal liner routes for a container shipping company. The problem is solved by generating a number of candidate routes for the different ships first, and then, the problem is formulated and solved as a linear programming model, where the columns represent the candidate routes. They extend this model to a mixed integer programming model that also considers investment alternatives to expanding fleet capacity. Bendall and Stent [3] also present a model for determining the optimal number of ships and fleet deployment plan. On the other hand, Nicholson and Pullen [39] were the first ones to propose dynamic programming application to ship fleet management. The problem they dealt with was to determine the sequence in which the currently owned ships should be sold and the extent to which charter ships should be taken on. They tackle the problem in two stages. The first stage determines a good priority ordering for selling the ships regardless of the rate at which charter ships are taken on. The second stage uses dynamic programming to determine an optimal level of chartering given the priority replacement order. This first stage priority ordering essentially reduces the dynamic programming calculation from a problem with as many as states as number of ships in fleet to a 1 state variable problem which is computationally manageable by dynamic programming methods. 30 Several authors use benchmark instances to compare the results of different strategies and heuristics. Gheysens et al. [23] define 20 test instances with 12100 nodes for the standard fleet size and mix vehicle routing problem. Wu et al. [49] deals with trucks that vary in capacity and age are utilized over space and time to meet customer demand. Operational decisions (including demand allocation and empty truck repositioning) and tactical decisions (including asset procurements and sales) are explicitly examined in a linear programming model to determine the optimal fleet size and mix. The method uses a time-space network, common to fleet-management problems, but also includes capital cost decisions, wherein assets of different ages carry different costs, as is common to replacement analysis problems. A two-phase solution approach is developed to solve large-scale instances of the problem. Phase I allocates customer demand among assets through Benders decomposition with a demand-shifting algorithm assuring feasibility in each subproblem. Phase II uses the initial bounds and dual variables from Phase I and further improves the solution convergence through the use of Lagrangian relaxation. A network optimization approach has been proposed by Bookbinder and Reece [4], where they formulate a multi-commodity capacitated distribution-planning problem as a non-linear mixed integer programming model, and solve it as a generalized assignment problem within an algorithm for the overall distribution/routing problem based on a Bender’s type decomposition. Lei et al. [33] proposes an approach to a bi-directional flow problem where each iteration starts with a given planning horizon, which is then partitioned into three planning 31 intervals, where each interval consists of consecutive time periods in the given planning horizon. Afterwards, some constraint relaxations are applied to the problem in which all the forward demand and all the backward demand of the time periods in the third planning interval are consolidated into a single forward demand and a single backward demand, which is an idea we use in one of our approaches. Choi et al. [14] focuses on minimizing total tardiness, rather than the operating costs, and the routes for vessels are observed under three different cases, one of them being arbitrary, just like in our problem. Later on, they talk about the other problems in the literature and how their approach is related to them. 3.2 Problem Definition This paper, brings together some of the ideas that were proposed in the literature before. We are given a fleet |V | of container vessels, v ∈ V that distributes the goods from a main distribution center to a number of customer ports over a |T |-period planning horizon. Each vessel has two loading capacities: the regular loading capacity u0v , and the maximum loading capacity umax so that carrying a load beyond u0v will impose an v over loading charge gv0 /unit and carrying a load beyond umax violates the feasibility. In v addition to this limitation, for every vessel there is total available time τv which is used up by the berthing time bv,p at ports which vary depending on vessel type. There are |p| customer ports on the network, each port p ∈ P has a demand, dp,t ≥ 0 in period t ∈ T. For every port, unsatisfied demand are penalized at pp,t /unit based on the unsatisfied 32 demand and no backlogging is allowed. On the other hand, end of period inventory incurs a holding cost of hn /unit. Let cfv denote the fixed cost if the vessel is being dispatched in a period. The problem is finding a feasible vessel dispatching schedule to minimize the total shortage and overage penalty plus the vessel overloading and fixed cost.The minimum cost flow network formulation proposed guarantees optimality when the number of vessels dispatched in every period is known. To define our problem more formally, we define the following set of variables: V ariable Group T ype p∈Pt∈T Z+ Qv,p,t v ∈ V p ∈ P t ∈ T Z+ Sp,t Explanation amount of shortage at port p in period t, amount of supply delivered to port p in period t via vessel v’s regular capacity Ov,p,t v ∈ V p ∈ P t ∈ T amount of supply delivered to port p in period Z+ t via vessel v’s overloading capacity Ip,t v∈Vt∈T Yv,p,t v ∈ V p ∈ P t ∈ T {0, 1} Z+ ending inventory at port p in period t Yv,n,t = 1 if vessel v delivers to port p in period t Zv,t v∈Vt∈T {0, 1} Zv,t = 1 if vessel v is dispatched in period t Based on this, the constraints to the problem will include the following: 33 A vessel must not be carrying anything if it’s not dispatched, nor visiting ports: Yv,p,t ∀ v ∈ V, p ∈ P, t ∈ T Qv,p,t + Ov,p,t ≤ umax v Yv,p,t ≤ Zv,t ∀ v ∈ V, p ∈ P, t ∈ T (3.1a) (3.1b) Vessels dispatched must not be used over their time and regular/maximum capacity: X bv,p Yv,p,t ≤ τv ∀ v ∈ V, t ∈ T (3.2a) Qv,p,t ≤ u0v ∀ v ∈ V, t ∈ T (3.2b) (Qv,p,t + Ov,p,t ) ≤ umax ∀ v ∈ V, t ∈ T v (3.2c) p∈P X p∈P X p∈P The last group of constraints is to help to formulate our objective, which is a compositions of all expenses (penalties, etc.). Vessel dispatching costs: cD = XX cfv Zv,t (3.3a) XX hp Ip,t (3.3b) v∈V t∈T Early arrival penalties: cH = p∈P t∈T 34 Unsatisfied demands’ penalties: cU = XX pp,t Sp,t (3.3c) p∈P t∈T Overloading penalties: cO = X XX gv0 Ov,p,t (3.3d) v∈V p∈P t∈T Then our problem is to minimize cD + cH + cU + cO , subject to the constraints (3.1)-(3.3) and the sign and type restrictions in the definitions of the decision variables. If the dispatching information is already available, i.e. |V1 | vessels for t=1, |V2 | for t=2, . . . , |VT | for t=T, then there becomes no need for the binary variables. In addition, define new variables, xv,k,n,t and rv,k,n,t , which are the normal and over flows shipped by vessel v dispatched in period k for port n to satisfy the demand on period t. Based on this definition, the following can be established: Qv,p,t = Ov,p,t = T X k=t T X xv,t,p,k ∀ v ∈ V, p ∈ P, t ∈ T (3.4a) rv,t,p,k ∀ v ∈ V, p ∈ P, t ∈ T (3.4b) k=t Sp,t = dp,t − X X (xv,k,p,t + rv,k,p,t ) ∀ p ∈ P, t ∈ T (3.4c) (xv,k,p,w + rv,k,p,w ) ∀ p ∈ P, t ∈ T (3.4d) k∈T v∈Vk Ip,t = T X X X k∈T w=t+1 v∈Vk 35 Based on the above assumptions and definitions, we get the following model: min XX pp,t (dp,t − p∈P t∈T + XX X (xv,k,p,t + rv,k,p,t )) + k=1 v∈Vk v∈V t∈T s.t. T X X cfv Zv,t + XX p∈P t∈T X XX v∈V p∈P t∈T hp T t X X X gv0 T X rv,t,p,k k=t (xv,k,p,w + rv,k,p,w ) (3.5a) k=1 w=t+1 v∈Vk bv,p Yv,p,t ≤ τv ∀ v ∈ V, t ∈ T (3.5b) p∈P T X X (xv,t,n,k + rv,t,n,k ) ≤ umax ∀ v ∈ V, t ∈ T v (3.5c) k=t n∈N T X xv,t,n,k + rv,t,n,k ≤ umax Yv,n,t ∀ v ∈ V, p ∈ P, t ∈ T v (3.5d) k=t T XX n∈N k=t t X X xv,t,n,k ≤ u0v ∀ v ∈ V, t ∈ T (3.5e) (xv,k,n,t + rv,k,n,t ) ∀ p ∈ P, t ∈ T (3.5f) k=1 v∈Vk Yv,n,t ≤ Zv,t ∀ v ∈ V, p ∈ P, t ∈ T (3.5g) Objective function is the same except that the last part is now a constant based on vessel dispatching information, i.e. Zv,t values are known. 3.5b and 3.5c assure that normal and over capacity are not exceeded, where as 3.5d prevents shipments for a specific demand to be more than the demand itself, therefore making the first part of the objective function always nonnegative. Lemma 3.2.1. The above problem can be reformulated without the berthing time constraint and solved as a minimum cost flow problem by assuming the knowledge of the number of vessels dispatched in each time period. Proof. First, we construct a dummy source node S, and a dummy sink node F . Associate 36 to each vessel v ∈ Vk , 2 nodes (v, k)P and (v, k)O , one for normal and other for over capacity. These nodes are connected to the source node with 0 and gv0 costs, a lower bound of 0 and an upper bound u0v and umax − u0v respectively. Add another set of |P | v nodes (pp ) for case of shortage at each port with 0 costs, 0 lower bounds and no upper bounds. Next, take care of the ports by adding |P | ∗ |T | nodes denoted (p, t) for each port n at every period t. The arcs between nodes corresponding to vessels and ports incur a holding cost of hp (t − k), has a lower bound of 0 and no upper bound. Also, there will be arcs between shortage nodes, (pp ), and ports, (p, t), where the shortage costs pp,t will be charged. Finally, add arcs between ports and the sink, with a lower and upper bound of dp,t and no cost.This network will have 2|V ||T | + |P | + |P ||T | many nodes, and |V ||P ||T 2 | + |P 2 | many arcs, making minimum cost flow approach practical for problems of reasonable size. An example network is shown in Figure 3.1. Lemma 3.2.2. The objective function values and constraints for both problems above are the same, assuming we guessed the right number of vessels. Proof. First of all, the fixed cost due to vessels for both problems will be the same. ∗ Next, assume x∗v,k,n,t and rv,k,n,t are the optimal flow vectors corresponding to the min- imum cost flow problem. Then, using the equalities corresponding the variables of two problems, the objective function value of the original problem becomes: P P p∈P P v∈V P P p∈P t∈T pp,t (dp,t p∈P P P t∈T hp − 0 t∈T gv P P k∈T P v∈Vk (xv,k,p,t PT k∈T k=t rv,t,p,k PT w=t+1 P + P v∈V + rv,k,p,t ))+ P f t∈T cv Zv,t + v∈Vk (xv,k,p,w + rv,k,p,w ) 37 Figure 3.1: An example of the network for |V1 | = 2, |V2 | = 1, |V3 | = 2, |N | = 3, and |T | = 3 = P p∈P P ∗ t∈T pp,t Sp,t P P v∈V + P v∈V f ∗ t∈T cv Zv,t + P P p∈P p∈P P P 0 ∗ t∈T gv Ov,k,p,t + ∗ t∈T hp Ip,t The first 2 lines of this objective function and the objective function of the minimum cost flow are exactly the same, which only leaves us with the inventory part. The w index is for shipments that are on a future date than current period t, and the k index is taking into account all shipments that have been made up to period t. Therefore, a shipment made on period k for period t will appear in the summation (t − k) many 38 times, allowing us to replace index w with t, remove the summation regarding w, and charge the holding cost as many times as necessary. This shows that both objective function values are the same. As far as the constraints are concerned, first realize that in the original problem, 3.2a is no longer required while berthing times are large enough. Similarly, 3.1a and 3.1b were associated with the fact that dispatching information was not available, so now, they could be dropped as well. 3.2b in the original problem is the same constraint as 3.5d in the reduced problem, and they are both concerned with normal capacity of a vessel. 3.5c and 3.5d of the reduced problem, added together, imply the same restriction on maximum vessel capacity as 3.2c of the original problem. On the other hand, the flow balance constraint in the original problem is taken care of by two means: 1) the new index k for the variables, tells us when shipment was made, so we now whether a shipment is held at inventory or used immediately, 2) in the reduced problem, shipment for a specific demand will not be more than the demand itself, therefore shortage never becomes negative according to the relation between Sn,t and xv,k,n,t , rv,k,n,t . Based on our minimum cost network flow approach, we propose the following two heuristics for no berthing time case: 3.2.1 Backward Heuristic 1. Divide the planning horizon into two groups, primary and secondary, for each port, the new demand is equal to sum of the individual demands in each group, holding 39 cost is the minimum and penalty cost is the maximum of individual penalties. 2. Start with |P | ∗ |T | vessels in that group in total, solve the minimum cost flow problem iterating through all vessel dispatching combinations available for a group. 3. Once the optimal number of vessels required for each group are determined, repeat the procedure of dividing into groups and solving as a minimum cost flow problem for the individual groups. Demand belonging to ports in the other individual group is also added to the demand of the ports in the secondary period of the group under consideration. 4. Once the primary group has only 1 period remaining, optimal number of vessels have been determined for that group, start over. 3.2.2 Greedy Heuristic 1. Start with no vessels assigned to each period. 2. Add a vessel to any period and solve the problem. Remove the vessel, and add to another period, and solve again. Once the best vessel addition has been determined, move on to next vessel addition. 3. Keep determining the best vessel to add until objective function no longer improves. Going back to the original problem with berthing constraints, we propose modified greedy heuristics and a decomposition based exact method. We first introduce the algorithms and then compare them with state of the are integer programming solver, XpressMP. 40 3.2.3 Improved Greedy Heuristic 1. Start with no vessels assigned to each period. 2. For each vessel type, compute maximal subsets of ports such that no further port can be added to a set due to berthing time constraint. 3. For each vessel type, in every period, sort the subset of ports in decreasing order based on P n∈M aximalv dn,t pn,t 4. Next, add any vessel to any period allowing it to only serve the top ranked subset of ports and solve the problem. Try different vessel types, for different periods in the same manner. Determine the best vessel to add to which period. 5. Once a vessel assignment has been determined, update remaining demand and sort subsets of ports accordingly. 6. Keep determining the best vessel to add until objective function no longer improves. 3.2.4 Bender’s Type Decomposition Approach 1. Choose a feasible assignment of ports to a vessel to start with. 2. Solve the master problem to obtain a new objective function value and new port assignments to other vessels. 3. Keeping port assignments fixed, solve the dual problem. 4. STOP, if master and dual objective values are close enough, otherwise go back to the master problem. 41 3.3 Computational Results With the minimum cost flow network formulation proposed, one question that arose was whether it was worth investing more resources into finding a dispatching information, and maybe go from there. However, as can be seen at 3.1, the minimum cost flow problem is not submodular. Case 1 and Case 2 are random dispatches, where as Case Int refers to the scenario where minimum of number of vessels dispatched in each latter case is used, and Case Union refers to of Case 1 and Case 2 are random dispatches, where as Case Int refers to the scenario, where minimum of number of vessels dispatched in each latter case is used, and Case Union refers to the scenario, where maximum of number of vessels dispatched in each latter case is used. Case 1 (3,2,2) (1,2,1) (2,2,1) (1,2,3) (1,3,2) Obj. Func. 170 840 580 350 390 Case 2 (2,1,3) (2,1,3) (2,1,3) (2,1,3) (2,1,3) Obj. Func. 390 390 390 390 390 Case Int (2,1,2) (1,1,1) (2,1,1) (1,1,3) (1,1,2) Obj. Func. 540 1140 880 650 800 Case Union (3,2,3) (2,2,3) (2,2,3) (2,2,3) (2,3,3) Obj. Func. 50 90 90 90 50 Comparison Lower Equal Equal Equal Lower Table 3.1: |N |=3, |T |=4, vessel type same, number of vessels in each period shown, as well as the optimal objective function value for each case, indicating that even a much simpler version of the original problem is not submodular Backward and Greedy Heuristic both performed well, however, as can be seen in Table 3.2, Backward Heuristic always takes shorter, where as Greedy Heuristic performs better by a slight margin. However, it must be kept in mind that Table 3.2 reflects results for the version of the problem with no berthing time constraint. Improved Greedy Heuristic designed to deal with this issue performs a bit slower than the previously mentioned heuristics, but gives 42 Xpress 7200 7200 7200 7200 7200 7200 7200 7200 7200 7200 Time (sec) Backward 412.50 420.40 393.30 332.30 289.90 310.30 345.2 429.50 421.8 306.70 Greedy 435.8 435.20 420.60 390.70 316.80 336.60 347.6 443.30 425.8 321.10 Objective Function Xpress Backward Greedy 7503 7799 7733 8141 8314 8298 8270 8483 8450 7759 7806 7800 7316 7395 7375 8412 8494 8487 8270 8356 8338 7918 8159 8112 7475 7825 7768 7448 7661 7600 Xpress 1.7384 2.1557 0.9988 0.3519 0.1414 3.2309 3.8278 1.1641 1.4574 2.4489 Gap (%) Backward 5.47 4.20 3.48 0.95 1.20 4.16 4.82 4.08 5.87 5.16 Greedy 4.66 4.01 3.11 0.88 0.94 4.09 4.62 3.53 5.17 4.40 Table 3.2: |N |=10, |T |=10, 3 different vessel types, number of vessels each method solves for vary for XpressMP, Backward and Greedy as 10, 100, 1 in respect. Runs are terminated after 2 hours or when 0.1% gap from the best bound is reached good bounds for the solution of the original problem as can be seen in Table 3.3. Time (sec) Xpress I. Greedy 7200 593.9 7200 582.2 7200 635.3 7200 678.2 7200 668.6 7200 594.4 7200 649.3 7200 667.3 7200 641.1 7200 662.1 Objective Function Xpress I. Greedy 7503 7533 8141 8226 8270 8381 7759 7921 7316 7436 8412 8548 8270 8353 7918 8092 7475 7513 7448 7565 Gap (%) Xpress I. Greedy 1.74 2.13 2.16 3.17 1.00 2.31 0.35 2.39 0.14 1.75 3.23 4.77 3.83 4.78 1.16 3.29 1.46 1.96 2.45 3.96 Table 3.3: |N |=10, |T |=10, 3 different vessel types, number of vessels each method solves for vary for XpressMP and Improved Greedy as 10 and 1 in respect. All runs are terminated after 2 hours or when 0.1% gap from the best bound is reached The formulation proposed for Bender’s type approach is computationally efficient, as can be on Tables 3.4 and 3.5. The running time is a bit longer, but we’re able to get exact solutions. Xpress 7200 7200 7200 7200 7200 7200 7200 7200 7200 7200 Time (sec) Decomposition 767.1 757.5 905.8 1291.4 1420.4 959 1183.5 907.4 1032.8 1349.4 Objective Function Xpress Decomposition 7503 7380 8141 7973 8270 8196 7759 7739 7316 7313 8412 8148 8270 7961 7918 7834 7475 7373 7448 7273 Xpress 1.74 2.16 1.00 0.35 0.14 3.23 3.83 1.16 1.46 2.45 Gap (%) Decomposition 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 Table 3.4: |N |=10, |T |=10, 3 different vessel types, number of vessels each method solves for vary for XpressMP and Improved Greedy as 10 and 1 in respect. All runs are terminated after 2 hours or when 0.1% percent gap from the best bound is reached 43 Xpress 7200 7200 7200 7200 7200 7200 7200 7200 7200 7200 Time (sec) Decomposition 805 777.9 990.5 1401.1 1544.6 1015.4 1216.9 984.8 1114.2 1413.5 Objective Function Xpress Decomposition 7764 7619 8652 8448 8326 8245 8001 7978 7464 7461 8789 8496 9423 9026 8128 8035 7849 7727 7652 7455 Xpress 1.98 2.46 1.07 0.39 0.15 3.43 4.30 1.24 1.66 2.67 Gap (%) Decomposition 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 Table 3.5: |N |=15, |T |=10, 3 different vessel types, number of vessels each method solves for vary for XpressMP and Improved Greedy as 10 and 1 in respect. All runs are terminated after 2 hours or when 0.1% percent gap from the best bound is reached 3.4 Conclusion In this study, we studied a difficult real life supply chain scheduling problem encountered in oil and petrochemical industry, which involves production, inventory, and distribution operations, and requires an integrated scheduling to minimize the total operation cost. We showed the hardness of this problem, and showed that some of its special cases are polynomial time solvable. A minimum cost flow based heuristic, motivated by the observations from one of the special cases, was proposed and demonstrated to have a promising performance under the set of test cases considered in this study. Also, a new formulation of the model was developed, which made Bender’s type decomposition method computationally efficient. Therefore, we’re now able to get really good(exact) results for big problems at a much faster fashion then solver XpressMP. There are several interesting extensions of the work presented here. These include integrating the inland production with single or multiple refineries at different locations on the network, and multiple products needed by the same customer port. This integration would cause the supply chain to become bigger, and therefore more complex, however closer to reality, as inland production and demand satisfaction are activities that need 44 synchronization. Furthermore, the involvement of multiple refineries and multiple products introduces the new optimization issues due to assigning refineries to customer ports and allocating vessel capacity for different products. This will make the modeling and the design of search procedures more interesting and challenging. Also, for the simplicity of modeling, in this study, we assumed a linear penalty function for vessel overloading. However, this penalty cost is in reality very complex and is affected by many factors such as the level of overloading and navigation conditions. A nonlinear cost function would be more meaningful in this case. 45 Chapter 4 Container Vessel Scheduling with Time Windows 4.1 Introduction The problem we are interested in is related to the movement of full and empty containers in between ports. In this study we take the viewpoint of a single shipping company who delivers full containers from a single source port to a number of ports. We assume that these ports can naturally be visited in a given cyclic order, though not all ships have to stop at each of the ports in this cycle. We also assume that at each of the ports there is a known demand or surplus of empty containers, and the ships can use their spare capacity to load these empty containers and deliver them to the other 46 ports on their route. We consider a single planning horizon with certain length in which we would like to optimize the schedule and loading plans for the vessels. In this planning cycle, vessels are dispatched from the source port, carrying full and empty containers, within their capacity, to the other customer ports to satisfy their needs. Travel time between ports are deterministic, and the deliveries of the full containers must be made within given fixed time windows, defined uniquely for every port. If the vessels arrive prior to delivery time window, they must wait till the window opens, and incurring a waiting cost, where as if the delivery is made after the time window, a late fee is charged. Once a vessel is at a port and able to deliver, it remains there for the duration of berthing time associated with the particular port, and then the vessels can move to their next destination. Any unsatisfied handling of full containers or empty containers at a customer port incurs a penalty as well. There is no limitation on how many vessels could visit a single customer port. In addition to the capacity restrictions of vessels, vessels only have a fixed amount of time to visit ports and get back to the source port eventually. Also, the fleet size is not fixed, and we must decide how many vessels should be dispatched in the considered planning horizon. 47 4.2 Literature As there are a number of existing ship routing and related scheduling studies, most of them are covered by the following three major review papers: Ronen [42], Ronen [44] and Christiansen et al. [16]. A vast and recent literature review on container scheduling has been made by Choi et al. [14], covering different types of problems, their complexities, common approaches and algorithms. However, there is not too much emphasis on empty container problems. Gavish [21] developed a decision support system for vehicle fleet management. Its adaptation to the container fleet problem is straightforward. In his study, prior to making a decision on empty container relocation, owned and leased containers are assigned to the demand points based on the marginal cost criterion. Florez [20] presents a profit optimization model for the problem of empty container repositioning and leasing for ocean-going ships. He discusses the sensitivity of the model to the length of the planning horizon, and finds that the solution of the case study is affected little by the changes in the length of planning horizon. However, there are no full containers, nor time windows considered with this model as well. Crainic et al. [17] propose models for empty container allocation and distribution between a land transportation system and international maritime shipping network. Our problem does not consider any inland traffic, and with the current computational advancements, problems of much bigger scale can be dealt with. Cheung and Chen [11] compare a two-stage stochastic model with a two-stage deterministic model for the dynamic empty container allocation problem. In 48 their model, once a container is used to meet the demand at a port, it leaves the network and will re-enter the network as a random supply in the future.They perform some experiments with rolling horizons and conclude that a longer planning horizon is not necessarily better than a shorter one. Imai and Rivera IV [27] deal with fleet size planning for refrigerated containers where they determine the necessary number of containers required to meet predicted future transportation demand. More recently, Choong et al. [15] developed an integer programming formulation for empty container relocation with use of both long and short-term leased containers. The studies of Imai and Rivera IV [27] and Choong et al. [15] deal with empty container distribution in a relatively broad geographical area. In contrast, two recent papers focus on empty container repositioning in the hinterland of a specific port, in spite of the many similarities that exist in theory and in practice with repositioning in a board area. Li et al. [35] study the empty container allocation in a port with the aim to reduce redundant empty containers. They consider the problem as a non-standard inventory problem with simultaneous positive and negative demand under a general holding cost function. Jula et al. [28] consider empty container repositioning, which they refer to as empty container reuse, from a different perspective from that of the above studies. Their aim is the reduction of the traffic congestion in the Los Angeles and Long Beach port area caused heavily by empty maritime container traffic. A network flow formulation is constructed, in order to optimize empty container movements from consignees to shippers directly and/or via inland depots. The problem is solved in two phases: the first phase deals with the model transformation to a bipartite 49 transportation network (i.e. a classical Transportation Problem) and the second phase in solving the Transportation Problem by Linear Programming. Thus, it can be attested that no container ship routing studies and empty container management studies deal with an integrated and simultaneous approach to determine the optimal fleet composition with corresponding routing characteristics and empty container repositioning, where time windows are also considered. Li et al. [36] considers a problem where allocation of empty containers from supply ports to demand ports is taken care of. In this paper, optimal pairs of critical policies, (U, D) for one port, which are importing empty containers up to U when the number of empty containers in the port is less than U , or exporting empty containers down to D when the number of empty containers is larger than D, doing nothing otherwise, are adapted to multi-port case so that decision-makers can make decisions about allocating the right amounts of empty containers to the right ports at the right time. Highly related to the empty container reposition problem is the empty vehicle redistribution problem, in which much more literature exists. Du and Hall [18] used a single-value threshold control policy to redistribute empties in a hub-and-spoke network with random demands and deterministic travel times. Hall and Zhong [25] extended the above model to more general networks, where a pull-type decentralized control policy is adopted. Köchel et al. [31] presented Genetic Algorithms based simulation model to seek optimal fleet size and repositioning policy by maximizing the gain in the steady state. Their empty redistribution policies were based on the traditional inventory control policies (s, Q) and (s, S). However, the majority research on empty vehicle redistribution does 50 not consider the vehicle leasing decision and often assumed that unsatisfied demands are lost. However, in container shipping business, it is rare for a shipping line to reject a demand because of the unavailability of the owned containers. The main reason is that it is always possible for a shipping line to lease containers from vendors. Another important difference is that container reposition route is often fixed and should follow the vessel schedule, whereas empty vehicle redistribution is not. There is also quite bit of work on personnel scheduling, and Brucker et al. [6] have an extensive literature review with different models used, and their complexities. They also discuss NP-hard special cases. According to their categorization, ”A problem with Flexible Demand”, where tasks must be performed by employees within a time window is the type of problem that is most similar to our discussion. Some equivalence relations may be captured as in employee vs. vessel, or task vs. demand, but vessel fixed costs, holding or late fees, penalties for unsatisfied demand or empty containers cannot be incorporated. Robinson et al. [41] formulated the intraday scheduling problem as a maximum flow problem, which is similar to our formulation introduced below, however they require that the shifts are given. Later on Brucker and Qu [5] extended this work where all employees cannot perform all tasks. On the other hand, Ernst et al. [19] has an older staff scheduling and rostering literature review. It is more extensive than that of Brucker et al. [6], but is older. 51 4.3 Problem Definition Let us denote by V = {1, 2, ..., V } the set of available vessels, by P = {0, 1, ..., P } the set of ports, 0 being the source, and 1, 2, ..., P being their natural order to get visited. Let us denote by T = {1, 2, ..., T } the set of time period, where time is measured in some natural unit, say in days. In other words, our planning horizon is T days, and all time windows are assumed to be specified as subintervals of consecutive days within T. In the sequel we shall refer to vessels by v, w, ..., to ports by p, q, ... and to time units by t. We assume that each vessel v ∈ V has a capacity Uv , and a fixed cost Cv incurred if the vessel is dispatched. Travel times between ports p, q ∈ P, p < q are deterministic, denoted by Tp,q . Each customer port p ∈ P have a known demand for full containers, denoted by Fp , as well as demand for empty containers, denoted by Ep . We assume that Fp ≥ 0 for all ports p ∈ P, while Ep may take both negative and positive values, where Ep < 0 indicates the presence of a surplus of empty containers to be taken away. All deliveries made to a customer port p ∈ P expected to be within a fixed time window, [Ap , Bp ] ⊂ T. However vessels may arrive early or late with respect to the target time window. We assume that vessels berth only when they load/unload. To simplify our notations, we consider a particular day the ”arrival” day for a vessel at a port if by that day it unloads/loads its cargo. We assume that the travel time Tp,q includes the berthing time at port q. Accordingly, if a vessel arrives early to a port p ∈ P, a holding cost of CpH /day/container is charged for the unloaded cargo, and similarly if a delivery is too late, a penalty cost of CpL /day/container is charged. Customer ports may be visited by more than one vessel, but in case full container demand, or loading/unloading 52 of empty containers are not totally fulfilled, then costs CpF /container and CpE /container are charged, respectively. Finally, a revenue of Rp is received for satisfied full container demand at port p ∈ P. Let us note that since the total demand for full containers is given input for the considered problem, we get a mathematically equivalent formulation by simply assuming that Rp = 0 for all ports p ∈ P. To be able to formulate our model, we need to introduce the 53 following decision variables: V ariable Group T ype Explanation xv,p v ∈ V p ∈ P {0, 1} xv,p = 1 if vessel v visits port p yv v∈V av,p v∈Vp∈P {0, 1} yv = 1 if vessel v is dispatched number of full containers unloaded from vessel Z+ v at port p bv,p v∈Vp∈P number of empty containers unloaded from Z+ vessel v at port p cv,p v∈Vp∈P number of empty containers loaded by vessel Z+ v at port p fv,p v∈Vp∈P number of full containers on vessel v when Z+ leaving port p ev,p v∈Vp∈P number of empty containers on vessel v when Z+ leaving port p hv,p v ∈ V p ∈ P number of holding days of vessel v at port p Z+ due to early arrival ℓv,p v∈Vp∈P Z+ number of late days of vessel v at port p tv,p v∈Vp∈P Z+ arrival day (within T) of vessel v at port p zpF p∈P Z+ unsatisfied full conatiner demand at port p zpE p∈P Z+ unsatisfied empty container demand/pickup at port p Let us note that the variables describing vessel content, fv,p and ev,p are defined for all 54 vessels and ports, even if vessel v is not dispatched, or even if it does not visit port p. We shall make sure that these variables take value 0 for un-dispatched vessels, and we can view them, for dispatched vessel and unvisited ports, as the content of a passing by ship. We can now start describing our model in terms of these decision variables and parameters. Let us start with a few necessary relations, stemming form the logical relations these quantities must satisfy. Un-dispatched vessels cannot visit ports, and no loading or unloading at unvisited ports: xv,p ≤ yv ∀ v ∈ V, p ∈ P (4.1a) av,p ≤ xv,p ∀ v ∈ V, p ∈ P (4.1b) bv,p ≤ xv,p ∀ v ∈ V, p ∈ P (4.1c) cv,p ≤ xv,p ∀ v ∈ V, p ∈ P (4.1d) Number of containers on board must obey the law of conservation, and cannot exceed the ship’s capacity: fv,p = fv,p−1 − av,p ∀ v ∈ V, p ∈ P \ {0} (4.1e) ev,p = ev,p−1 − bv,p + cv,p ∀ v ∈ V, p ∈ P \ {0} (4.1f) ev,p + fv,p ≤ Uv ∀ v ∈ V, p ∈ P (4.1g) Let us note that we defined all decision variables to be nonnegative. Thus, the above relations imply also that we cannot unload more than what we have on a ship, and we 55 cannot load more than the ship’s free capacity. The next group of constraints make sure that we load/unload empty containers only where there is a demand/surplus: bv,p = 0 ∀ p ∈ P with Ep ≤ 0, and ∀ v ∈ V (4.1h) cv,p = 0 ∀ p ∈ P with Ep ≥ 0, and ∀ v ∈ V (4.1i) Finally, we prescribe that container demands are met: Fp = zpF + X av,p ∀p∈P (4.1j) ∀p∈P (4.1k) v∈V Ep =   X   E  bv,p + z  p  if Ep ≥ 0 v∈V  X   E  − cv,p −z  p  if Ep ≤ 0 v∈V Note that we defined Ep to take both positive (demand) and negative (surplus) values, and hence we had to formulate the above constraints accordingly. The next group of constraints will help to get the right values assigned to the time related decision variables. If the same vessel visits two ports, then the arrival times must differ by at least the travel time between these ports: tv,p ≥ tv,q + Tq,p + T (xv,p + xv,q − 2) ∀ v ∈ V, p, q ∈ P, q < p (4.2a) 56 Note that since T is a maximum time difference in our problem, the last term makes the above inequality irrelevant, unless vessel v visits both ports p and q. Ap ≤ hv,p + tv,p + T (1 − xv,p ) ∀ v ∈ V, p ∈ P (4.2b) Bp ≥ tv,p − ℓv,p − T (1 − xv,p ) ∀ v ∈ V, p ∈ P (4.2c) Note that the last term makes the above constraints trivial, unless vessel v visits port p. Note also that both hv,p and ℓv,p are limited form below by these constraints. This, together with the fact that they will have nonnegative coefficients in our minimization objective, assures that our objective computes penalties for the true earliness/lateness. The last group of constraints is to help to formulate our objective, which is a compositions of all expenses (penalties, etc.). Vessel dispatching costs: cD = X C v yv (4.3a) v∈V Early and late arrival penalties: cH = X CpH X CpL p∈P cL = p∈P X hv,p (av,p + bv,p ) (4.3b) X ℓv,p (av,p + cv,p ) (4.3c) v∈V v∈V 57 Unsatisfied demands’ penalties: cU = X CpF zpF + CpE zpE (4.3d) p∈P Then our problem is to minimize cD + cH + cL + cU , subject to the constraints 4.1-4.3 and the sign and type restrictions in the definitions of the decision variables. Note that constraints 4.3b and 4.3c are quadratic in terms of our decision variables. All other constraints are linear. Note also, that if variables x and y are binary, then the integrality of the other decision variables follow from the structure of the constraints and the integrality of the input parameters. 4.4 Minimum Cost Flow Network Formulation We propose the following setup to make the handling of the problem easier. Initially, we only take into consideration full containers. We create a source port node and a penalty node for every day in our planning horizon, starting from day 1, resulting in overall 2T nodes. The next set of nodes are associated with every port and every day of the planning horizon, so a total of P ∗ T are required. The arcs emanating from the source nodes are connected to the port nodes, based on T0,p , i.e. if T0,1 = 3, then there would be arcs between every source node, (0, t), and nodes associated with port 1, (1, t + 3). These arcs have 0 lower bound, no upper bound, and Cv /Uv cost, where v is such that Cv /Uv is maximum ∀v ∈ V. The outgoing arcs from penalty nodes are connected to 58 the nodes associated with day Bp for p ∈ P, have 0 lower bound, no upper bound, and a unit cost of CpF . The (p, t) nodes, if t < Ap , have outgoing arcs to (p, t + 1), with 0 lower bound, no upper bound and a cost of CpH . The (p, t) nodes, if Ap ≤ t < Bp , have 2 different categories of outgoing arcs. One group consists of the flow within the same port, from node (p, t) to (p, t + 1), with 0 lower bound, no upper bound, and no cost. The other group of arcs represent the flows to other ports and again, arcs exist between two nodes associated with ports p and q based on Tp,q . (p, t), where t = Bp have one additional outgoing arc with respect to the set of nodes that we just discussed, and these nodes are connected to the sink node, with a lower bound and upper bound of Fp , and 0 cost. This ensures that either via the flow through the source port or through penalty nodes, demand constrained will be satisfied. Finally, for set of (p, t) nodes, where t > Bp , there is only one outgoing arc to (p, t − 1), with 0 lower bound, no upper bound, and a unit cost of CpL . Finally, we need to somehow consider vessel times, and in order to do that we make sure that Bp for p ∈ P, is such that Bp + Tp,0 ≤ T . Figure 4.1 shows an example of the network for |T | = 7, |P | = 2, Tp,q = 1 ∀p, q ∈ P, and Ap , Bp ∀p ∈ P are indicated by brackets around (p, t) nodes. 4.5 Minimum Cost Flow Heuristic With the above formulation, we propose to solve the above problem in 3 stages. In the first stage, we only worry about taking care of full container demand. The second stage is handled with the same minimum cost flow formulation, where flows on arcs are converted to vessels. Finally, on the third stage, based on the paths vessels have been assigned to at the end of second stage, another minimum cost flow problem is solved, 59 Figure 4.1: An example of the network for |T | = 7, |P | = 2, Tp,q = 1 ∀p, q ∈ P, and time windows for delivery are indicated by brackets around the nodes where we try to satisfy the empty container demand to our best with the remaining capacities on the vessels. 4.5.1 Stage 1 The solution to this minimum cost flow problem yields flow of full containers, that are carried from source to ports, and in between ports. Once these flow values are obtained, we move to the second stage. 4.5.2 Stage 2 Keeping the same network structure, all arc costs are removed. We divide the flow on the arcs by Uv and round it up, and set this quantity as the lower and upper bound on the capacity of the arcs. Once we solve this problem, we are able to tell exactly how 60 many vessels were dispatched, moreover, greedily, we can determine unique paths for each of the vessels, based on the flow on the arcs. Once the path of each vessel is known, we move to the third stage. 4.5.3 Stage 3 A new network is formed, where arcs between ports only exist based on the information coming from the second stage. The arc capacities are equal to the remaining capacities of the vessels based on their full container load at the time of departure from a port. If the ports have surplus of empty containers, then we have additional nodes that carry this surplus to the nodes, and if the ports are in demand of empty containers, then the outgoing arcs from such ports have a lower bound that is equal to the port’s demand. Once this problem is solved, the empty containers are taken care off as well, and the heuristic terminates. The results obtained with this approach will be displayed after the next section, as well as the discussion on its performance. 4.6 Decomposition Approach The proposed heuristic maybe misleading in terms of the number of vessels used, as during Stage 1, vessel fixed cost is incorporated as cost per unit carried, and also during Stage 2, the rounding up for the number of vessels used contributes to this. Therefore, 61 a better way of fleet management is necessary, and we propose the following formulation in correspondence with the network structure described above for the full container problem. The proposed approach will replace the first two stages of our heuristic. The main difference of the formulation below is that the source port nodes are replaced by vessel nodes, based on the dispatch info. To be able to formulate our model, we need to introduce the following decision variables: V ariable Group T ype Explanation number of full containers car- xv,t,p,t′ v ∈ V t ∈ T p ∈ P t′ ∈ T ried by vessel v, which was disZ+ patched at period t, to unload at port p arriving at t′ yv,t v∈Vt∈T {0, 1} yv,t = 1 if vessel v is dispatched at period t number of full containers car- xp,t,q,t′ p ∈ P t ∈ T q ∈ P t′ ∈ T ried shipped from port p at Z+ time t arriving to port q at time t′ number of full container de- zp,t p∈Pt∈T Z+ mand not satisfied at port p due at time t We can now start describing our model in terms of these decision variables and parameters. 62 Vessel dispatching costs: cD = XX Cv yv,t (4.4a) X X CpH xp,t,p,t+1 (4.4b) X X CpL xp,t,p,t−1 (4.4c) v∈V t∈T Early and late arrival penalties: cH = p∈P t<Ap cL = p∈P t>Bp Unsatisfied demands’ penalties: cU = XX CpF zp,t (4.4d) p∈P t∈T Constraint related to dispatching of vessels will be: xv,t,p,t′ ≤ yv,t Uv ∀ v ∈ V, t ∈ T, p ∈ P, t′ > t (4.5) 63 Next set of constraints are flow-balance constraints. X xv,t′ ,p,t + X xq,t′′ ,p,t + xp,t−1,p,t = xp,t,p,t+1 (4.6a) q<p v∈V ∀ p ∈ P, t < Ap , t′ = t − T0,p , t′′ = t − Tq,p X xv,t′ ,p,t + X xq,t′′ ,p,t + xp,t−1,p,t = xp,t,p,t+1 + q<p v∈V X xp,t,r,t′′′ (4.6b) r>p ∀ p ∈ P, Ap ≤ t < Bp , t′ = t − T0,p , t′′ = t − Tq,p , t′′′ = t + Tp,r X xv,t′ ,p,t + X q<p v∈V xq,t′′ ,p,t + xp,t−1,p,t + zp,t = X xp,t,r,t′′′ + xp,t,0,t′′′′ (4.6c) r>p ∀ p ∈ P, t = Bp , t′ = t − T0,p , t′′ = t − Tq,p , t′′′ = t + Tp,r , t′′′′ = t + Tp,0 X v∈V xv,t′ ,p,t + X xq,t′′ ,p,t + xp,t+1,p,t = xp,t,p,t−1 (4.6d) q<p ∀ p ∈ P, t > Bp , t′ = t − T0,p , t′′ = t − Tq,p Finally, to ensure that the demand is satisfied via vessels used or penalty nodes: xp,t,0,t′ = Fp ∀ p ∈ P, t = Bp , t′ = t + Tp,0 (4.7) This formulation makes decomposition rather easy to implement as the binary variable, yv,t appears only in one constraint, 4.5, and in one part of the objective function, 4.4a. Also, realize that, all equations except 4.5 and 4.7 have 0 right hand side. 64 Xpress 7200 7200 7200 7200 7200 7200 7200 7200 7200 7200 Time (sec) Heuristic 161.56 127.06 150.61 164.87 160.14 107.42 135.95 95.68 149.17 109.10 Decomp 1634.85 1294.33 1025.03 1191.42 1708.97 984.12 1324.92 1680.62 1256.40 933.26 Objective Function Xpress Heuristic Decomp 16259 16397 15954 9523 9578 9298 8823 8827 8737 15112 15425 15067 16431 16506 16423 14187 14389 13714 14672 14724 14054 9364 9472 9257 11307 11543 11130 10285 10387 10020 Xpress 1.97 2.46 1.07 0.39 0.15 3.43 4.30 1.24 1.66 2.67 Gap (%) Heuristic 2.80 3.02 1.11 2.41 0.60 4.78 4.64 2.37 3.67 3.63 Decomp 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 Table 4.1: |N | = 20, |V | = 40, |T | = 30, 0.1% gap termination condition or 2 hrs for both XpressMP and Decomposition Xpress 7200 7200 7200 7200 7200 7200 7200 7200 7200 7200 Time (sec) Heuristic 168.33 316.45 315.44 202.37 242.31 193.69 302.12 169.61 171.10 203.57 Decomp 1904.57 1788.67 1450.06 1782.84 2126.91 1152.38 1649.85 2161.24 1912.81 1299.79 Objective Function Xpress Heuristic Decomp 23817 24494 23426 16179 16570 15846 20663 21120 20477 15909 16400 15869 17306 17809 17299 21797 22199 21114 20302 20594 19544 15361 15779 15197 17840 17899 17598 14137 14544 13805 Xpress 1.74 2.16 1.00 0.35 0.14 3.23 3.83 1.16 1.46 2.45 Gap (%) Heuristic 4.45 4.47 3.14 3.34 2.96 4.98 5.19 3.78 1.78 5.18 Decomp 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10 Table 4.2: |N | = 40, |V | = 80, |T | = 30, 0.1% gap termination condition or 2 hours for both XpressMP and Decomposition 4.7 Results Both 4.1 and 4.2 suggest that decomposition approach performs in general better than XpressMP for the full container problem. The accuracy of heuristic on the other hand varies, however it takes significantly shorter for it to run in both cases. The weakness of the heuristic is mainly due to the setup on Stage 1. Unfortunately, vessel fixed costs are incorporated as unit carriage costs to simplify formulation, and therefore the heuristic is not able to realize that an additional vessel may be used, where only a small portion of its capacity is utilized. As far as the number of vessels used are concerned, the network heuristic usually requires 1 or 2 more vessels than the number of vessels used with XpressMP. 65 Chapter 5 Conclusion This dissertation addressed different problems, first one being related to the issue of managing inventories of medical supplies and especially the critical ones in hospitals under surge (pandemic) scenarios. Inventory management basically deals with how much of each unit to maintain on hand by deciding on when to order and how many to order for each item under consideration. This is a critical issue since too much of medical supply inventory can cost a lot and perhaps perish, while short inventories will not satisfy the demand. Finding the right amount of stock is always a challenge in inventory management due to uncertainties involved in demand. We have developed a sound and practical approach that combines epidemiologic modeling techniques with simulation and optimization modeling to provide the best strategy for managing inventories. It involves a high-fidelity Disease Progress Module Influenza Pandemic like scenarios using already validated data from the historical epidemiological literature. For the proposed simulation framework, a Virtual Hospital Module was developed to capture resource 66 consumption in healthcare settings during an Influenza Pandemic. A dynamic programming optimization model was constructed to optimally manage the inventory of critical medical supplies (that is the decisions regarding when to order and how much to order) in hospital settings. A number of numerical scenarios were analyzed results are obtained. The approach is quite practical and readily implementable in hospital settings. In addition, demand forecasting tools using a hidden Markov model and a deterministic model have been developed, enabling us to deal with unknown pandemic as well. The strengths and weaknesses of these forecasting methods are displayed, and discussion on how these weaknesses can be overcome is made. As mentioned earlier, this work aims to provide guidelines to implement a formal procedure to help decision makers managing inventories in the health care industry. Due to that fact that shortcomings in medical supplies may end up in dire consequences, inventories are typically held at higher levels as compared to other industries. Clearly, this has consequences that appear as significant contributions to health care costs. A serious effort to reduce health care related costs is to implement effective inventory management policies so that supply stocks are kept at reasonable levels during times when there is no mass demand. During emergency scenarios such as pandemic situations the system should be effective enough to respond quickly to build right amount of inventories. This research can be extended to the cases where the demand is following another distribution, or where products perish, or by adding another layer of the supply chain to 67 the problem, as well as investing more in the demand forecasting. Next study was involved with the study of a difficult real life supply chain scheduling problem encountered in oil and petrochemical industry, which involves production, inventory, and distribution operations, and requires an integrated scheduling to minimize the total operation cost. We showed the hardness of this problem, and showed that some of its special cases are polynomial time solvable. A minimum cost flow based heuristic, motivated by the observations from one of the special cases, was proposed and demonstrated to have a promising performance under the set of test cases considered in this study. Also, a new formulation of the model was developed, which made Bender’s type decomposition method computationally efficient. Therefore, we’re now able to get really good results for big problems at a much faster fashion then solver XpressMP. There are several interesting extensions of the work presented here. These include integrating the inland production with single or multiple refineries at different locations on the network, and multiple products needed by the same customer port. For the simplicity of modeling in this study, we assumed a linear penalty function for vessel overloading. However, this penalty cost is in reality very complex and is affected by many factors such as the level of overloading and navigation conditions. A nonlinear cost function would be more meaningful in this case. Furthermore, the involvement of multiple refineries and multiple products introduces the new optimization issues due to assigning refineries to customer ports and allocating vessel capacity for different products. This will make 68 the modeling and the design of search procedures more interesting and challenging. Last of all, we have dealt with a very general model for container vessel scheduling with time windows constraints. We have managed to devise a fast heuristic, as well as, an efficient formulation, making Bender’s type approach suitable. Bender’s type approach performs much better than solver XpressMP, both in terms of time, and error bound. Overall, for problems of smaller size, the heuristic we proposed yields fast solutions with an acceptable gap from the optimal, where as, as the problem gets bigger, decomposition approach becomes the better option. The contribution of this work is the formulation making Bender’s type approach applicable in a useful way. Further study can be developed, where inland shipment of containers, thus becoming empty is included as well. However, there is a shortcoming of the proposed Bender’s type approach. The empty containers are not dealt with until the full containers are taken care of. Therefore, we are not able to guarantee optimality within the big picture, which is not surprising in the end, as with full and empty containers, the problem becomes a variant of multicommodity flow problems, which are hard in general. Empty and full containers, because of the way they are handled, on the other hand, may allow an easier formulation. 69 Bibliography [1] R.M. Anderson, R.M. May, and B. Anderson. Infectious diseases of humans: dynamics and control, volume 28. Wiley Online Library, 1992. [2] H. Andersson, A. Hoff, M. Christiansen, G. Hasle, and A. Løkketangen. Industrial aspects and literature survey: Combined inventory management and routing. Computers & Operations Research, 37(9):1515–1536, 2010. [3] HB Bendall and AF Stent. A scheduling model for a high speed containership service: A hub and spoke short-sea application. International Journal of Maritime Economics, 3(3):262–277, 2001. [4] J.H. Bookbinder and K.E. Reece. Vehicle routing considerations in distribution system design. European Journal of Operational Research, 37(2):204–213, 1988. [5] P. Brucker and R. Qu. Network flow models for intraday personnel scheduling problems. [6] P. Brucker, R. Qu, and E. Burke. Personnel scheduling: Models and complexity. European Journal of Operational Research, 210(3):467–473, 2011. [7] JA Buzacott and JG Shanthikumar. Safety stock versus safety time in mrp controlled production systems. Management science, 40(12):1678–1689, 1994. [8] M. Chan. World now at the start of 2009 influenza pandemic. World Health Organization, 11, 2009. [9] C.A. Chang. The interchangeability of safety stocks and safety lead time. Journal of Operations Management, 6(1):35–42, 1985. [10] L. Cheng and M.A. Duran. Logistics for world-wide crude oil transportation using discrete event simulation and optimal control. Computers & chemical engineering, 28(6):897–911, 2004. [11] R.K. Cheung and C.Y. Chen. A two-stage stochastic network model and solution methods for the dynamic empty container allocation problem. Transportation Science, 32(2):142–162, 1998. [12] S.C. Cho and AN Perakis. Optimal liner fleet routeing strategies. Maritime Policy and Management, 23(3):249–259, 1996. 70 [13] S.C. Cho and A.N. Perakis. An improved formulation for bulk cargo ship scheduling with a single loading port. Maritime Policy & Management, 28(4):339–345, 2001. [14] B.C. Choi, K. Lee, J.Y.T. Leung, M.L. Pinedo, and D. Briskorn. Container scheduling: Complexity and algorithms. Production and Operations Management, 21(1): 115–128, 2012. [15] S.T. Choong, M.H. Cole, and E. Kutanoglu. Empty container management for intermodal transportation networks. Transportation Research Part E: Logistics and Transportation Review, 38(6):423–438, 2002. [16] M. Christiansen, K. Fagerholt, and D. Ronen. Ship routing and scheduling: Status and perspectives. Transportation Science, 38(1):1–18, 2004. [17] T.G. Crainic, M. Gendreau, and P. Dejax. Dynamic and stochastic models for the allocation of empty containers. Operations Research, 41(1):102–126, 1993. [18] Y. Du and R. Hall. Fleet sizing and empty equipment redistribution for centerterminal transportation networks. Management Science, 43(2):145–157, 1997. [19] A.T. Ernst, H. Jiang, M. Krishnamoorthy, and D. Sier. Staff scheduling and rostering: A review of applications, methods and models. European journal of operational research, 153(1):3–27, 2004. [20] H. Florez. Empty-containers repositioning and leasing: An optimization model. Dissertation Abstracts International Part B: Science and Engineering[DISS. ABST. INT. PT. B- SCI. & ENG.],, 46(12), 1986. [21] B. Gavish. A decision support system for managing the transportation needs of a large corporation. AIIE Transactions, 13(1):61–85, 1981. [22] S. Geary, S.M. Disney, and D.R. Towill. On bullwhip in supply chainshistorical review, present practice and expected future impact. International Journal of Production Economics, 101(1):2–18, 2006. [23] F. Gheysens, B. Golden, and A. Assad. A comparison of techniques for solving the fleet size and mix vehicle routing problem. OR Spectrum, 6(4):207–216, 1984. [24] S.C. Graves. A single-item inventory model for a nonstationary demand process. Manufacturing & Service Operations Management, 1(1):50–61, 1999. [25] R.W. Hall and H. Zhong. Decentralized inventory control policies for equipment management in a many-to-many network. Transportation Research Part A: Policy and Practice, 36(10):849–865, 2002. [26] A. Hoff, H. Andersson, M. Christiansen, G. Hasle, and A. Løkketangen. Industrial aspects and literature survey: Fleet composition and routing. Computers & Operations Research, 37(12):2041–2061, 2010. [27] A. Imai and F. Rivera IV. Strategic fleet size planning for maritime refrigerated containers. Maritime Policy & Management, 28(4):361–374, 2001. 71 [28] H. Jula, A. Chassiakos, and P. Ioannou. Port dynamic empty container reuse. Transportation Research Part E: Logistics and Transportation Review, 42(1):43–60, 2006. [29] S. Karlin. Dynamic inventory policy with varying stochastic demands. Management Science, 6(3):231–258, 1960. [30] S. Karlin and H. Scarf. Inventory models of the arrow-harris-marschak type with time lag. Studies in the Mathematical Theory of Inventory and Production, pages 155–178, 1958. [31] P. Köchel, S. Kunze, and U. Nieländer. Optimal control of a distributed service system with moving resources: Application to the fleet sizing and allocation problem. International Journal of Production Economics, 81:443–459, 2003. [32] Y. Le Strat and F. Carrat. Monitoring epidemiologic surveillance data using hidden markov models. Statistics in medicine, 18(24):3463–3478, 1999. [33] L. Lei, H. Zhong, and W.A. Chaovalitwongse. On the integrated production and distribution problem with bidirectional flows. INFORMS Journal on Computing, 21(4):585–598, 2009. [34] R. Levi, G. Janakiraman, and M. Nagarajan. A 2-approximation algorithm for stochastic inventory control models with lost sales. Mathematics of Operations Research, 33(2):351–374, 2008. [35] J.A. Li, K. Liu, S.C.H. Leung, and K.K. Lai. Empty container management in a port with long-run average criterion. Mathematical and Computer Modelling, 40 (1):85–100, 2004. [36] J.A. Li, S.C.H. Leung, Y. Wu, and K. Liu. Allocation of empty containers between multi-ports. European Journal of Operational Research, 182(1):400–412, 2007. [37] L.G. Mallon and J.P. Magaddino. An integrated approach to managing local container traffic growth in the long beach–los angeles port complex, phase ii. Technical report, 2001. [38] T.E. Morton. Bounds on the solution of the lagged optimal inventory equation with no demand backlogging and proportional costs. SIAM review, 11(4):572–596, 1969. [39] TAJ Nicholson and RD Pullen. Dynamic programming applied to ship fleet management. Operational Research Quarterly, pages 211–220, 1971. [40] A. Riley. The coming of the russian gas deficit: Consequences and solutions. CEPS Policy Briefs, (1-12):1, 2006. [41] R. Robinson, R. Sorli, and Y. Zinder. Personnel scheduling with time windows and preemptive tasks. In Proceedings of the 5th International Conference on the Practice and Theory of Automated Timetabling, 18th August-20th August, pages 561–566, 2004. 72 [42] D. Ronen. Cargo ships routing and scheduling: Survey of models and problems. European Journal of Operational Research, 12(2):119–126, 1983. [43] D. Ronen. Short-term scheduling of vessels for shipping bulk or semi-bulk commodities originating in a single area. Operations Research, 34(1):164–173, 1986. [44] D. Ronen. Ship scheduling: The last decade. European Journal of Operational Research, 71(3):325–333, 1993. [45] R. Salomon and R.G. Webster. The influenza virus enigma. Cell, 136(3):402–410, 2009. [46] T. Schoenmeyr and S.C. Graves. Strategic safety stocks in supply chains with evolving forecasts. Manufacturing & Service Operations Management, 11(4):657– 673, 2009. [47] J.S. Song and P. Zipkin. Inventory control in a fluctuating demand environment. Operations Research, 41(2):351–370, 1993. [48] A.F. Veinott Jr. Optimal policy in a dynamic, single product, nonstationary inventory model with several demand classes. Operations Research, 13(5):761–778, 1965. [49] P. Wu, J.C. Hartman, and G.R. Wilson. An integrated model and solution approach for fleet sizing with heterogeneous assets. Transportation Science, 39(1):87–103, 2005. [50] X. Xinlian, W. Tengfei, and C. Daisong. A dynamic model and algorithm for fleet planning. Maritime Policy & Management, 27(1):53–63, 2000. [51] Y.S. Zheng and A. Federgruen. Finding optimal (s, s) policies is about as simple as evaluating a single policy. Operations Research, 39(4):654–665, 1991.