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 1, a general distribution is converted to a balanced hyper-exponential distribu- tion, HE2. The density function of HE2 is described by f(x) = p1 1e 1x +p2 2e 2x;x 0; where p1 = 12 1 + q C2 1 C2+1 , p2 = 1 p1, 1 = 2p1=E[X] and 2 = 2p2=E[X]. If C2 = 1, we use the exponential distribution, with density function, f(x) = e x;x 0; 44 where = 1=E[X]. 3.3 Exponential service times In this section, we assume all servers in a workstation are identical and have the same exponential distribution of processing times. Whitt (1999) suggested that the waiting time distribution in such a system follows an Erlang-(k + 1) distribution with mean 1=c when the system has c servers and there are k customers ahead. For now, we are interested in the complete sojourn time distribution (waiting plus processing), so we develop an approximation model for multi-server queueing systems with exponential service times. 3.3.1 Approximation model The state dependent sojourn time distribution does not rely on the interarrival time distribution because the waiting time distribution is determined only by the number of servers, the number of orders ahead and the service rate. Also, we do not need to consider the remaining service time for exponential service, due to the memoryless property. Every order arriving to the system has Erlang waiting time in queue and is processed in exponential service time. We approximate the sojourn time distribution of an order in the system based on the convolution property of the phase-type distribution, because the waiting times and service time are phase-type. Neuts (1981) introduced the convolution property of the phase-type distribution. If H( ) and I( ) are both continuous phase-type distributions with representa- tions ( ;T) and ( ;S) of orders m and n, then the convolution of two distributions, H I( ) is also a phase-type distribution with representation ( ;Q) and the in nitesimal generator and the initial probability vector are given by Q = 2 64 T T0 0 S 3 75; 45 = [ ; m+1 ]: We construct the initial probability vector and the in nitesimal generator of the system based on this result. The in nitesimal generator is composed of the waiting time represen- tation ( ;W) and service time representation ( ;S). The sizes of the matrix W and vector are determined by the number of orders ahead in the queue. For example, if there are k orders ahead, the size of the in nitesimal generator and initial probability vector of waiting time are given by W(k+1) (k+1) and 1 (k+1). Thus, the phase-type representation of the sojourn time distribution is determined by K = 2 64 W W0 0 S 3 75; = [1;0; ;0]: We can generate the cumulative distribution function (CDF) of the system based on the phase-type representation ( ;K) F(t) = P(T t) = 1 eKte;(t 0): The probability density function (PDF) is given by f(t) = eKtK0 = eKt( K)e;(t 0): 3.3.2 Numerical results To test the approximation model, we consider two factors|the number of servers c and the number of orders ahead k. We use c =f 2, 3, 5, 10, 20, 50, 100, 200g, and k =f5, 10, 20g for c 100, f40, 60, 80g for c = 200. 46 Table 3.1: A comparison of mean sojourn times for exponential service times. Servers 2 3 5 Orders ahead 5 10 20 5 10 20 5 10 20 Approximation 20.00 32.50 57.28 15.00 23.33 40.00 11.00 16.00 26.00 % di erence 0.08 0.20 -0.01 -0.63 -0.56 -0.44 0.52 0.59 0.58 Simulation 19.98 32.44 57.28 15.09 23.46 40.18 10.94 15.91 25.85 Servers 10 20 30 Orders ahead 5 10 20 5 10 20 5 10 20 Approximation 8.00 10.50 15.50 6.50 7.75 10.25 5.99 6.81 8.47 % di erence 0.17 0.31 0.36 -0.04 -0.01 -0.01 0.04 0.08 0.17 Simulation 7.99 10.47 15.44 6.50 7.75 10.25 5.98 6.80 8.46 Servers 50 100 200 Orders ahead 5 10 20 5 10 20 40 60 80 Approximation 5.61 6.09 7.10 5.29 5.56 6.04 6.01 6.54 7.01 % di erence 0.34 -0.04 -0.06 -0.05 0.31 -0.18 -0.18 0.30 -0.14 Simulation 5.59 6.10 7.10 5.30 5.55 6.05 6.02 6.52 7.02 We test possible combinations of these two factors under the same capacity (E[Ts] = 5) and compare with the results of a simulation built in Arena 10.0. We assume processing time in the simulation model is exponentially distributed. Table 3.1 shows the mean sojourn time results from the approximation model and the simulation model. Di erences in the means are all less than 1 percent. Figures 3.1{3.3 show the PDF and CDF of the approximation and simulation results. As with the mean sojourn time, we see nice agreement. To make a more scienti c statement, we chose three cases and checked the agreement using the Anderson-Darling (A-D) test. We nd that all statistics of the tests are less than the critical value with = 5% (see Table 3.2). The percent di erences of the 90th and 95th percentile between the two models are also small. 3.4 General service times In addition to the exponential service time case, Whitt (1999) developed an approxima- tion model for general service times in single stage systems. For a large number of waiting 47 20 40 60 80 100 Time 0.005 0.01 0.015 0.02 0.025 0.03 Density PDF 20 40 60 80 100 Time 0.2 0.4 0.6 0.8 1 Probability CDF Approximation Simulation Figure 3.1: Comparison of the PDF and CDF of the approximation model and the simulation model (2 servers and 20 orders ahead). 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 Approximation Simulation Figure 3.2: Comparison of the PDF and CDF of the approximation model and the simulation model (50 servers and 10 orders ahead). 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 Approximation Simulation Figure 3.3: Comparison of the PDF and CDF of the approximation model and the simulation model (100 servers and 5 orders ahead). Table 3.2: Anderson-Darling tests for the exponential service time of single-stage queue. Servers Orders ahead A-D = 5% Decision % di erence90th 95th 2 20 1.405 2.492 accept 0.66 2.44 50 10 1.329 accept 3.10 0.31 100 5 0.631 accept 0.82 0.17 48 customers k, a Normal approximation can be used. The mean waiting time is given by E[W] k + 1 c 1 + 12c ; and the full distribution of waiting time is described by a Normal approximation. Whitt states that his approximation is appropriate only when the number of servers is substantially larger than the number of customers ahead, as is commonly the case in call centers. In contrast to Whitt (1999), we approximate the sojourn time distribution using phase- type distributions. The procedure is similar to the exponential service times case, except for considering the remaining service times of the orders in service. The procedure is: 1. Approximate the service time distribution as a corresponding phase-type distribution ( ;S) based on the C2s. 2. Approximate the rst waiting time distribution as a corresponding phase-type distri- bution ( 1;W1) based on the Markov process. 3. Approximate the second waiting time distribution as a corresponding phase-type dis- tribution ( 2;W2) based on the Markov process and the rst initial probability vector 1. 4. Approximate successive waiting time distributions as corresponding phase-type distri- butions ( i;Wi), i 3, according to the same procedure. 5. Generate the initial probability vector and in nitesimal generator ( ;K) for the system using the convolution property of the phase-type distribution. 6. Solve F(t) = P(T t) = 1 eKte;(t 0) to obtain the CDF of the sojourn time distribution. 49 To analyze a system using continuous time Markov chains, we need to approximate a general service time distribution as a phase-type distribution. We generate the phase-type representation of service time ( ;S) using the method in Section 3.2.2. The number of phases and the transition rate are determined by the SCV and the number of servers. If a server is idle, an arriving order enters service immediately, and we do not need to consider waiting time to calculate the sojourn time|the sojourn time is equal to the service time. If there is no idle server and the arriving order nds k orders ahead in the queue, its sojourn time consists of three times | the waiting time for the rst order to depart, the waiting time for the next k orders to depart, and the service time. In this condition, if one of the servers nishes its order, the order can go one step forward and then all servers are working immediately. So the order waits k + 1 \sub-waiting times" to enter service and departs the system after receiving service. Hereafter, we refer to a \sub-waiting time" as an epoch. 3.4.1 Approximating the rst epoch Each epoch is the time an arriving order spends in queue until one of the servers nishes its order. That is, an epoch starts when all servers are busy (all-busy) and ends when one of the servers nishes its order (partial-busy). To estimate the distribution of waiting time, we introduce a continuous time Markov process fN(t);t 0g with some absorbing states to model the all-busy period, where the system-state of the Markov chain is the number of servers in each phase. Asmussen and O?Cinneide (1998) showed that the waiting time distribution in a GI/PH/c queue is phase- type, and that the number of phases is 0 B@ m+c 1 c 1 CA; where c is the number of servers and m is the number of phases of each server. 50 1 2 3 Server 1: Server 2: System-state: (2, 0) Server-state (a) Server 1: Server 2: System-state: (1, 0)(b) Absorbing state Figure 3.4: System-state and server-state. Suppose there are c homogeneous servers, each with m phases, and an order nds k orders ahead in the queue. There are m + 1 server-states, and each server can be in one of those states. A order arrives to nd k orders ahead and servers in one of m states (none are idle). We describe the system-state as an m-vector of server-states. The ith element in the vector records the number of servers in state i. Example 3.1 An order arrives to a system with 2 homogeneous servers and nds 3 orders ahead in the queue. We wish to estimate the sojourn time distribution of this order given E[Ts] = 2 and C2s = 0:8. First, we approximate the general service time distribution as a corresponding distribu- tion using the tting distribution method in Section 3.2.2. We approximate with an Erlang (2;1) because m = l 1 C2s m = 2; = mE[Ts] = 1. Figure 3.4 shows two examples of the relationship between server-state and system- state: (a) shows two servers working in the rst server-state, and this is system-state (2, 0); (b) shows server 1 working in the rst server-state and server 2 has nished its order, and this is system-state (1, 0). There are 5 system-states: (2, 0), (1, 1), (0, 2), (1, 0), (0, 1). We call (2, 0), (1, 1), (0, 2) all-busy states and (1, 0), (0, 1) partial-busy states. The time between partial-busy states is an epoch, except for the rst epoch which begins upon arrival to the system and therefore to an all-busy state (otherwise a server was empty and the order entered service immediately). 51 We are interested in the distribution of the rst epoch, which is the time until the process enters the rst partial-busy state, given the process started from an all-busy state. The in nitesimal generator W1 of the rst epoch includes state changes from the all-busy states to partial-busy states. In Example 3.1 , W1 = 2 66 66 4 2 2 0 0 2 1 0 0 2 3 77 77 5 : We also need to compute the initial probability vector 1 of the the rst epoch, for which we use the stationary distribution of the all-busy states, because an arriving order nds the system in one of these all-busy states. In our approximation model, if the system reaches a partial-busy state, the system-state is changed immediately to an all-busy state because there is another order in the queue. Thus we should consider two kinds of state changes for 1. One is state changes from the all-busy states to partial-busy states, and the other is state changes from partial-busy states to all-busy states. If we sequence the states such that the all-busy states precede the partial-busy states, the former state changes have in nitesimal generator W1, which has zero value for the lower triangle below the diagonal. The latter state changes have transition rate matrix H, which has the same size as W1 and has zero value for the upper triangle and diagonal. The transition rate matrix Q = W1 +H includes all possible state changes when an order arrives at a certain position in queue. In Example 3.1, H and Q are H = 2 66 66 4 0 0 0 1 0 0 0 2 0 3 77 77 5 ;Q = 2 66 66 4 2 2 0 1 2 1 0 2 2 3 77 77 5 : 52 Now we can compute 1 using the transition rate matrix Q. The transition probability matrix function P(t) of the continuous time Markov chain is given by P(t) = eQt = I + 1X n=1 1 n!Q ntn; where I is the identity matrix. The stationary distribution is given by j = X i iPij(t): In Example 3.1, the initial probability vector 1 of the all-busy states in the rst epoch is 1 =f jg= ( 20; 11; 02) = (0:25;0:5;0:25): 3.4.2 Approximating the remaining epochs The in nitesimal generators Wk of the following epochs are the same as the in nitesimal generator of the rst epoch, because the state changes from the all-busy to partial-busy states are the same. But the initial probability vectors k of the following epochs are not the same. Rather, they must be derived successively from the initial probability in the previous epoch, because these a ect which partial-busy state was reached. As we mentioned above, if the system reaches a partial-busy state, the system state is changed immediately to an all-busy state. This means the initial probability vectors k of the all-busy states in epoch k come directly from the stationary distribution k 1 of the partial-busy states in the former epoch k 1. Now we introduce a method to compute the stationary distribution k 1 of partial-busy states in epoch k 1. Every epoch starts from one of the all-busy states and ends in one of the partial-busy states, so if we know the stationary distribution of the all-busy states and the absorbing probability from the all-busy states to partial-busy states, then the stationary 53 probability of the partial-busy state j j = X all busy i iuij; where i is stationary probability of all-busy state i and uij is the absorbing probability from all-busy state i to partial-busy state j. From the Chapman-Kolmogorov equations, uij = X all busy h Pihuhj: Finally, we get the elements of k directly from the corresponding elements in k 1. In Example 3.1, the initial probability vector 2 of the all-busy states in the second epoch is 2 = ( 10; 01;0) = (0:375;0:625;0): Note that elements of the initial probability vector k converge to the same value as the epoch k increases. 3.4.3 Approximating the sojourn time distribution of the system The sojourn time of an order with c+k orders ahead in the system has k + 1 di erent sub-waiting times and one service time. So, the state-dependent sojourn time distribution is given by the convolution of these k + 2 distributions. We construct the initial probability vector and the in nitesimal generator of the queue- ing system based on the work of Neuts (1981). The in nitesimal generator is composed of the waiting time representation ( k;W) of each order?s waiting time in the system and the service time representation ( ;S). If there are k orders ahead in the queue, the in nitesimal 54 generator and initial probability vector are given by K = 2 66 66 66 66 66 4 W W0 2 0 0 0 0 W W0 3 0 0 ... ... ... ... ... ... 0 0 0 W W0 k+1 0 0 0 0 S 3 77 77 77 77 77 5 ; = [1;0; ;0]: The CDF and PDF are given by F(t) = P(T t) = 1 eKte;(t 0); f(t) = eKtK0 = eKt( K)e;(t 0): 3.4.4 Normal Approximation Model One drawback of our approximation model is computation time for very large systems. For example, models with up to 30 servers and fewer than, say, 20 orders ahead, generally solve within 1 minute; problems with 50 servers can take 3 minutes; problems with 100 servers can take 8.5 hours and problems with 200 servers, more than a day. To address these larger problems, we introduce a second approximation called the Normal Approximation Model (NAM). The NAM uses the same phase-type representation of waiting and processing times as the approximation we have already developed, except that the CDF and PDF are not calculated with the matrix exponential expression in Equation 3.1. Instead, we assume the nal distribution is Normally distributed, and that we can compute the required rst and second moments with Equation 3.2, which requires only a matrix power computation. 55 3.4.5 Numerical results Consider the example in Section 3.3.2, but now with general service times and C2 = 0:5. We compare our results with Whitt?s method and with simulation in Table 3.3. For the comparison, we add the mean processing time to the waiting time of Whitt?s results because he computes only mean waiting time. In the simulation model, we used the Gamma distribution for the service time because it is often an appropriate model for task completion time (Law and Kelton, 2000). Table 3.3 shows the mean sojourn time results from the approximation model, Whitt?s model and the simulation model. As expected, Whitt?s model shows good results when the number of servers is high (in this case, greater than 30). With a few exceptions, the approximation model appears to perform well over a wide range of problem instances. We should note however, that the largest problems in this table are very di cult to compute using our approximation model. For example, a problem with 100 servers and 20 orders ahead takes 8.3 hours, and 200 server problems take about 1 day. We record computation time using a 2.4 GHz Intel Core 2 Duo processer, and most of the problems except for those mentioned above are computed in a few seconds. Therefore, one might view Whitt?s method and ours as complements|ours performs well for small and medium sized systems, his for large. Figures 3.5{3.7 compare the PDF and CDF of the approximation model with the simu- lation model under di erent numbers of servers and orders ahead in the queue. We use the Anderson-Darling (A-D) test to check the agreement between two distributions. As shown in Table 3.4, all test statistics are signi cant with = 5%. 3.4.6 Testing the Normal Approximation Model Table 3.5 shows the percent di erences of percentiles between the NAM and the sim- ulation when c = 2;20 and 100. We see that di erences at the 95th percentile are greater than 5 percent for cases k < c, where k is the number of orders ahead and there are c 56 Table 3.3: A comparison of mean sojourn times for general service times. Servers 2 3 5 Orders ahead 5 10 20 5 10 20 5 10 20 Whitt 18.75 34.38 65.63 11.67 21.39 40.83 11.60 17.10 28.10 % di erence -4.95 6.51 15.65 -20.99 -7.06 5.21 5.23 6.43 8.27 Approximation 19.38 31.88 56.87 14.71 23.20 40.17 10.97 16.24 26.82 % di erence -1.79 -1.24 0.21 -0.39 0.79 3.51 -0.53 1.06 3.34 Simulation 19.73 32.27 56.75 14.77 23.01 38.81 11.02 16.07 25.95 Servers 10 20 30 Orders ahead 5 10 20 5 10 20 5 10 20 Whitt 8.15 10.78 16.03 6.54 7.82 10.38 6.02 6.86 8.56 % di erence 2.25 2.81 3.42 0.81 1.20 1.49 0.68 0.93 1.32 Approximation 8.02 10.79 16.45 6.48 7.81 10.67 5.99 6.84 8.63 % di erence 0.62 2.99 6.19 -0.04 1.04 4.26 0.15 0.53 2.19 Simulation 7.96 10.48 15.43 6.46 7.70 10.26 5.97 6.79 8.43 Servers 50 100 200 Orders ahead 5 10 20 5 10 20 40 60 80 Whitt 5.61 6.11 7.12 5.30 5.55 6.06 6.03 6.53 7.03 % di erence 0.61 0.52 0.82 0.21 0.26 0.27 0.16 0.22 0.25 Approximation 5.59 6.09 7.12 5.30 5.55 6.04 6.02 6.51 7.02 % di erence 0.39 0.21 0.73 0.19 0.19 0.09 -0.01 -0.01 0.05 Simulation 5.56 5.96 7.06 5.28 5.53 6.03 5.99 6.49 6.99 20 40 60 80 100 Time 0.01 0.02 0.03 0.04 Density PDF 20 40 60 80 100 Time 0.2 0.4 0.6 0.8 1 Probability CDF Approximation Simulation Figure 3.5: Comparison of the PDF and CDF of the approximation model and the simulation model (2 servers and 20 orders ahead). Table 3.4: Anderson-Darling tests for general service time of the single-stage queue. Servers Orders ahead A-D = 5% Decision % di erence90th 95th 2 20 1.755 2.492 accept 0.92 1.78 50 10 1.828 accept 1.56 0.51 100 5 2.327 accept 0.22 0.27 57 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 Approximation Simulation Figure 3.6: Comparison of the PDF and CDF of the approximation model and the simulation model (50 servers and 10 orders ahead). 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 Approximation Simulation Figure 3.7: Comparison of the PDF and CDF of the approximation model and the simulation model (100 servers and 5 orders ahead). 58 servers, which suggests that the model does not perform well under these conditions. The explanation has two parts: First, the Normal approximation will perform well, in general, when (1) successive waiting times have the same distribution, and (2) there are many of them. Second, successive waiting times in our case are not identically distributed because of the initial set of remaining times (Section 3.4.2). However, successive waiting time distri- butions do converge as memory of the initial remaining times is \lost," which is re ected by successive initial probability vectors k converging to the same set of values. Systems with fewer servers converge faster, so, in general, the NAM performs better when k c: Table 3.5: The percent di erences of percentiles between NAM and Simulation for the single stage queue. Servers 2 20 100 Orders ahead 5 10 20 5 10 20 5 10 20 90th -1.69 0.46 1.48 -1.47 -0.53 2.53 -9.23 -1.89 -1.76 95th -2.73 0.24 2.34 -7.33 -5.84 -2.00 -14.06 -8.16 -7.74 3.5 Multi-Stage Systems 3.5.1 Approximation model We extend our work to serial lines and acyclic queueing networks based on the single stage approximation model. Consider rst a simple, 3-station serial line, and suppose an order is in queue at the rst workstation. The order will experience three sojourn times, one at each queue, and the total sojourn time is their convolution. After each of the rst two sojourn times, the status of the remaining queues changes, so we must estimate how many orders are in these queues when the order of interest arrives. After each sojourn time, the remaining queues change in two ways: orders arrive from the upstream workstation, and some orders are completed and depart. During each stage, then, we must estimate how many orders arrive to and depart from each remaining workstation (Figure 3.8). 59 5 th 4 th 2 nd Sojourn time 1 Sojourn time 2 Sojourn time 3 after E[S 1 ] after E[S 2 ] 1 1 ( ,K)? 2 2 ( ,K)? 3 3 ( ,K)? 1 1 2 2 3 3 ( ,K) ( ,K) ( ,K) ( ,N)? ? ? ?+ + = Figure 3.8: State-dependent queueing network model. Black discs represent workers; squares represent orders. For the rst workstation, we simply apply the single-stage model, giving us a sojourn time distribution and its mean E[S1]. Now, let qij be the estimated queue length of work- station j after the ith sojourn time. Waiting time at the second workstation is based on the starting queue length q02, plus arriving orders, minus departing orders. We assume that during E[S1], all q01 + c1 orders in front of the order of interest arrive to workstation 2. If all servers were busy during this time, then we would estimate the number processed at workstation 2 by c2 E[S1]=E[T2]. However, it is possible that some servers could go idle during this time, so we correct for this by using a oor function, bc2 E[S1]=E[T2]c. The expected number of orders ahead upon arrival to the second workstation becomes, q12 = max(0;q02 + (q01 +c1) bc2 E[S1]=E[T2]c): During each stage of calculation, we must revise the number in queue at remaining workstations. For the third and following workstation, this leads to qij = max 0;qi 1j + cj 1 E[Si]E[T j 1] cj E[Si]E[T j] ; for i = 1;:::;j 2: 60 1 2 3 Figure 3.9: A serial line with 3 stations. The nal revision (j 1) for workstation j is qj 1j = max 0;qj 2j +qj 2j 1 +cj 1 cj E[Sj 1]E[T j] : These approximations assume heavy tra c and are far from exact. However, they are su cient for our purposes, as we are about to show. For large serial systems, an alternative approach is to modify the NAM by \collapsing" the system into a single-stage queue. If we wish to know the state-dependent sojourn time distribution of an order in front of the rst station in a serial line with 3 stations, we assume our system is a single stage queue with c3 servers, and the number of orders ahead is q1 +c1 +q2 +c2 +q3, assuming all servers are busy. Now we can generate the state-dependent sojourn time distribution using the single stage model directly. 3.5.2 Numerical results We apply the approximation model to a serial line with 3 stations (Figure 3.9) and an acyclic queueing network with 4 stations (Figure 3.10). The black disc in each station represents an occupied server, and in the case of the acyclic queueing network, an order departing from the rst station selects its follow on station with probability p. We test state-dependent sojourn time distributions for both systems when the order is located in front of the rst station. We also test the acyclic queueing network when the order is located in front of the second station. The system information is described in Table 3.6. 61 0.67p= 0.33p= 1 2 3 4 Figure 3.10: An acyclic queueing network with 4 stations. Black discs represent workers. Table 3.6: The system information of the serial line and acyclic queueing network. Serial line E[T] C2 Servers Interarrival 1.05 0.7 Station 1 5 0.8 5 0.95 Station 2 3 0.8 3 0.95 Station 3 2 0.8 2 0.95 Acyclic queueing network E[T] C2 Servers Interarrival 1.05 0.7 Station 1 5 0.8 5 0.95 Station 2 4.47 0.8 3 0.95 Station 3 9.09 0.8 3 0.95 Station 4 2 0.8 2 0.95 62 Table 3.7: The comparison of the mean sojourn time of the serial line. Orders in queues Approximation % di erence NAM % di erence Simulation 11-12-13 47.98 0.11 46.75 -2.46 47.93 6-12-13 42.70 -0.12 41.75 -2.34 42.75 16-12-13 53.27 0.50 51.75 -2.37 53.01 11-20-13 55.13 -1.62 54.75 -2.30 56.04 17-12-17 58.33 1.15 56.75 -1.60 57.67 7-20-4 41.90 -2.28 41.75 -2.64 42.88 31-22-29 93.32 0.30 92.75 -0.31 93.04 We compare the mean sojourn time and the distribution among the approximation, the NAM, and the simulation model under di erent numbers of orders in each station. We assume service times and interarrival times in the simulation model are Gamma distributed. To estimate a sojourn time distribution in simulation for problems given in Table 3.7, we collect sojourn times only for orders seeing the required state conditions upon arrival. When an order arrives to the system, we check the number in queue at each workstation. If this vector of values corresponds to the state of interest problem condition, we record the sojourn time; otherwise we do not. This \brute force" method alleviates the problem of estimating upon arrival the remaining times of orders in process. However, because of occurrences of a particular state are rare, these simulations can take a very long time to run. For example, the rst problem in Table 3.7 (11-12-13) takes 2 days for 50 runs, with 40 million (simulated) hours in each run. Table 3.7 shows the results. In the table, the column \ Orders in queues" represents the number of orders in the respective queues. The approximation model is extremely close to the simulation results for nearly every case. The NAM is acceptable, but not quite as good. We believe the NAM underestimates the mean (notice the negative di erences) because it does not account for potential starving situations. Figures 3.11 and 3.12 compare the PDF and CDF of the approximation model with the simulation model under di erent numbers of servers and orders ahead in the queue. Again, the Anderson-Darling (A-D) tests reveal that there is no signi cant di erence between the 63 10 20 30 40 50 60 70 80 Time 0.01 0.02 0.03 0.04 0.05 0.06 0.07 Density pdf 10 20 30 40 50 60 70 80 Time 0.2 0.4 0.6 0.8 1 Probability CDF Approximation Simulation Figure 3.11: Comparison of the PDF and CDF of the serial line (11-12-13). 10 20 30 40 50 60 70 80 Time 0.01 0.02 0.03 0.04 0.05 0.06 0.07 Density pdf 10 20 30 40 50 60 70 80 Time 0.2 0.4 0.6 0.8 1 Probability CDF Approximation Simulation Figure 3.12: Comparison of the PDF and CDF of the serial line (6-12-13). simulation and approximation results. However distributions generated by the NAM do not exactly t the simulation (Table 3.9). Most of the percent di erences are greater than 5%, except when the numbers of orders at each station are 31, 22 and 29. This suggests that the NAM works well only when the total orders ahead is signi cantly greater than the number of servers at the last station, because of reduced chances of starving. Table 3.10 shows the mean sojourn time of the approximation model and simulation for an acyclic queueing network. We test two cases: the order is in front of the rst station Table 3.8: Anderson-Darling tests for the serial line. Orders in queues A-D = 5% Decision % di erence90th 95th 11-12-13 2.22 2.492 accept -3.08 -4.13 6-12-13 1.79 accept -1.61 -2.52 16-12-13 0.86 accept -2.29 -2.95 17-12-17 1.85 accept -0.84 -2.14 64 Table 3.9: The percent di erences of percentiles between the NAM and simulation for the serial line. Orders in queues Total orders ahead 90th 95th 11-12-13 44 -5.24 -7.83 6-12-13 39 -6.64 -9.24 16-12-13 49 -5.28 -7.16 17-12-17 54 -4.59 -6.82 31-22-29 90 -2.43 -3.98 Table 3.10: The comparison of the mean sojourn time of the acyclic queueing network. Order location Orders in queues Approximation % di erence Simulation The rst station 11-12-8-13 60.87 1.86% 59.77 6-4-16-13 55.59 0.76% 55.18 16-12-15-13 72.26 0.84% 71.66 11-20-12-13 73.10 0.73% 72.57 17-12-7-17 70.54 2.78% 68.64 The second station 11-12-8-13 39.56 -1.28% 40.07 6-8-12-13 33.48 -1.76% 34.08 16-12-15-13 39.56 -1.81% 40.29 11-20-12-13 51.72 -0.53% 51.99 17-12-7-17 43.56 2.25% 42.60 and in front of the second station. The approximation model appears to work well for both cases. Figures 3.13{3.16 compare the PDF and CDF of the approximation model with the simulation model under di erent locations of the order of interest. Notice the unusual distri- bution in Figure 3.14. This is caused by the mixture of the distributions of the two potential serial line paths (path 1: 1-2-4 and path 2: 1-3-4) because an arriving can take either path to depart the system. If the di erence between the two mean sojourn times is large (E[S] of the two paths are 44.44 and 75.06 hours), the sojourn time distribution has two peaks, as in Figure 3.14. Notice that our model re ects this because 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 65 20 40 60 80 100120140 Time 0.01 0.02 0.03 0.04 Density PDF 20 40 60 80 100120140 Time 0.2 0.4 0.6 0.8 1 Probability CDF Approximation Simulation Figure 3.13: Comparison of the PDF and CDF of the acyclic queueing network (The rst station, 11-12-8-13). 20 40 60 80 100120140 Time 0.01 0.02 0.03 0.04 Density PDF 20 40 60 80 100120140 Time 0.2 0.4 0.6 0.8 1 Probability CDF Approximation Simulation Figure 3.14: Comparison of the PDF and CDF of the acyclic queueing network (The rst station, 6-4-16-13). paths. We use the Anderson-Darling (A-D) test to check the agreement between the two distributions. As shown in Table 3.11, all test statistics are signi cant with = 5%. 3.6 Chapter conclusions We have developed an approximation model for state-dependent sojourn time distribu- tions of queueing systems with multiple servers, using the characteristics of the phase-type distribution. We develop an approximation model for single-stage queues with exponential and general service time distributions and an approximation model for the serial lines and acyclic queueing networks with general service time distributions. The waiting time of the single-stage G/M/c follows the Erlang distribution and has only one phase-type representation. So, it is easy to generate the CDF of the state-dependent 66 10 20 30 40 50 60 70 80 Time 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 Density pdf 10 20 30 40 50 60 70 80 Time 0.2 0.4 0.6 0.8 1 Probability CDF Approximation Simulation Figure 3.15: Comparison of the PDF and CDF of the acyclic queueing network (The second station, 11-12-8-13). 10 20 30 40 50 60 70 80 Time 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 Density pdf 10 20 30 40 50 60 70 80 Time 0.2 0.4 0.6 0.8 1 Probability CDF Approximation Simulation Figure 3.16: Comparison of the PDF and CDF of the acyclic queueing network (The second station, 6-8-12-13). Table 3.11: Anderson-Darling tests for the acyclic queueing network. Order location Orders in queues A-D = 5% Decision % di erence90th 95th The rst station 11-12-8-13 2.01 2.492 accept -2.53 -4.40 6-4-16-13 1.65 accept 1.75 -0.39 16-12-15-13 1.21 accept -0.91 -1.75 11-20-12-13 2.31 accept -3.72 -2.99 17-12-7-17 1.74 accept -1.98 -3.63 The second station 11-12-8-13 2.01 2.492 accept -3.50 -4.15 6-8-12-13 2.36 accept -4.85 -6.61 16-12-15-13 2.46 accept -3.89 -5.43 11-20-12-13 1.35 accept -2.25 -2.65 17-12-7-17 1.02 accept -0.81 -0.98 67 sojourn time distribution, whereas the G/G/c single-stage queue has many phase-type rep- resentations of waiting time. For the waiting time distribution of the G/G/c queue, we introduce a Markov process fN(t);t 0g to describe the waiting time distribution consid- ering remaining service time. The approximation model for a network was generated by convoluting sojourn time distributions in the routing of the order. In contrast with Whitt?s model, our approximation model performs well over a wide range of problem instances, whereas Whitt?s model shows good results when the number of servers is high. However our single-stage approximation model is limited to problems with fewer than about 200 servers, because the computation time of our model increases as the number of servers and orders ahead increases. This problem is caused by the matrix exponential calculation. The approximation allows us to estimate the probability that the state-dependent so- journ time will be less than a desired time t, ps = P(T < t). This result can be used in various ways based on the state of the system. For example, an order ful llment system can do dynamic order promising to customers, at the point of order. During times of heavy loading, when the probability of an order making it onto a premium shipping mode is low, the retailer might declare to o er this mode to the customer. Service systems could inform customers of their wait time with an estimated probability. For example, in an amusement park such as Disney world, we see the phrase \From this point, expected time is 1 hour" in a waiting line. Disney could develop this time at a higher probability of success and therefore \exceed expectations." Also, it could be used as a criterion for sta scheduling. If a service system has non-stationary arrival streams, we can schedule sta based on the sojourn time distribution of a peak time and non-peak time. 68 Chapter 4 Dynamic Worker Allocation Policies to Improve Service in an Order Ful llment System We propose several dynamic worker allocation policies to maximize service performance 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 distributions. Our approach di ers from other work on dynamic worker allocation in that we 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. 69 4.1 Introduction The most important function for most distribution centers is order ful llment, which is the demand-side activity that satis es customer requests. In many distribution centers, there are three primary activities: picking, packing, and shipping. Orders often arrive to such systems throughout the day and are prepared for shipping through picking and packing processes by workers. 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. Several years ago, we noticed from data of a large distributor in San Diego that on average, 5 percent of all orders arrived at the shipping dock in the 30 minutes following the departure of the last truck each night. In this chapter we ask, could workers have been redeployed somehow, perhaps in the few hours before the truck departure, to \ ush" the system of orders that were almost completed, thereby increasing the service performance of the system? For example, workers could move from picking to shipping in order to ush the system of nearly completed orders right before the truck departs. The problem, of course, is that this causes an imbalance of the system and reduces the throughput of the picking area temporarily, which has reverberating e ects later. Whether or not the system can be gamed in this way is what we explore in this chapter. Our goal is to maximize a performance metric called Next Scheduled Departure (NSD), which was introduced by Doerr and Gue (2006) to measure service performance in a distri- bution 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, whereas existing metrics such as WIP, cycle time and throughput are from the distributor?s perspective. 70 To introduce the dynamics of a worker allocation policy, consider a na ve policy, which moves a xed number of workers at a xed time close to the deadline everyday. We model the order ful llment system as a serial line of 3 workstations (picking, packing, shipping) with 10, 12, and 9 workers per workstation. We move 8 workers from picking to shipping one hour before the deadline everyday, in an attempt to ush the system of almost completed orders. On some days, the policy is e ective, but on others it is not. Because the na ve policy ignores the state of the system when moving workers, it is possible that workers are moved from a busy work center (picking) to a work center with little or no work to be done (shipping). This is a double error: (1) busy workers are made idle after the shift and (2) orders in picking back up during the shift, thereby increasing their chances of being late the next day. Simulation experiments on this simple policy 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 a system will make it on the next departing truck, if the allocation of workers does not change. When this probability drops below a threshold value, then some workers in picking should be moved so as to improve the probability of success for other, more promising orders. We introduce several dynamic worker allocation policies designed to maximize customer satisfaction without disturbing system balance. In the next section, we describe previous research on dynamic worker allocation policies. Section 4.3 presents the meaning of the Next Scheduled Departure (NSD) metric in terms of the sojourn time distribution. We also de ne the probability of success (ps), which is the crucial criterion for our worker allocation policies. In Section 4.4 and 4.5, we introduce and evaluate several policies. We o er conclusions in Section 4.6. 71 4.2 Literature Review There is a rich literature on workers moving among tasks in a dynamic fashion in man- ufacturing and, to a lesser extent, warehousing. Research in dynamic worker allocation falls under several names, such as worksharing systems, worker cross-training or cross-utilization, collaborative and non-collaborative work systems, and agile workforces. The use of cross- trained workers to improve the performance of manufacturing systems has been investigated by several researchers. Askin and Chen (2006) categorize worksharing policies as two types: Dynamic assembly-Line Balancing (DLB) and Moving Worker Modules (MWM). MWM systems have fewer workers than machines, so workers share every machine and move ev- erywhere within a zone, whereas DLB has the same number of machines and workers. DLB has xed tasks and shared tasks; the former are assigned to a designated worker, and the latter can be carried out by either of an adjacent pair of workers. In this review, we classify the oating worker system studies by MWM or DLB based on the number of workers per machine, the extent of skill per worker, and the amount of WIP. Toyota Sewn-products management System (TSS) and the Bucket Brigade System are the representative MWM (Askin and Chen, 2006). These systems have fewer workers than machines. In a TSS system, a worker moves downstream doing his or her job at every machine within a restricted zone until either the job is nished or the downstream worker takes over the job. Bischak (1996) studied a U-shaped manufacturing line under the TSS system. She developed a simulation model and showed that the throughput of moving worker modules is higher than the xed worker system with and without the use of bu ers. Zavadlav et al. (1996) also examined a U-shaped serial line under the TSS system using a Markov decision process and simulation. They compared three systems: xed assignments, overlapping (shared) worker assignments, and free- oating worker assignments based on the same constraints. They showed that the free- oating system is the most e cient way to lower WIP. 72 Bartholdi and Eisenstein (1996) and Bartholdi et al. (2001) developed the bucket brigades protocol based on TSS. They showed that sequencing workers from slowest to fastest can make a balanced partition of work emerge without management intervention, and that it produces the highest possible production rate. McClain et al. (2000) showed that bucket brigades does not work well with random processing times and nearly-equal worker velocities, because these conditions can lead to idle time for workers when the bucket brigade rule is followed strictly. To increase system performance, they suggested relaxing bucket brigade?s \wait" restriction and allowing small amounts of inventory. In the DLB problem, Gel et al. (2007) developed a zone policy in CONWIP produc- tion systems with hierarchical cross-training. They suggest the \ xed-before-shared" policy, which states that a cross-trained worker should do his unique task before helping a xed worker with a shared task. They showed that partial cross-training in hierarchical patterns may yield substantial bene t if cross-trained workers are signi cantly faster than static workers under the \ xed-before-shared" policy. Gel et al. (2002) studied factors a ecting the opportunity for worksharing under the same conditions. They showed that the ability to preempt the shared task, granularity of the shared task, and less variability of the task times increases the system e ectiveness. The half full bu er (HFB) control policy, which provides enough work for the down- stream worker not to starve and enough empty space for the upstream worker not to be blocked, was studied by Ostalaza et al. (1990), McClain et al. (2000), Gel et al. (2002), Chen and Askin (2006) and Askin and Chen (2006). Ostalaza et al. (1990) suggested an SPT/RSPT rule, and showed that the system performance increases when the number of shared tasks is controlled by the HFB control policy. McClain et al. (2000) considered zoned worksharing and bucket brigades with a variety of situations. Using simulation, they tested the HFB control policy to move exible workers for shared tasks. Following this research, Gel et al. (2002) proposed a more delicated HFB control policy considering general processing times and di erent worker speeds. They also suggested the 50-50 work content heuristic rule 73 that an upstream machine transfers the shared task to the downstream machine in the case where the ratio of work content at the upstream machine to the total work content in the system is greater than 50 percent. Askin and Chen (2006) and Chen and Askin (2006) con- rmed again the e ect of the HFB control policy and suggested another threshold heuristic rule, smallest R no starvation (SRNS), to calculate the threshold value R under CONWIP levels. They showed that cross-training can increase system e ectiveness with a work content heuristic rule such as SRNS and HFB. In contrast to existing research, we o er a dynamic worker allocation policy with the following distinctives: First, the purpose of most former research on worker allocation has been to maximize system throughput, to minimize cycle time, or to decrease WIP in a manufacturing system. Our objective is to improve customer service in an order ful llment system, which we measure with the Next Scheduled Departure (NSD) performance metric (Doerr and Gue, 2006). Second, most of the existing research is limited to small tandem lines, such as 2 or 3-station tandem lines with one or two servers in the system. Existing results show that dynamic capacity allocation between adjacent workstations works well based on those simple conditions. However, real systems may not be so simple: there may be many servers in each workstation and perhaps many workers need to move among workstations to improve system performance. For example, a typical distribution center can have dozens of workers in each functional area (picking, packing, shipping). Our policies are e ective in these larger systems. Third, most of the methods which have been used by existing research are Markov decision processes, heuristics, or simulation. We compute the probability that an order in the system will make it on the next departing truck without worker allocation based on its sojourn time distribution in a multi-server queueing network. When this probability drops below a threshold value, we move some workers from one workstation to another to improve the probability of success for more promising orders. 74 Cutoff t c 0 Cutoff t c 1 Deadline t d 0 Deadline t d 1 Day 0 Day 1 Arriving orders in this period Are due by this deadline Figure 4.1: The relationship between the cuto time and NSD (Doerr and Gue, 2006). 4.3 Preliminaries Before designing worker our allocation policies, we introduce three fundamental concepts integral to their development: Next Scheduled Departure (NSD), state-dependent sojourn time distributions, and the probability of success. 4.3.1 Next Scheduled Departure (NSD) Doerr and Gue (2006) de ne NSD as the fraction of orders arriving in 24 hours between cuto times tc on successive days that make it onto the truck in the current day. An earlier cuto time means more orders arriving before the cuto time will make it on the departing truck and therefore NSD will be higher; thus, NSD depends directly on the cuto time. Figure 4.1 illustrates the relationship between the cuto time and NSD, where tnc is the cuto time and tnd is the deadline of day n. Expected NSD for a system can be computed based on the steady-state sojourn time distribution. Every arriving order sees the same steady-state sojourn time distribution to depart the system. If an order arrives to a system close to the deadline, the probability that it will make it on the next departing truck will be low. Those arriving early in the 24 hour window have the highest chance of making it on the truck. As shown in Figure 4.2, given cuto time tc and deadline td, the time average of the probability of being completed on time for orders arriving in this 24 hour period is the expected NSD. We can compute this 75 5 10 15 20 25 30 Time 0.2 0.4 0.6 0.8 1 Probability CDF 24 hours ? ?+24 Figure 4.2: The meaning of NSD based on the steady-state sojourn time distribution. steady-state sojourn time distribution for multi-server queueing networks using methods we develop in Gue and Kim (2009a). Formally, if P(t) = P[ event occurs at time t ], then the time average probability is given by p = 1t Rt0 P(t)dt. Proposition 4.1 For an order ful llment system with deadline td and cuto time tc, NSD = 124 Z +24 P[T < >: P[T te] if an order is in service (4.2) where tr the remaining is the remaining time, t? is the duration of time between entry into service and the deadline, and te is the elapsed time since an order entered service. We can compute ps for an order in service from the sojourn time CDF as follows. P[T te] = 1 P[T >t?jT >te] = 1 R1 t? fT(u)duR 1 te fT(u)du = R1 te fT(u)du R1 t? fT(u)duR 1 te fT(u)du = Rt? te fT(u)duR 1 te fT(u)du = F(t?) F(te)1 F(t e) : To illustrate, we record ps for a particular order in a simulated system. Example 4.1 Consider a single-stage queue with 30 identical servers, with processing times having mean 5 hours and squared coe cient of variation 0.8. An order arrives at 10:00 and sees 19 orders ahead. The last truck departs the system at 17:00. What is the probability of success, and how does it change with time? 79 To illustrate how ps changes over time, we run a simulation model and compute ps and the mean sojourn time E[S] at 15 minute increments, based on a recalculated state- dependent sojourn time distribution. For example, to compute the probability of success of this order at 10:00 (remaining time tr=7 hours), we generate the phase-type representation of the system ( ;K) using the approximation model of Gue and Kim (2009b) based on the system state. Then the probability of success is computed using Equation 4.2, ps = P(T tr = 7) = 0:41: Figure 4.3 shows ps for the order over time. In this case, the order is still in the system at td=17:00, when its ps=0. We see that ps increases or decreases according to the condition of the system. However, once an order enters service, its ps can only decrease (see Equation 4.2). Once an order is seized by a worker, there is no way to help its process time by worker allocation because we do not consider worker collaboration on an order. Our worker allocation policies are based on the insight that as time progresses, the probability of success for orders far from completion goes down and so e ort devoted to them might be better spent on orders \on the bubble" of making the next truck. We formalize this idea below. 4.4 Worker allocation policies We introduce dynamic worker allocation policies that assign cross-trained workers to tasks in an order ful llment system such that performance against the Next Schedule De- parture (NSD) metric is improved. We suggest two types of policies: a ushing policy and a cascade policy. We further classify a ushing policy as either single or multi- ush in accordance with the number of switching events. A single ush policy moves workers based on the state of the shipping area at a de ned time in order to ush the system of nearly completed orders right before the truck departs 80 11:00 12:00 13:00 14:00 15:00 16:00 17:00 0.1 0.2 0.3 0.4 p s Time In queue In service Figure 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. (Figure 4.4). For this reason, this policy only focuses on the state of the shipping area and the sojourn time distribution of orders in it. One limitation of the single ush policy is that we have only one opportunity per day to correct an overloaded shipping area. In the multi- ush policy, we move workers periodically during a day, considering the state of the shipping area. Multi- ush is the same as single ush, except for the number of switching events. A cascade policy moves workers sequentially from picking to packing, then from packing to shipping based on the state of the packing and shipping areas at a certain times. The cascade policy is an attempt to ush across more than one workstation. Figure 4.5 shows the concept. At a xed time, say, 15:00 each day, move workers from picking to packing. At a later time, after ushing orders in the queue of the packing area, move workers from packing to shipping. The goal of the cascade policy is to increase the ps of an order in the system regardless of its location in a relatively short time. 81 9 th 4 th p s =5% p s =60% t d =17:00 t s =16:00 + 16:00 17:00 Figure 4.4: Flushing policy 5 th 2 nd 5 th 1 st Phase I Phase II t d =17:00 t s 1 =15:00 t s 2 =16:00 t d =17:00 15:00 + 16:00 + Figure 4.5: Cascade policy 82 Table 4.1: The characteristics of 12 systems. High High Low Low High SCV Low SCV High SCV Low SCV Small system with short E[S] System 1 System 2 System 3 System 4 Large system with short E[S] System 5 System 6 System 7 System 8 Small system with long E[S] System 9 System 10 System 11 System 12 4.5 Experiments To make de nitive statements about which policy is \best" would require that we test it on an in nite number of potential systems. We settle for a more limited experiment on 12 carefully constructed systems, which we believe stress the model along important dimensions. All of the example systems are 3-stage serial lines, in keeping with the picking- packing-shipping structure of many order ful llment systems. The systems vary according to size (number of workers), mean sojourn time, variability of processing times, and utilization. Table 4.1 summarizes the systems. We select three basic candidate systems according to size and mean sojourn time, then divide each system into four di erent systems based on two levels of the squared coe cient of variation (SCV: 0.5, 0.9) and utilization ( =0.85, 0.95). For the squared coe cient of variation, we deal exclusively with the case SCV<1 because it is common in practice. We consider three basic serial lines. The rst is a small system (total 31 workers) with short mean sojourn time (6.52 hours); the second is a large system (total 126 workers) with short mean sojourn time (7.31 hours); and the third is a small system (total 44 workers) with long mean sojourn time (19.76 hours). Table 4.2 shows system information of the representative systems 1, 5 and 9. To implement a worker allocation policy correctly, we need to answer the questions: how many workers should we move and when? Switching workers too early and switching 83 Table 4.2: System information for Systems 1, 5, and 9. System 1 E[T] C2 Workers Interarrival 0.117 0.75 Picking 1.07 0.9 0.91 10 Packing 1.3 0.9 0.93 12 Shipping 1.0 0.9 0.95 9 System 5 E[T] C2 Workers Interarrival 0.053 0.75 Picking 2.7 0.9 0.92 56 Packing 1.5 0.9 0.95 30 Shipping 2.0 0.9 0.95 40 System 9 E[T] C2 Workers Interarrival 0.23 0.75 Picking 3.3 0.9 0.96 15 Packing 4.3 0.9 0.93 20 Shipping 2.0 0.9 0.97 9 too many workers can cause imbalance in the system, whereas switching too late or too few might not have an e ect on NSD. The exact number of switching workers is decided by the target ps (the probability of success of an order of interest that we want to reach by worker allocation), and switching time ts, so we test the 12 systems to nd best ts and the target ps for each policy. For example, in Figure 4.4, the probability of success ps for the last order in the shipping queue is 5% based on the number of orders ahead 8, the number of workers 6, and the remaining time tr=1 hour (td ts). To reach the target ps=60% for this order, we add 5 workers to shipping, then the system state is changed: the number of orders ahead is 4, and the number of workers 11. The number of workers required is computed by simple search, based on the method of Gue and Kim (2009b) Target ps = 60% = P(T tr = 1) = 1 eKtre: where ( ;K) depends on the system state. 84 We use 6 di erent levels of the target ps; 10%, 30%, 50%, 60%, 70%, 90%, and use switching time ts from early in the morning to close to the deadline. We do not include a switching time that makes the remaining time tr less than the mean processing time of the last workstation because this situation limits the e ect of the policy (more on this below). We run simulation models built in Arena 12.0. Each scenario runs 100 times and each run has 50 days? NSD results. For easy implementation of the simulation model, in advance, we compute the exact number of workers to reach the given target ps based on current ps for possible combinations between target ps and switching time ts for each system. Then, we just move required workers based on the number in queue in the shipping area based on this given information. For example, in the case with the target ps=60% and switching time ts=16:00 for System 1, suppose there are 4 orders in the shipping queue at 16:00, then we just move 4 workers from picking to shipping based on the given information. In our simulation model, we assume two important points: We do not consider worker transition time when we move workers from one workstation to other workstations. This is not true in practice, of course, but transition times are highly application speci c, so we ignore them for this analysis. We also assume that the picking process is not a batch process. For many order ful llment systems, this will not be true, but for others (part-to-picker systems, in particular) it is reasonable. 4.5.1 Single Flush The truck departure time for all experiments is 17:00. To nd the best ts and target ps for the single ush policy, we use Target ps: f10%, 30%, 50%, 60%, 70%, 90%g, and Switching time ts: f16:00, 15:00, 14:00, 13:00, 07:00g for Systems 1{4, f15:00, 14:00, 13:00, 12:00, 07:00g for Systems 5{8, f15:00, 14:00, 13:00, 12:00, 07:00, 02:00g for Systems 9{12. 85 We choose candidate ts from early in the morning (02:00 or 07:00) to close to the deadline (16:00 or 15:00). We do not consider ts later than 16:00 (Systems 1{4) or 15:00 (Systems 5{12) because mean service times of the shipping area E[T3] are 1 or 2 hours. If the remaining time tr is shorter than E[T3], then ps of the switching time is too low to have an e ect on NSD. The algorithm is: Algorithm 1 Single Flush (SF) Step 1: Decide the cuto time tc using the steady-state sojourn time distribution, based on the baseline NSD. Step 2: Check the number of orders in front of the shipping area at a given switching time ts. Step 3: Compute ps of the last order in the shipping queue using the state-dependent sojourn time distribution. Step 4: If ps is below a given target, compute the number of switching workers required to reach the target ps and move these workers from the picking to the shipping area. Otherwise, do not move workers. Step 5: Return workers to their original place when clock time reaches the deadline td. First, we decide cuto time tc using the steady-state sojourn time distribution model under baseline NSD 77% ( xed model without worker allocation of System 1): Baseline NSD = 77% = 124 Z 24+ P[T