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