Modeling Service Performance and Dynamic Worker Allocation Policies for Order Ful llment Systems by Hyun Ho Kim A dissertation submitted to the Graduate Faculty of Auburn University in partial ful llment of the requirements for the Degree of Doctor of Philosophy Auburn, Alabama December 18, 2009 Keywords: Multi-server queues, approximation models, phase-type distributions, steady state sojourn time distributions, state-dependent sojourn time distributions, order ful llment systems, probability of success, dynamic worker allocation Copyright 2009 by Hyun Ho Kim Approved by: Kevin R. Gue, Chair, Associate Professor of Industrial & Systems Engineering Ming Liao, Professor of Mathematics & Statistics Emmett J. Lodree, Jr., Assistant Professor of Statistics & Management Science, University of Alabama Abstract In this dissertation, we propose several dynamic worker allocation policies based on the sojourn time distribution in due-date order ful llment systems. To establish a sojourn- based policy, we must be able to compute the sojourn time distribution for an arriving order and each order in the system. We require the distribution because we must make a probabilistic statement about the completion time of an order. We need two types of sojourn time distribution for the system. The rst is a steady-state distribution to estimate the probability of success for an arriving order. The second is the state-dependent sojourn time distribution, which estimates the probability of success for an order in system. First, we develop an approximation model for the sojourn time distribution of customers or orders arriving to an acyclic multi-server queueing network. The model accepts general interarrival times and general service times, and is based on the characteristics of phase- type distributions. Distributions produced by the model agree well with those produced by simulation for a variety of serial and general queueing networks. Second, we develop an approximation model for the state-dependent sojourn time dis- tribution of customers or orders in a multi-server queueing system, when interarrival and service times can take on general distributions. The model is based on the characteristics of phase-type distributions, and uses a Markov process to represent the all-busy period for a waiting time distribution. We also discuss how the model might be used to execute dynamic delivery promises in an order ful llment system. Third, we propose several dynamic worker allocation policies to maximize service perfor- mance of an order ful llment system. Our policies make use of the state-dependent sojourn time distribution for an order, which we compute with a model based on phase-type dis- tributions. Our approach di ers from other work on dynamic worker allocation in that we ii focus on service performance as perceived by the customer, instead of traditional system performance measures such as cycle time, throughput, and WIP. Our results suggest that order ful llment systems can signi cantly improve their service performance by moving the right number of workers to the right place, at the right time. iii Acknowledgments I would like to acknowledge the Republic of Korea Army for giving me the opportunity to study and for the support during the course of my research. When I look back on my past, I have received many bene ts from the Army. I have had chances to earn an education from Korea Military Academy to a Ph.D. and it has created the essence of my character as a man and soldier. From now on, I think it is time to pay back my debts more than I received from ROKA. My rst thanks must go to my advisor, Dr. Kevin R. Gue. He was always willing to provide advice on academic work and help me along my way. I have to attribute most of my success to Dr. Gue?s guidance. Through most of my years here, I also have learned about life from him. He thoughtfully considers others, tries as best as he can to live a moderate life, and has a serious attitude on life. I am so lucky to have him as my mentor. I am deeply grateful to Dr. Ming Liao and Dr. Emmett J. Lodree, Jr for their help and support, especially I was able to gure out a lot of problems for my dissertation with Dr. Liao?s classes and his help. I shall never forget attending a benevolent professor?s classes. I also give thanks to an excellent batch of fellow graduate students, who I will not mention by name for fear of leaving someone out. They made the social, academic and all other aspects of my years at Auburn University the most enjoyable of my life. Here I am in mid-life without quite realizing it. I am absolutely delighted that I did something to delight my parents as a son. I could not nish my work without their devoted love. This dissertation is dedicated to my wife Hyesook. I owe much of my success to my wife and I thank my wife, Hyesook, for all of her love and support during the last four years. In an unfamiliar foreign land, my family (Hyesook and Siheon) was a great comfort and help for me. Life in Auburn will stay with our family as a happy memory. iv I aspire to guide my son as well as my parents guided me! v Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv List of Illustrations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xii 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 The order ful llment system . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Problem Statement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.3 Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.4 Organization of the dissertation . . . . . . . . . . . . . . . . . . . . . . . . . 7 2 An Approximation Model for Sojourn Time Distributions in Acyclic Multi-server Queueing Networks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.2 Related literature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.3 Phase-type distributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.4 Intuitive model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.4.1 The ow rate and SCV of the arrival process . . . . . . . . . . . . . . 16 2.4.2 Waiting time for a single stage queue with multiple servers . . . . . . 17 2.4.3 Approximation of sojourn time distribution of the queueing network . 19 2.4.4 Numerical results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.5 Interpolation model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 2.5.1 Validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.6 Determining an optimal cuto time . . . . . . . . . . . . . . . . . . . . . . . 34 2.7 Chapter conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 vi 3 Predicting Departure Times in Multi-Stage Queueing Systems . . . . . . . . . . 39 3.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 3.2 Preliminaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 3.2.1 Phase-type distributions . . . . . . . . . . . . . . . . . . . . . . . . . 43 3.2.2 Fitting distribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 3.3 Exponential service times . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 3.3.1 Approximation model . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 3.3.2 Numerical results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 3.4 General service times . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 3.4.1 Approximating the rst epoch . . . . . . . . . . . . . . . . . . . . . . 50 3.4.2 Approximating the remaining epochs . . . . . . . . . . . . . . . . . . 53 3.4.3 Approximating the sojourn time distribution of the system . . . . . . 54 3.4.4 Normal Approximation Model . . . . . . . . . . . . . . . . . . . . . . 55 3.4.5 Numerical results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 3.4.6 Testing the Normal Approximation Model . . . . . . . . . . . . . . . 56 3.5 Multi-Stage Systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 3.5.1 Approximation model . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 3.5.2 Numerical results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 3.6 Chapter conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 4 Dynamic Worker Allocation Policies to Improve Service in an Order Ful llment System . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 4.2 Literature Review . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 4.3 Preliminaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 4.3.1 Next Scheduled Departure (NSD) . . . . . . . . . . . . . . . . . . . . 75 4.3.2 State-dependent sojourn time distribution . . . . . . . . . . . . . . . 77 4.3.3 The probability of success . . . . . . . . . . . . . . . . . . . . . . . . 78 vii 4.4 Worker allocation policies . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 4.5 Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 4.5.1 Single Flush . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 4.5.2 Multi-Flush . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 4.5.3 Cascade policy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 4.6 Chapter conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 5 Conclusions and Outlook . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 Appendices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 viii List of Illustrations 1.1 A nominal order ful llment system . . . . . . . . . . . . . . . . . . . . . . . 1 2.1 Sojourn time distributions for Example 1, when all 1=C2 are integers. . . . . 21 2.2 A comparison of sojourn time distributions between the intuitive model and simulation when all 1=C2 are not integers. . . . . . . . . . . . . . . . . . . . 22 2.3 A comparison of mean sojourn times among three models( =0.5). Dots indi- cate data; lines are inserted for visual clarity only. . . . . . . . . . . . . . . . 25 2.4 A serial line with 3 workstations. . . . . . . . . . . . . . . . . . . . . . . . . 25 2.5 A general network with 4 workstations. . . . . . . . . . . . . . . . . . . . . . 26 2.6 CDF and PDF of the serial line ( =0.85, C2 = 0:75; 0:75; 0:75; 0:75). . . . 27 2.7 CDF and PDF of the serial line ( =0.85, C2 = 0:57; 0:7; 0:6; 0:9). . . . . . 28 2.8 CDF and PDF of the serial line ( =0.85, C2 = 0:4; 0:6; 0:26; 0:8). . . . . . 28 2.9 CDF and PDF of the serial line ( =0.5, C2 = 0:75; 0:75; 0:75; 0:75). . . . . 28 2.10 CDF and PDF of the serial line ( =0.5, C2 = 0:57; 0:7; 0:6; 0:9). . . . . . . 28 2.11 CDF and PDF of the serial line ( =0.5, C2 = 0:4; 0:6; 0:26; 0:8). . . . . . . 29 2.12 CDF and PDF of the general network( =0.85,C2 = 0:75;0:75;0:75;0:75;0:75). 30 2.13 CDF and PDF of the general network ( =0.85, C2 = 0:57; 0:7; 0:6; 0:6; 0:9). 30 2.14 CDF and PDF of the general network ( =0.85, C2 = 0:4; 0:6; 0:26; 0:26; 0:8). 31 2.15 CDF and PDF of the general network ( =0.5,C2 = 0:75;0:75;0:75;0:75;0:75). 31 2.16 CDF and PDF of the general network ( =0.5, C2 = 0:57; 0:7; 0:6; 0:6; 0:9). 31 2.17 CDF and PDF of the general network ( =0.5, C2 = 0:4; 0:6; 0:26; 0:26; 0:8). 31 2.18 PDF and CDF of the Example 2.5. . . . . . . . . . . . . . . . . . . . . . . . 33 ix 2.19 PDF and CDF for the long serial line in Example 2.6. . . . . . . . . . . . . . 34 2.20 Relationship between cuto time and probability of success of an arriving order. 35 2.21 Pro t curve for a range of cuto times. . . . . . . . . . . . . . . . . . . . . . 37 3.1 Comparison of the PDF and CDF of the approximation model and the simu- lation model (2 servers and 20 orders ahead). . . . . . . . . . . . . . . . . . 48 3.2 Comparison of the PDF and CDF of the approximation model and the simu- lation model (50 servers and 10 orders ahead). . . . . . . . . . . . . . . . . . 48 3.3 Comparison of the PDF and CDF of the approximation model and the simu- lation model (100 servers and 5 orders ahead). . . . . . . . . . . . . . . . . 48 3.4 System-state and server-state. . . . . . . . . . . . . . . . . . . . . . . . . . . 51 3.5 Comparison of the PDF and CDF of the approximation model and the simu- lation model (2 servers and 20 orders ahead). . . . . . . . . . . . . . . . . . . 57 3.6 Comparison of the PDF and CDF of the approximation model and the simu- lation model (50 servers and 10 orders ahead). . . . . . . . . . . . . . . . . . 58 3.7 Comparison of the PDF and CDF of the approximation model and the simu- lation model (100 servers and 5 orders ahead). . . . . . . . . . . . . . . . . 58 3.8 State-dependent queueing network model. Black discs represent workers; squares represent orders. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 3.9 A serial line with 3 stations. . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 3.10 An acyclic queueing network with 4 stations. Black discs represent workers. . 62 3.11 Comparison of the PDF and CDF of the serial line (11-12-13). . . . . . . . . 64 3.12 Comparison of the PDF and CDF of the serial line (6-12-13). . . . . . . . . . 64 3.13 Comparison of the PDF and CDF of the acyclic queueing network (The rst station, 11-12-8-13). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 3.14 Comparison of the PDF and CDF of the acyclic queueing network (The rst station, 6-4-16-13). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 3.15 Comparison of the PDF and CDF of the acyclic queueing network (The second station, 11-12-8-13). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 3.16 Comparison of the PDF and CDF of the acyclic queueing network (The second station, 6-8-12-13). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 x 4.1 The relationship between the cuto time and NSD (Doerr and Gue, 2006). . 75 4.2 The meaning of NSD based on the steady-state sojourn time distribution. . . 76 4.3 A change of ps of an order while in the system. The probability of success ps increases or decreases with time because the system state, the number of orders ahead, the number of servers, and the service time distribution change; but ps of an order always decreases in service because it is a conditional probability of elapsed time te. . . . . . . . . . . . . . . . . . . . . . . . . . . 81 4.4 Flushing policy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 4.5 Cascade policy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 4.6 The maximum ps that an order can reach by worker allocation: If we add 6 workers to shipping, the ps of the last order in the shipping queue increases to 39%. If we add 9 workers, then the ps of the last order reaches 60%. Adding more than 9 workers creates idle workers. . . . . . . . . . . . . . . . . . . . . 88 4.7 Best scenarios for all 12 Systems: Open circles represent good scenario; gray circles represent best scenario; empty space has the same NSD with the next lower ps. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 4.8 An approximate pattern for scenarios in Figure 4.7 . . . . . . . . . . . . . . 90 1 Best scenarios based on box and whisker chart of System 1 4 using SF . . . 112 2 Best scenarios based on box and whisker chart of System 5 8 using SF . . . 113 3 Best scenarios based on box and whisker chart of System 9 12 using SF . . 114 4 Best scenarios based on box and whisker chart of System 1 4 using MF . . . 121 5 Best scenarios based on box and whisker chart of System 5 8 using MF . . . 122 6 Best scenarios based on box and whisker chart of System 9 12 using MF . . 123 xi List of Tables 2.1 Shipping policies of America?s top 10 retail business. . . . . . . . . . . . . . 10 2.2 A comparison of the mean sojourn times between the intuitive model and simulation for a 3-workstation serial line, when all 1=C2 are integers. . . . . 20 2.3 A comparison of the mean sojourn times between the intuitive model and simulation when all 1=C2 are not integers. . . . . . . . . . . . . . . . . . . . 21 2.4 The weight coe cient of serial lines with 2, 3 and 4 workstations. . . . . . . 24 2.5 Mean sojourn times for the interpolation and simulation models of the serial line. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 2.6 Anderson-Darling tests for the serial line test case. . . . . . . . . . . . . . . . 29 2.7 Mean sojourn times for a 4-workstation general network. . . . . . . . . . . . 30 2.8 Anderson-Darling tests for the general queueing network. . . . . . . . . . . . 32 2.9 The system information of the Example 2.5 . . . . . . . . . . . . . . . . . . 32 2.10 Mean sojourn time for the Example 2.5. . . . . . . . . . . . . . . . . . . . . 32 2.11 The system information for Example 2.6. . . . . . . . . . . . . . . . . . . . . 33 2.12 Mean sojourn time for Example 2.6. . . . . . . . . . . . . . . . . . . . . . . . 34 3.1 A comparison of mean sojourn times for exponential service times. . . . . . . 47 3.2 Anderson-Darling tests for the exponential service time of single-stage queue. 48 3.3 A comparison of mean sojourn times for general service times. . . . . . . . . 57 3.4 Anderson-Darling tests for general service time of the single-stage queue. . . 57 3.5 The percent di erences of percentiles between NAM and Simulation for the single stage queue. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 3.6 The system information of the serial line and acyclic queueing network. . . . 62 xii 3.7 The comparison of the mean sojourn time of the serial line. . . . . . . . . . . 63 3.8 Anderson-Darling tests for the serial line. . . . . . . . . . . . . . . . . . . . . 64 3.9 The percent di erences of percentiles between the NAM and simulation for the serial line. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 3.10 The comparison of the mean sojourn time of the acyclic queueing network. . 65 3.11 Anderson-Darling tests for the acyclic queueing network. . . . . . . . . . . . 67 4.1 The characteristics of 12 systems. . . . . . . . . . . . . . . . . . . . . . . . . 83 4.2 System information for Systems 1, 5, and 9. . . . . . . . . . . . . . . . . . . 84 4.3 Results of the single ush policy for System 1. . . . . . . . . . . . . . . . . . 87 4.4 Change of NSD according to switching times ts. . . . . . . . . . . . . . . . . 89 4.5 Multi- ush results for System 1. . . . . . . . . . . . . . . . . . . . . . . . . . 92 4.6 Cascade policy results and comparison among three policies. . . . . . . . . . 95 4.7 The system performance according to the system characteristic. . . . . . . . 96 1 Single ush of the System 1 and 2 . . . . . . . . . . . . . . . . . . . . . . . 106 2 Single ush of the System 3 and 4 . . . . . . . . . . . . . . . . . . . . . . . 107 3 Single ush of the System 5 and 6 . . . . . . . . . . . . . . . . . . . . . . . 108 4 Single ush of the System 7 and 8 . . . . . . . . . . . . . . . . . . . . . . . . 109 5 Single ush of the System 9 and 10 . . . . . . . . . . . . . . . . . . . . . . . 110 6 Single ush of the System 11 and 12 . . . . . . . . . . . . . . . . . . . . . . 111 7 Multi- ush of the System 1 and 2 . . . . . . . . . . . . . . . . . . . . . . . . 115 8 Multi- ush of the System 3 and 4 . . . . . . . . . . . . . . . . . . . . . . . . 116 9 Multi- ush of the System 5 and 6 . . . . . . . . . . . . . . . . . . . . . . . . 117 10 Multi- ush of the System 7 and 8 . . . . . . . . . . . . . . . . . . . . . . . . 118 11 Multi- ush of the System 9 and 10 . . . . . . . . . . . . . . . . . . . . . . . 119 12 Multi- ush of the System 11 and 12 . . . . . . . . . . . . . . . . . . . . . . 120 xiii Chapter 1 Introduction 1.1 The order ful llment system An order ful llment system receives orders from points of sale in a continuous stream, processes them and delivers the nished goods to the customer at a speci c time or times during the day. Manufacturing and production systems and distribution centers are the most typical order ful llment systems. Consider a distribution center in which orders move through three stages | picking, packing and shipping. Order picking is the task of retrieving items from their storage lo- cations, and typically involves a worker traveling through the warehouse or an automated system retrieving the item and delivering it to the next stage via a conveyor or automated guided vehicle. Packing is the task of preparing items for delivery, and typically involves consolidating items in an order, packaging them for delivery and labeling the packages. Ship- ping is the task of staging orders and loading them onto trucks, and typically involves sorting orders (often on pallets) and building loads on outbound trucks. Orders arrive in this system throughout the day and are prepared for shipping through the picking and packing processes by workers. Figure 1.1 illustrates a nominal order ful llment system, where circles represent workers. Picking Packing Shipping Orders Completed Orders Figure 1.1: A nominal order ful llment system 1 In this system, we can view each activity (picking, packing and shipping) as a worksta- tion and workers in each workstation as servers. Thus we can model the order ful llment system as an open queueing network to analyze system performance. Customers expect a fair price, fast on-time delivery, and high quality in their purchase of products. Among those customer expectations, we focus on the delivery lead time because it is the most important function of the order ful llment system. Customers usually want to receive their orders as soon as possible. Thus, how companies respond to customer orders is a matter of primary concern for this system. In general, when we model a manufacturing or distribution center, we often use cycle time, throughput and WIP as system performance measures. That is, the objective of most research is to minimize cycle time, maximize throughput, or perhaps to minimize WIP. However, these metrics do not directly re ect customer satisfaction and two of them (throughput and WIP) are of no concern to the customer at all. Throughout our dissertation we use a performance metric called Next Scheduled Departure (NSD), which was introduced by Doerr and Gue (2006) to measure customer service in a distribution center. NSD records the percentage of orders arriving in the 24 hours between cuto times on successive days that make it onto the truck in the current day. The cuto time is established by managers in the distribution center, and any order arriving before this time is due in the current day. Doerr and Gue (2006) claim that an increase in the NSD means a direct increase in customer satisfaction. We might say that NSD is a system performance measure from the customer?s perspective, and existing metrics such as WIP, cycle time and throughput are from the manufacturer?s perspective. 1.2 Problem Statement Our research considers ways in which workers can be allocated dynamically in the few hours before the truck departure to \ ush" the system of orders that are almost completed, thereby increasing the service performance of the system. Our discussions with distribution 2 center managers suggest that they are increasingly considering cross trained workers to achieve operational exibility, as is often done in manufacturing. To implement worker allocation correctly, we need to answer the following questions: how many workers should we move and when? Simulation experiments on simple policies have led us to consider more complex, but intuitively appealing, policies based on sojourn time. The idea is to consider the likelihood (probability) that an order in the system will make it on the next departing truck if the allocation of workers does not change. When this probability drops below a threshhold value, then some workers in picking should be moved so as to improve the probability of success for other, more promising orders. To establish a sojourn-based policy, we must be able to compute the sojourn time distri- bution for an arriving order and for each order in the system. Note that this is di erent from, and more di cult than, computing the expected sojourn time. We require the distribution because we must make a probabilistic statement about the completion time of an order. For example, if an order arrives 3 hours before the truck departs, what is the probability that the order will make it on that truck? After a time lapse of 1 hour, what is the probability that the order will make it on that truck? If the probability is below the threshhold, then we \abandon" some orders and recon gure the system to help those with a higher chance of making it on time. To implement the dynamic worker allocation policies we propose in the dissertation, we need to know the probability that an order in the system will make it onto the next departing truck. Because we must be able to compute this probability at any time, we require the state-dependent distribution of sojourn time, where the state of the system is de ned by the number of busy servers and the number of orders in queue at any workstation. We develop this model in Chapter 3. As an intermediate step in the development, we develop an approximation model for the steady-state distribution of sojourn time for acyclic queueing networks with multiple servers. 3 In Chapter 2, we show how an order ful llment system might use this model to establish a cuto time that maximizes revenue for premium shipping. In a competitive setting, a rm wishes to set the cuto time as late as possible, in order to maximize the number of orders that could have the option to choose premium shipping; but not so late that the promised delivery cannot be made. We use our models to show how to make this tradeo . In the nal major section of the dissertation (Chapter 4), we use the state-dependent sojourn time distribution models from Chapter 3 to design dynamic worker allocation poli- cies to improve service performance according to the NSD metric. We test these policies extensively with simulation to show that, despite intuition to the contrary, it is possible to improve NSD by moving workers shortly before truck departure. We use the insights gained from this modeling to suggest a simple, easy-to-implement policy, which is also very e ective. 1.3 Contributions Contribution 1 We develop the rst approximation model for steady state sojourn time distributions in multi-server acyclic queueing networks. Our steady state sojourn time distribution model extends the work of Asmussen and M ller (2001) and You et al. (2002). Given the general interarrival times and general service times of the system, Asmussen and M ller (2001) compute waiting time distributions of multi- server single-stage queues, and You et al. (2002) approximate sojourn time distributions of single-server serial lines based on the characteristics of phase-type distributions. Because our system is a multi-server queueing network, we follow the method of Asmussen and M ller (2001) for the waiting time distribution of a single-stage with multiple servers, and convolve all waiting times and service times in the system using the method of You et al. (2002). Distributions produced by the model agree well with those produced by simulation for a variety of serial and general queueing networks. Research on steady state sojourn time distributions for queueing networks with general interarrival and service times has been addressed by Shanthikumar and Sumita (1988), You 4 et al. (2002) and Yoon (1994), but their studies were limited to the single-server case. Man- delbaum et al. (1998) develop di usion approximations for multi-server queueing networks, but their models assume exponential interarrival and service times. Of course, we can nd an approximation model for the G=G=c queueing network (Curry and Feldman, 2009), but such models are limited to the expected value of sojourn time. Contribution 2 We develop the rst approximation model for state-dependent sojourn time distributions in multi-server acyclic queueing networks. We develop our model based on the characteristics of phase-type distributions and a Markov processfN(t);t 0gto model the all-busy period for a waiting time distribution, where the state is de ned as the number of servers in each phase. The state-dependent sojourn time distribution is determined by the number of other orders ahead of an order, the number of servers, and the processing rate of the servers. In contrast to the steady state sojourn time distribution, we do not need to consider the arrival process for the state-dependent sojourn time distribution because it does not a ect the sojourn time distribution of an order already in the system. If there are k orders ahead of a particular order in a multi-server queueing system, it can enter service after k+1 di erent waiting times, so the state-dependent sojourn time distribution is the sum of all k+ 1 waiting times and the service time. We approximate each waiting time as a corresponding phase-type distribution based on the Markov process fN(t);t 0g, and nally all waiting times and the service time are convolved for the single- stage queue. We extend our work to acyclic G=G=c queueing networks based on the single-stage queue model. First, the state-dependent sojourn time distribution of each workstation is generated after revising the number of orders in front of it at each workstation based on the expected sojourn time of each workstation. And then all state-dependent sojourn time distributions are convolved for the queueing network. We compare distributions generated by the model with simulation. 5 Little research has been conducted on state-dependent sojourn time distributions for single-stage systems. Whitt (1999) and Nakibly (2002) proposed the method of estimating the waiting time distribution of each customer in a single-stage G=G=c queue. However, they focus on the waiting time distribution using a normal approximation, and their ap- proximation models perform well only for a small queue and a large number of servers. Our single-stage state-dependent sojourn time distribution model focuses on the sojourn time distribution, not waiting time distribution, and can be applied well to any condition. Our approximation model is also di erent from existing methods because we use the character- istics of the phase-type distribution. Contribution 3 We show how to allocate workers dynamically in an order ful llment sys- tem to improve its service performance. All of the policies we develop are based on the observation that, without dynamic worker allocation, some orders will just miss a departing truck, and that shifting workers (in a ware- housing context) from picking to shipping shortly before the last scheduled truck departure, could cause some of those orders to ship an entire day earlier than they otherwise would. An early question in our research was, \By shifting workers away from picking late in the day, would we hurt performance the next day by reducing capacity in picking?" Or, informally, are we \robbing Peter to pay Paul?" In Chapter 4 we show that, indeed, we are robbing Peter to pay Paul; But Peter still has enough cash: it is possible to improve performance| signi cantly, in some conditions|by moving the right number of workers at the right time everyday. We test three di erent policies that use sojourn time distributions, each for a wide range of parameters, and nd that the simplest policy performs just as well as more complicated ones. Generally speaking, it seems more important to do anything reasonable, rather than worrying about details of the allocation. We use this observation to construct a very simple policy based on a simple rule of thumb. This simple policy, performs nearly as well as the best policy, and it is much easier to implement in practice. 6 There is a rich literature on workers moving among tasks in a dynamic fashion in manufacturing and, to a lesser extent, warehousing. It is important to note that in all former papers that we know of the objective is to minimize cycle time, maximize throughput, or perhaps to minimize WIP, but these are not directly related to customer satisfaction. We use NSD as a system performance measure to establish the e ect of dynamic worker allocation, and the test results show that worker allocation can improve customer service. We believe an important feature of our research is the application of worker allocation models to systems that are due-date or deadline driven. 1.4 Organization of the dissertation The remainder of this dissertation is organized as follows. In the next chapter, we in- troduce an approximation model of steady state sojourn time distributions for multi-server acyclic queueing networks. In Chapter 3, we present an approximation model of state- dependent sojourn time distributions for multi-server single-stage queues and acyclic queue- ing networks. In Chapter 4, we describe dynamic worker allocation policies based on sojourn time distribution models which were developed in Chapters 2 and 3. In Chapter 5, we o er conclusions. 7 Chapter 2 An Approximation Model for Sojourn Time Distributions in Acyclic Multi-server Queueing Networks We develop an approximation model for the sojourn time distribution of customers or orders arriving to an acyclic multi-server queueing network. The model accepts general interarrival times and general service times, and is based on the characteristics of phase- type distributions. Distributions produced by the model agree well with those produced by simulation for a variety of serial and general queueing networks. 8 2.1 Introduction In many order ful llment systems, the rm makes a delivery promise based on the time of order receipt. Those ordering before a cuto time are \guaranteed" delivery by a particular mode. For example, Amazon?s web site famously asks if the customer wants his or her order the next day, and gives the condition that if so, the order must be placed within the next so many hours and minutes. This deadline counts down, in real time, as the customer makes a decision on mode of shipment. Typically, the current time plus this remaining time creates a cuto time of about 17:30 (at least in the author?s time zone), suggesting that after this time Amazon is no longer con dent it can make good on its promise. Amazon?s service is called \Guaranteed Accelerated Delivery" because if they miss the service promise, they refund the shipping cost. Promise-to-deliver service is bene cial to the customer because it reduces the uncertainty of delivery time. Most people have experienced waiting for a much anticipated package wondering, \Will it arrive today?" Promise-to-deliver service removes such anxiety. But notice that in such a system the service provider, in fact, manages the expectation of the customer by establishing the cuto time. Those ordering before the cuto time expect a certain delivery time, for sure, but those who just miss it have no expectation of that delivery. If the adage \service equals perception minus expectation" is true, this puts the retailer or distributor in control of perceived service with respect to delivery time|a very powerful position. Table 2.1 shows cuto times (if any) and shipping modes o ered by the ten largest internet retailers (Internetretailer, 2009). Six of the ten make an explicit promise on some type of shipping mode. For highly competitive industries such as books or electronics, a later cuto time gives the rm an advantage over its rivals; however, if the cuto time is too late, then a signi cant number of orders might not be processed in time to make the last package carrier pickup that evening (assuming next day service, for now). If the cuto time is too early, the rm loses both its competitive advantage and the pro t on customers who 9 Table 2.1: Shipping policies of America?s top 10 retail business. Firm Fastest Shipping Method Cuto time Compensation 1. Amazon.com Inc. One-day shipping 6:30 p.m. Refund 2. Staples Inc. Standard delivery 5 p.m. local None 3. Dell Inc. Next business day None None 4. O ce Depot Inc. Same day delivery 10 a.m. local None 5. Apple Inc. 2-Day Shipping 5 p.m. ET None 6. O ceMax Inc. Next business day 5 p.m. local None 7. Sears Holding Corp. Overnight air 6 p.m. CT None 8. CDW Corp. Next business morning None None 9. Newegg.com Next business day None None 10. Best Buy 2-Day Shipping None None (http://www.internetretailer.com, 2009) otherwise would have ordered premium shipping. How then to establish the cuto time as late as possible, but not \too late?" The tradeo we describe above is similar in spirit to the newsvendor problem, with a potential pro t (for premium shipping) and a penalty cost for non-performance. The penalty cost could be in terms of goodwill, or, in Amazon?s case, in direct monetary terms. We will assume the penalty cost can be explicitly determined by the rm. Still required for newsvendor-like analysis are the probabilities of performance and non-performance, which in our case require the expected sojourn time distribution for an arriving order. That is, if we know the distribution of sojourn time, then we can determine for any arriving order its probability of meeting or missing the service deadline. The challenge, then|and the major contribution of our chapter|is to develop a steady state sojourn time distribution for orders arriving to an order ful llment system, given in- formation about the order stream and processing times for workstations. After covering pre- liminaries on phase-type distributions in Section 3.2, we describe an approximation model for sojourn time distributions in Section 2.4 and 2.5. With the sojourn time distribution in hand, we can easily calculate the time at which the expected pro t for an order equals the expected loss, which then establishes the cuto time having the highest pro t to the rm. We validate our models with simulation in Section 2.5. 10 To summarize, we make two contributions in this chapter: First, we develop what we believe is the rst approximation model of sojourn time distributions for queueing networks with multiple servers per workstation and general interarrival and service times. The approx- imation we develop in this chapter is quite general, and could be used for any reasonably-sized system for which a queueing network is a reasonable model. Second, we show how an order ful llment rm might use such an approximation to establish a pro t-maximizing cuto time for promised, premium shipping. 2.2 Related literature Much of the literature on sojourn time modeling is focused on calculating the mean ow time (or simply the waiting time) in single-stage queues or in queueing networks (Whitt, 1983). Workstations in these systems may have one or many servers. The subject of our work is calculating the distribution of sojourn time, a more di cult task. We divide the literature on sojourn time modeling into four categories, according to whether the models address single-stage systems or networks, and whether they address single or multiple servers per workstation. The simplest case is single-stage, single-server systems. Neuts (1981) describes a matrix-geometric method to calculate the sojourn time distribution for a GI=G=1 system using phase type distributions. Luh and Zheng (2005) implemented Neuts? method in Mathematica. Sengupta (1989) used a bivariate Markov process to model waiting time and queue length distributions in a GI=PH=1 queue. His method is a \continuous analog" of the matrix-geometric method in Neuts (1981). Sengupta also showed that if the interarrival and service time distributions in a single-stage, single- server queue are phase-type, then the waiting time distribution is also phase-type. Asmussen and O?Cinneide (1998) verify the same result for a single-stage, multi-server queue, in a paper extending Sengupta?s work to the GI=PH=c case. Their model admits heterogeneous servers. Asmussen and M ller (2001) show how to calculate the waiting time distribution in a GI=PH=c and MAP=PH=c queue, for both homogeneous and heterogeneous 11 servers. Whitt (1999) also addresses single-stage, multi-server systems, but examines state dependent waiting time distributions. For example, he shows how to compute the waiting time distribution for the k-th customer in line in a c-server, single-stage system. Rueda (2003) develop an approximation for the waiting time distribution of single-stage queues with non-stationary interarrival and processing time distributions. Sojourn time distributions for queueing networks of single-servers have been addressed by Shanthikumar and Sumita (1988), You et al. (2002), and Yoon (1994). Shanthikumar and Sumita (1988) approximate the sojourn time distribution of an M=G=1 queueing system as one of three phase-type distributions (generalized Erlang, exponential, and hyperexponen- tial), based on a \service index," which is de ned as the squared coe cient of variation of the total service time of an arbitrary order. Neuts (1981) showed that the convolution of phase type distributions is also phase type. You et al. (2002) used this observation to show how to calculate the sojourn time distribution for queueing networks, with general interarrival and service times. In a paper published in Korean, Yoon (1994) developed a method very similar to that of You et al. Mandelbaum et al. (1998) develop di usion approximations for multi-server queueing networks, but their models assume exponential interarrival and service times. Missing from the literature are models of sojourn time distributions for queueing networks with multiple servers, when interarrival and service times can take general distributions. We ll that void here. We believe ours is the rst approximation model of sojourn time distributions for queueing networks of multi-server workstations. Our model is based on the work of Neuts (1981), Asmussen and M ller (2001), and You et al. (2002): We use the bivariate Markov process of Asmussen and M ller to extract the mean and variance of waiting times for each multi-server workstation. We then con- struct phase type distributions for waiting time and service time based on those means and variances. Then we use the method of You et al. to construct an in nitesimal generator and initial probability vector for the network. With these we can calculate a sojourn time 12 distribution for the network of multi-server queues. For reasons we discuss below, this model produces good results only when the inverse of the squared coe cient of variation (1=C2) is close to an integer. We correct this weakness with a simple but e ective interpolation scheme. We demonstrate the validity of our model by comparing it with results from a simulation model. In the next section, we provide some known results for phase type distributions. We introduce an intuitive model in Section 2.4. In Section 2.5, we develop an interpolation model and test it against a simulation model. We o er conclusions in Section 2.7. 2.3 Phase-type distributions Consider a Markov process on states f1;:::;m+ 1g, having in nitesimal generator T = 2 64 Q Q0 0 0 3 75; where Q is an m m matrix satisfying Qii < 0, for 1 i m, and Qij 0, for i6= j (For consistency, we follow the notation of Neuts, 1981). Q0 is a column vector of size m such that, Qe + Q0= 0; where 0 is a row vector of zeros and e is a column vector of ones. The initial probability vector of Q is = ( 1; 2;:::; m) such that e = 1. We are interested in the time until absorption into state m+1, whose distribution Neuts (1981) de nes as phase-type. Lemma 2.1 (Neuts, 1981) The probability distribution F( ) of the time until absorption in the state m+ 1, corresponding to the initial probability vector ( ; m+1) is given by F(x) = 1 eQxe: 13 Because F( ) is completely speci ed by and Q, the pair ( ;Q) is called a representation of F( ). Neuts (1989) provides the density function on (0;1), f(x) = eQxQ0 = eQx( Q)e: For each distribution of service or waiting times in a network, we seek a phase-type approximation for which we can generate a matrix-analytic model of the CDF. We use the fact that nite convolutions of phase-type distributions are also phase-type (Neuts, 1981) to generate solutions for a queueing network, as in You et al. (2002). The phase-type distribution is used to t a general distribution based on the in ow rate and the squared coe cient of variation C2 of a given positive random variable X. Sauer and Chandy (1975), You et al. (2002), and Tijms (1994) showed di erent tting methods that can convert a general distribution to a corresponding phase-type distribution based on the C2 of the distribution. A general distribution can be approximated with the Erlang distribution, the exponential distribution, or the hyperexponential distribution, based on its rst and second moments. We follow the rule of You et al. (2002): When C2 < 1, we convert a general distribution to an Erlang distribution, Erlang(m; ), of order m. The in nitesimal generator is Q = 2 66 66 66 64 0 0 0 0 0 0 ... ... ... ... ... ... 0 0 0 0 3 77 77 77 75 ; Q0 = 2 66 66 66 64 0 0 ... 3 77 77 77 75 ; where m = 1C2 and = m . The phase-type representation of a general distribution is represented by = (1;0;:::;0) and Q. 14 When C2 > 1, we use the hyperexponential distribution with balanced means, HE2, of order 2. We use the following normalization for this process, p 1 = q 2: The phase-type representation of a general distribution is given by = (p; q) and Q = 2 64 1 0 0 2 3 75; Q0 = 2 64 1 2 3 75; where p = 12 1 + q C2 1 C2+1 , q = 1 p, 1 = 2p and 2 = 2q . When C2 = 1, we have an exponential distribution (Tijms, 1994) represented by = 1 and Q, Q = ; Q0 = : 2.4 Intuitive model We now describe our model, which extends the work of Asmussen and M ller (2001) and You et al. (2002) to produce an approximation of the sojourn time distribution for a network of multi-server queues. For the waiting time distribution of a single-stage, multi- server system, we follow Asmussen and M ller?s method, then we convolve waiting times and service times in the system based on the method of You et al. The model requires only the following input, which should be available in most real systems: mean and variance of processing times for each workstation and mean and variance of interarrival times to the system. The procedure is: 1. Compute the arrival rate and SCV of an arrival process at each workstation using the Queueing Network Analyzer (QNA) method of Whitt (1983). 15 n the number of workstations in the general network i the total in ow rate into workstation i i the external in ow rate into workstation i pji the probability that a product ows from workstation j to i C2d(i) the SCV of the departure stream at workstation i C2a(i) the SCV of the interarrival stream at workstation i C2a(0;i) the SCV of external arrivals at workstation i C2s(i) the SCV of the service time at workstation i i the utilization of workstation i ci the number of servers at workstation i 2. Approximate an interarrival time and service time distribution of each workstation i as a corresponding phase-type distribution ( i;Ti) and ( i;Si) based on C2a and C2s. 3. Compute the mean and variance of waiting time of each workstation i using the method of Asmussen and M ller. Approximate the waiting time distribution of each worksta- tion i as a corresponding phase-type distribution ( i;Wi) based on C2w. 4. convolve all waiting times and service times sequentially to generate a phase-type representation of the sojourn time distribution ( ;K) using the convolution property of the phase-type distribution. 5. Solve F(t) = P(T t) = 1 eKte;(t 0) to obtain the cumulative distribution function (CDF) of the sojourn time distribution. 2.4.1 The ow rate and SCV of the arrival process We adopt the notation of Curry and Feldman (2009): To compute the waiting time distribution for each workstation, we must approximate the interarrival time and service time distribution as a corresponding phase-type distribution. Thus, the ow rate i and C2a(i) of an interarrival stream at each workstation i must be calculated rst, because the tting distribution is based on the rate and C2 of a general distribution. 16 We compute the in ow rate and SCV of the interarrival stream at each workstation i using Whitt?s approximations. The total in ow rate into workstation i is i = i + nX j=1 pji j: The SCV of the departure stream is approximately, C2d(i) = 1 + (1 2i)(C2a(i) 1) + 2i (C 2 s(i) 1)p ci : Finally, the SCV of an interarrival stream is approximately, C2a(i) = i i C2a(0;i) + nX j=1 jpji i [pjiC 2 d(j) + 1 pji]: From these expressions, we can approximate the interarrival time distribution and ser- vice time distribution of each workstation i as a corresponding phase-type distribution ( i;Ti) and ( i;Si) based on the ow rate and SCV. 2.4.2 Waiting time for a single stage queue with multiple servers Asmussen and M ller (2001) proposed a method to compute the steady-state waiting time distribution in a multi-server queue, with each server having a phase-type service time distribution. They used matrix-analytic methods based on the method of Sengupta (1989) and Asmussen and O?Cinneide (1998). Asmussen and O?Cinneide (1998) showed that the GI=PH=c waiting time is always phase-type, and that the number of phases for the homo- geneous case should be 0 B@ m+c 1 c 1 CA; where m is the number of phases of the service time distribution. 17 They solve the system using a bivariate Markov process fXt;Nt;t 0g to model the all busy period, where Xt is the time since arrival of the last customer to enter service in the all-busy-period, and Nt is the current phase in which the server is working. They proved that the steady state density of the Markov process fXt;Ntg is matrix exponential: (x) = (0) exp(Tx); x> 0; where T is the rate matrix. Their procedure is as follows: The rate matrix T satis es the non-linear integral equation T = S + Z 1 0 exp(Tu)A(jump)H(du); where the matrix S contains transition rates without service completion, A(jump) contains transition rates with service completion in fNtg, and H(du) is the distribution function of the interarrival times. The matrix T is solved by iterative procedure and the phase-type representation of the waiting time of the single-stage is de ned by ( ; G); i = i?i ?; Gij = jTji i ; where = (0)T 1 and ? = (T S)e. Finally, the mean and variance of the waiting time distribution are computed from the rst and second moment, where the n-th moment is E(Xn) = ( 1)nn! G ne: 18 2.4.3 Approximation of sojourn time distribution of the queueing network Now we are ready to apply the convolution property of phase-type distributions. Neuts (1981) showed that if F( ) and G( ) are both continuous phase-type distributions with repre- sentations ( ;T) and ( ;S) of orders m and n, then the convolution F G( ) is a phase-type distribution with representation ( ; L), where the in nitesimal generator L and initial prob- ability vector are L = 2 64 T T0 0 S 3 75; = [ ; m+1 ]: We construct the initial probability vector and the in nitesimal generator of a queueing network by convolving all waiting times ( i;Wi), and service times ( i;Si) by order in the network. For example, for a k station serial line, the in nitesimal generator K and initial probability vector are K = 2 66 66 66 66 66 4 W1 W01 1 0 0 0 0 S1 S01 2 0 0 ... ... ... ... ... ... 0 0 0 Wk W0k k 0 0 0 0 Sk 3 77 77 77 77 77 5 ; =[ 1 0 0 ]: The cumulative distribution function (CDF) and the probability density function (PDF) are F(t) = P(T t) = 1 eKte;(t 0); f(t) = eKtK0 = eKt( K)e;(t 0): 19 Table 2.2: A comparison of the mean sojourn times between the intuitive model and simu- lation for a 3-workstation serial line, when all 1=C2 are integers. E[T] C2 Mean sojourn time (hrs) % di erenceSimulation Intuitive Interarrival 0.5 0.5 5.79 5.77 -0.23Workstation 1 1.8 0.5 0.6Workstation 2 2.2 0.5 0.73 Workstation 3 1.5 0.5 0.5 In a more general queueing network, ow between workstations depends on product routings. We estimate the sojourn time distribution for a random order by approximating the CDF of all possible \serial lines" (paths in the network) and mixing those CDFs according to the probabilities of taking those paths. Because our model is based on serial line analysis, we are able to model only acyclic queueing networks. We give an example in a later section. 2.4.4 Numerical results To validate the intuitive model, we compare it with results of a simulation built in Arena 10.0. We assume processing and interarrival times in the simulation model are Gamma distributed. Example 2.1 Consider a serial line of 3 workstations, with 6 servers per workstation. Information on the interarrival and processing times are given in Table 2.2. We can approximate the interarrival time and three service time distributions as Erlang(m; ) distributions according to the rule of You et al. (2002) because all C2 are less than 1. Note that in this case all 1C2 and 1=C2 have the same value of 2. The mean sojourn times in Table 2.2 show that the intuitive model is close to the simulation result. Figure 2.1 shows the sojourn time distributions from the intuitive model and from the simulation, and we see nice agreement between them. 20 5 10 15 20 Time 0.025 0.05 0.075 0.1 0.125 0.15 0.175 Density PDF 5 10 15 20 Time 0.2 0.4 0.6 0.8 1 Probability CDF Intuitive Simulation Figure 2.1: Sojourn time distributions for Example 1, when all 1=C2 are integers. Table 2.3: A comparison of the mean sojourn times between the intuitive model and simu- lation when all 1=C2 are not integers. E[T] C2 Mean sojourn time (hour) % di erenceSimulation Intuitive Interarrival 0.5 0.8 5.97 5.77 -3.35Workstation 1 1.8 0.6 0.6Workstation 2 2.2 0.5 0.73 Workstation 3 1.5 0.9 0.5 Example 2.2 Consider the same 3-workstation serial line with C2a = 0:8, C2s(1) = 0:6, C2s(2) = 0:5 and C2s(3) = 0:9. The mean sojourn time and sojourn time distribution of this system are shown in Table 2.3 and Figure 2.2 respectively. The percent di erence of mean sojourn time between the two models is much higher, and the distributions appear not to be so similar. The intuitive model produces the same result when 1=C2 falls into the interval 1C2 ; 1C2 , but we get a new simulation result with a new C2. For example, although C2a changed from 0.5 to 0.8, the interarrival time distribution is still tted with the same Erlang(m; ) distribution, because 1C2 is the same for both cases. This example illustrates a characteristic of the model that we have con rmed with more thorough testing: the intuitive model produces good results only when 1=C2 is close to an integer value (under the C2 < 1 condition). When 1=C2 is between integer values (say, for example, when C2 = 0:75), the intuitive model produces less satisfactory results. In fact, 21 5 10 15 20 Time 0.025 0.05 0.075 0.1 0.125 0.15 0.175 Density PDF 5 10 15 20 Time 0.2 0.4 0.6 0.8 1 Probability CDF Intuitive Simulation Figure 2.2: A comparison of sojourn time distributions between the intuitive model and simulation when all 1=C2 are not integers. the model produces the same sojourn time distribution for all values of C2 having a common value of 1C2 . 2.5 Interpolation model A na ve attempt to correct for this de ciency by replacing the ceiling operator with a round( ) function produced even worse results. We propose, then, the following interpola- tion scheme: 1. Produce CDF F(t) using 1C2 phases for each workstation. 2. Produce CDF G(t) using 1C2 phases for each workstation. 3. Compute the mixed CDF, H(t) = G(t) + (1 )F(t), where is an interpolation coe cient. The interpolation coe cient is based on two observations. First, if 1=C2 for interarrival and processing times is close to 1C2 , then we should expect to be close to zero, and if 1=C2 1C2 , should be close to one. Second, interarrival stream and workstation processing times will have di erent levels of in uence on the sojourn time distribution, depending on whether they appear early or late in the (serial) process. These observations lead us to a regression model of the squared coe cients of variation for each distribution in the process. To address the rst observation, we de ne a tting 22 coe cient for each interarrival and processing time distribution. The tting coe cient is an interpolation coe cient re ecting the nearness of the squared coe cient of variation to successive values of C2 when 1=C2 takes on an integer value. For example, if the SCV of the arrival process C2a = 0:8, we compute a tting coe cient 0 = 0:6, because 0:5+0:6(1 0:5) = 0:8, or said another way, 0.8 is 60% of the way between 0.5 and 1. We can compute tting coe cients k for each workstation k in a similar fashion. Next, we determine the di erent levels of in uence of the tting coe cients with a regression model. We use the intuitive model to establish the dependent variable (mean sojourn time) under possible combinations of SCVs, which function as the independent variables. We test the intuitive model for all possible combinations of three SCVs (1, 1/2, 1/3), each corresponding to an integer value of 1=C2. For example, the regression for two workstation serial lines must consider three distributions|the interarrival times and service times for two workstations, and there are 33 = 27 combinations. We excluded SCVs less than 1/3 because they have very little e ect on mean sojourn time. There is one last detail: we must consider the possible e ects of utilization on the mean sojourn times and therefore on the coe cients !i produced by the regression. We ran two cases ( = 0:5; 0:85 for each workstation in the system), and found that the resulting coef- cients !i are di erent. Which values to use, given that we seek a single set of coe cients for any serial line of a speci ed length? We answered this question by using for each distri- bution the average of the two values, for two reasons: (1) For most systems in practice, the utilization of workstations in a system are di erent anyway, and (2) the results of the model are not very sensitive to the precise values of the !i?s. Table 2.4 contains results of the regression analysis, which are the weight coe cients !i for 2, 3, and 4 station serial lines. In general, the interarrival stream and the rst workstation account for the largest in uence on the sojourn time distribution and the in uence decreases as it goes downstream. 23 Table 2.4: The weight coe cient of serial lines with 2, 3 and 4 workstations. Workstation !0 !1 !2 !3 !4 2 0.493 0.313 0.194 3 0.412 0.292 0.174 0.123 4 0.353 0.262 0.180 0.123 0.081 Now that we have the tting coe cients i and weight coe cients !i. We calculate the interpolation coe cient = kX i=0 !i i; where index zero represents interarrival times. We should emphasize the conditions for which the !i?s in Table 2.4 are relevant: These coe cients should be appropriate for any serial lines of these lengths for which the Gamma distribution is a reasonable approximation of processing times. Before fully testing the interpolation model, we need to show that it is better than the simpler, intuitive model. Figure 2.3 shows the di erences in mean sojourn times among the simulation and the two approximation models when each workstation has the same utilization =0.5 and SCVs. The interpolation and intuitive models are the same when all 1=C2 are integers, but when they are not, the interpolation model is closer to the simulation. The intuitive model also corrects the de ciency noted in Table 2.3: for this problem, the di erence in means is reduced to -1.26% (mean sojourn time = 5.90 hours). 2.5.1 Validation To test the performance of the model, we consider four queueing networks: The rst is a small serial line with 3 workstations (Figure 2.4), the second a simple acyclic queueing network with 4 workstations (Figure 2.5), the third is a large-scale system with almost 100 servers in 3 workstations, and the last is a long serial line with 8 workstations. We wish to con rm the e ectiveness of the interpolation model on a basic serial line and general queueing network through the rst two experiments, and wish to know computation time and accuracy 24 0.2 0.4 0.6 0.8 1 C 2 4.52 4.54 4.56 4.58 4.60 4.62 4.64 Mean sojourn time Simula tion In ter pola tion In tuitiv e Figure 2.3: A comparison of mean sojourn times among three models( =0.5). Dots indicate data; lines are inserted for visual clarity only. 1 2 3 Figure 2.4: A serial line with 3 workstations. of the interpolation model as we increase the number of servers and workstations. In Figures 2.4 and 2.5, the black circle in each workstation represents a server, and in the case of the general network, a customer or order departing from the rst workstation selects its follow on workstation with probability p. Example 2.3 Consider a 3-workstation serial line with 6 servers per workstation (Figure 2.4). All servers have the same capacity (E[T] = 1:5), and interarrival times are adjusted to produce the appropriate utilization ( = 0:5; 0:85). We compare simulation and model results over di erent values of SCV for interarrival and processing times. Table 2.5 shows the mean sojourn times for the interpolation model and simulation. 25 0.67p= 0.33p= 1 2 3 4 Figure 2.5: A general network with 4 workstations. Table 2.5: Mean sojourn times for the interpolation and simulation models of the serial line. C2a C2s(1) C2s(2) C2s(3) Interpolation Simulation % di erence 0.85 1 1 1 1 7.63 7.57 0.83 0.75 0.75 0.75 0.75 6.77 6.77 -0.04 0.57 0.7 0.6 0.9 6.44 6.43 0.05 0.45 0.75 0.4 0.9 6.13 6.19 -0.35 0.4 0.6 0.26 0.8 5.88 5.89 -0.13 0.33 0.33 0.33 0.33 5.57 5.48 1.71 0.5 1 1 1 1 4.65 4.62 0.51 0.75 0.75 0.75 0.75 4.53 4.58 -0.30 0.57 0.7 0.6 0.9 4.57 4.55 0.29 0.45 0.75 0.4 0.9 4.54 4.54 0.04 0.4 0.6 0.26 0.8 4.53 4.53 -0.03 0.33 0.33 0.33 0.33 4.52 4.53 -0.23 26 5 10 15 20 Time 0.025 0.05 0.075 0.1 0.125 0.15 Density PDF 5 10 15 20 Time 0.2 0.4 0.6 0.8 1 Probability CDF Interpolation Simulation Figure 2.6: CDF and PDF of the serial line ( =0.85, C2 = 0:75; 0:75; 0:75; 0:75). The percent di erence of the mean sojourn time between the two models under high utilization is a little bit higher than in the low utilization case. We conjecture that this is due to the amount of waiting time at each workstation. In the low utilization case, there is very little waiting in front of each workstation, so inaccuracies in the waiting time model do not have much e ect on the estimate of the mean sojourn time. In the high utilization case, model inaccuracies have a greater e ect. However, we see that the percent di erences between the interpolation model and the simulations are still less than 2 percent. Notice that this statement only pertains to the mean, when what we seek is the distribution; that is, agreement in the means is only an indication that the model is e ective, not a proof. Figures 2.6{2.11 show the distributions for representative cases. We use the Anderson- Darling test (Law and Kelton, 2000) to check the graphical agreement of the distributions between simulation and approximation results. We also include the percent di erence of the 90th and 95th percentiles. As shown in Table 3.8, we accept that there is no di erence between the distributions based on the A-D test, and the percent di erences of the 90th and 95th percentiles are less than 3 percent. Figures 2.6{2.11 indicate that the model consistently has a higher peak in the PDF than does the simulation model. We believe this is due to our approximation of waiting time distributions, which for all cases we tested had C2 > 1, and therefore led us to the hyper- exponential distribution as an approximation. Because the hyperexponential distribution is 27 5 10 15 20 Time 0.025 0.05 0.075 0.1 0.125 0.15 0.175 Density PDF 5 10 15 20 Time 0.2 0.4 0.6 0.8 1 Probability CDF Interpolation Simulation Figure 2.7: CDF and PDF of the serial line ( =0.85, C2 = 0:57; 0:7; 0:6; 0:9). 5 10 15 20 Time 0.05 0.1 0.15 0.2 Density PDF 5 10 15 20 Time 0.2 0.4 0.6 0.8 1 Probability CDF Interpolation Simulation Figure 2.8: CDF and PDF of the serial line ( =0.85, C2 = 0:4; 0:6; 0:26; 0:8). 5 10 15 20 Time 0.05 0.1 0.15 0.2 Density PDF 5 10 15 20 Time 0.2 0.4 0.6 0.8 1 Probability CDF Interpolation Simulation Figure 2.9: CDF and PDF of the serial line ( =0.5, C2 = 0:75; 0:75; 0:75; 0:75). 5 10 15 20 Time 0.05 0.1 0.15 0.2 Density PDF 5 10 15 20 Time 0.2 0.4 0.6 0.8 1 Probability CDF Interpolation Simulation Figure 2.10: CDF and PDF of the serial line ( =0.5, C2 = 0:57; 0:7; 0:6; 0:9). 28 5 10 15 20 Time 0.05 0.1 0.15 0.2 Density PDF 5 10 15 20 Time 0.2 0.4 0.6 0.8 1 Probability CDF Interpolation Simulation Figure 2.11: CDF and PDF of the serial line ( =0.5, C2 = 0:4; 0:6; 0:26; 0:8). Table 2.6: Anderson-Darling tests for the serial line test case. C2a C2s(1) C2s(2) C2s(3) A-D = 5% Decision % di erence90th 95th 0.85 0.75 0.75 0.75 0.75 1.691 2.492 accept 0.44 -0.32 0.57 0.7 0.6 0.9 1.726 accept 0.25 -0.30 0.4 0.6 0.26 0.8 1.529 accept -1.15 -0.38 0.5 0.75 0.75 0.75 0.75 0.639 2.492 accept 1.92 -0.26 0.57 0.7 0.6 0.9 1.366 accept -2.35 -1.99 0.4 0.6 0.26 0.8 1.542 accept -1.65 -1.88 skewed to the left, it may cause our distributions of sojourn time also to be overly skewed to the left. Example 2.4 Consider the general acyclic queueing network with 4-workstations in Figure 2.5. Workstations 1 and 4 have 6 servers, and workstations 2 and 3 have 4 and 2 servers respectively. All servers have the same capacity (E[T] = 1:5), and interarrival times are adjusted to produce the appropriate utilization ( = 0:5; 0:85). Table 2.7 shows the mean sojourn time of the interpolation model and the simulation model in a 4-workstation general queueing network. Similar to the serial line case, the interpolation model gives better results under low utilization than under high utilization. The percent di erences are all less than 5 percent. Figures 2.12{2.17 show distributions from the model and simulation. The Anderson- Darling tests are all acceptable and values for the 90th and 95th percentiles show fairly nice agreement. 29 Table 2.7: Mean sojourn times for a 4-workstation general network. C2a C2s(1) C2s(2) C2s(3) C2s(4) Interpolation Simulation % di erence 0.85 1 1 1 1 1 8.99 8.91 0.89 0.75 0.75 0.75 0.75 0.75 7.78 7.98 -2.61 0.57 0.7 0.6 0.6 0.9 7.31 7.49 -2.42 0.45 0.75 0.4 0.4 0.9 6.73 7.05 -4.48 0.4 0.6 0.26 0.26 0.8 6.38 6.68 -4.46 0.33 0.33 0.33 0.33 0.33 6.11 6.27 -2.62 0.5 1 1 1 1 1 4.76 4.81 -1.11 0.75 0.75 0.75 0.75 0.75 4.73 4.74 -0.17 0.57 0.7 0.6 0.6 0.9 4.68 4.69 -0.36 0.45 0.75 0.4 0.4 0.9 4.60 4.66 -1.12 0.4 0.6 0.26 0.26 0.8 4.59 4.64 -0.94 0.33 0.33 0.33 0.33 0.33 4.58 4.62 -0.84 5 10 15 20 25 30 Time 0.02 0.04 0.06 0.08 0.1 0.12 0.14 Density PDF 5 10 15 20 25 30 Time 0.2 0.4 0.6 0.8 1 Probability CDF Interpolation Simulation Figure 2.12: CDF and PDF of the general network( =0.85,C2 = 0:75;0:75;0:75;0:75;0:75). 5 10 15 20 25 30 Time 0.025 0.05 0.075 0.1 0.125 0.15 Density PDF 5 10 15 20 25 30 Time 0.2 0.4 0.6 0.8 1 Probability CDF Interpolation Simulation Figure 2.13: CDF and PDF of the general network ( =0.85, C2 = 0:57; 0:7; 0:6; 0:6; 0:9). 30 5 10 15 20 25 30 Time 0.025 0.05 0.075 0.1 0.125 0.15 0.175 Density PDF 5 10 15 20 25 30 Time 0.2 0.4 0.6 0.8 1 Probability CDF Interpolation Simulation Figure 2.14: CDF and PDF of the general network ( =0.85, C2 = 0:4; 0:6; 0:26; 0:26; 0:8). 5 10 15 20 Time 0.05 0.1 0.15 0.2 Density PDF 5 10 15 20 Time 0.2 0.4 0.6 0.8 1 Probability CDF Interpolation Simulation Figure 2.15: CDF and PDF of the general network ( =0.5,C2 = 0:75;0:75;0:75;0:75;0:75). 5 10 15 20 Time 0.05 0.1 0.15 0.2 Density PDF 5 10 15 20 Time 0.2 0.4 0.6 0.8 1 Probability CDF Interpolation Simulation Figure 2.16: CDF and PDF of the general network ( =0.5, C2 = 0:57; 0:7; 0:6; 0:6; 0:9). 5 10 15 20 Time 0.05 0.1 0.15 0.2 Density PDF 5 10 15 20 Time 0.2 0.4 0.6 0.8 1 Probability CDF Interpolation Simulation Figure 2.17: CDF and PDF of the general network ( =0.5, C2 = 0:4; 0:6; 0:26; 0:26; 0:8). 31 Table 2.8: Anderson-Darling tests for the general queueing network. C2a C2s(1) C2s(2) C2s(3) C2s(4) A-D = 5% Decision % di erence90th 95th 0.85 0.75 0.75 0.75 0.75 0.75 2.343 2.492 accept -4.02 -1.33 0.57 0.7 0.6 0.6 0.9 2.401 accept -2.94 0.22 0.4 0.6 0.26 0.26 0.8 2.046 accept -5.70 -5.27 0.5 0.75 0.75 0.75 0.75 0.75 1.475 2.492 accept -3.57 -3.28 0.57 0.7 0.6 0.6 0.9 2.086 accept -3.98 -4.33 0.4 0.6 0.26 0.26 0.8 1.454 accept -2.93 -3.05 Table 2.9: The system information of the Example 2.5 E[T] C2s C2a Number of workers Interarrival 0.23 0.75 Workstation 1 7.80 0.70 0.85 40 Workstation 2 4.50 0.70 0.78 25 Workstation 3 5.80 0.80 0.84 30 Example 2.5 Consider a serial line with 3-workstations, and 95 workers in the system. We assume that both the interarrival and service times follow general distributions. Table 2.9 shows the system information in detail. We compare the mean sojourn time and distribution between the interpolation and simulation model. We also record computation time using a 2.4 GHz Intel Core 2 Duo processer, and it takes 1.55 hours. Table 2.10 and Figure 2.18 show mean sojourn time and sojourn time distribution. We see that the interpolation model also works well for the large scale system, but the computation time has increased signi cantly. As mentioned before, the size of the in- nitesimal generator K and the initial probability vector , which decide the sojourn time distribution, are determined by the size of the matrix of the waiting time distribution W Table 2.10: Mean sojourn time for the Example 2.5. Interpolation % di erence Simulation % di erence90th 95th Mean sojourn time 18.60 0.30 18.54 0.01 0.01 32 20 40 60 80 Time 0.01 0.02 0.03 0.04 0.05 Density PDF 20 40 60 80 Time 0.2 0.4 0.6 0.8 1 Probability CDF Interpolation Simulation Figure 2.18: PDF and CDF of the Example 2.5. Table 2.11: The system information for Example 2.6. E[T] C2s C2a Number of workers Interarrival 0.23 0.75 Workstation 1 2.90 0.70 0.84 15 Workstation 2 2.50 0.60 0.91 12 Workstation 3 1.40 0.80 0.87 7 Workstation 4 1.50 0.65 0.82 8 Workstation 5 2.75 0.50 0.91 13 Workstation 6 2.20 0.85 0.87 11 Workstation 7 3.10 0.70 0.79 17 Workstation 8 1.50 0.60 0.82 8 and service time distribution S. Here, W is determined by 0 B@ m+c 1 c 1 CA, where m is the number of phases of the service time distribution and c is the number of servers, and S is determined by only m. So if we increase c, it a ects mainly to the size of matrix W, and this creates a signi cant computational burden. Example 2.6 Consider a serial line with 8-workstations, and 91 workers in the system. We assume that both the interarrival and service times follow general distributions. Table 2.11 shows the system information in detail. We compare the mean sojourn time and distribution between the interpolation and simulation model. Table 2.12 and Figure 2.19 show mean sojourn time and sojourn time distribution, but the results are not good. The percent di erence in the means is greater 33 Table 2.12: Mean sojourn time for Example 2.6. Interpolation % di erence Simulation % di erence90th 95th Mean sojourn time 20.91 -7.71 22.65 7.51 -6.31 10 20 30 40 50 Time 0.01 0.02 0.03 0.04 0.05 0.06 0.07 Density PDF 10 20 30 40 50 Time 0.2 0.4 0.6 0.8 1 Probability CDF Interpolation Simulation Figure 2.19: PDF and CDF for the long serial line in Example 2.6. than 5%, and the percentiles of 90th and 95th also are greater than 5%. However, it takes just about 13 minutes because the number of servers in each workstation is not too large (around 15), and matrix size of each waiting time distribution is of reasonable size. 2.6 Determining an optimal cuto time Consider an order ful llment center with three major functional areas|picking, packing, and shipping, all of which function as multi-server queues. Assume that each functional area processes orders individually. For most warehouses, this is a reasonable assumption for the packing and shipping areas, but for manual warehouses, picking is often done in batches. We assume for purposes of illustration a suitable part-to-picker automated system for which a multi-server queue is an appropriate model. Nearly every order moves from picking to packing to shipping, sequentially, so we can model the process as a serial line, with multiple workers in each stage. Assume the rm realizes marginal shipping pro t r for an order successfully delivered for premium shipping. If the rm promises premium service but fails to deliver, it realizes marginal cost c. For most applications, we would expect r c (in Amazon?s case, c is 34 Optimal cutoff t c * Deadline t d Delivery date Shipping time p s 1-p s Arriving time t a Figure 2.20: Relationship between cuto time and probability of success of an arriving order. the entire cost of shipping). We seek a cuto time t c, which gives the maximum expected shipping pro t. Let probability of success ps be the probability that an arriving order will make it on the next departing truck. We can compute ps of an arriving order from the steady state sojourn time distribution. Let a continuous random variable T be the sojourn time of an order. Then, ps = P[T