Wave Release Strategies for Order Ful llment Systems with Deadlines by Erdem C even 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 August 3, 2013 Keywords: Order release control, wave processing, time guaranteed delivery, cuto time, uid approximations Copyright 2013 by Erdem C even Approved by Kevin R. Gue, Tim Cook Associate Professor of Industrial and Systems Engineering Jorge Valenzuela, Professor and Chair of Industrial and Systems Engineering Chase Murray, Assistant Professor of Industrial and Systems Engineering Abstract A distribution center is ultimately in the business of shipping orders. Because every customer?s experience is determined by how timely and accurately an order is received, distribution centers should focus not only on how fast they process orders, but also on keeping their operations synchronized and delivering orders on time. In traditional distribution centers, including internet and catalog ful llment systems, synchronization and economies of scale are achieved with \waves." A wave is a large group of orders that are picked and packed together, sorted for individ- ual customers, and shipped. Problems arise when those internal processes are not properly coordinated. When operational parameters, such as number and timing of waves, are not planned systematically, missed orders and upset customers are the result. We present the rst systematic investigation of wave planning in order ful llment sys- tems. We develop continuous uid models of work content in an order ful llment system, in order to determine the optimal timing and number of waves with the goal of maximizing service performance against a deadline-oriented metric. ii Acknowledgments I would like to express my deepest gratitude to my advisor, Dr. Kevin R. Gue, for his excellent guidance, patience, brilliant ideas and providing me with an excellent atmosphere for doing research. I am most grateful for his help and support. I would also like to thank my committee members, Dr. Chase Murray and Dr. Jorge Valenzuela for their thoughtful criticism, time and suggestions. I am grateful to the Naval Postgraduate School, which partially supported this research under Grant NPS-BAA-11-002 at Acquisition Research Program, Auburn University Industrial and Systems Engineering Department and Graduate School, which also funded my graduate work. I would also like to thank my parents, H useyin and S uheyla and my sister Z umr ut. Without their presence, endless support, great understanding and love, my research would not have been possible. They always support me and encourage me with their best wishes. I would like to thank my all friends who always supported me: Eren Sak n c, Kaan Berbero glu, Polly Chen. They were always there cheering me up and stood by me through the good times and bad. This dissertation is dedicated to my parents. iii Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Order Flow in Distribution Centers . . . . . . . . . . . . . . . . . . . . . . . 2 1.2 Wave Processing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.3 Next Scheduled Deadline Metric . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.4 Literature Review . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.5 Problem Statement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 1.6 Organization of the Dissertation . . . . . . . . . . . . . . . . . . . . . . . . . 12 2 Optimal Release Strategies for Order Ful llment Systems with Deadlines . . . . 14 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.2 Single Wave Systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.3 Multiple Wave Systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.4 Optimal Number of Waves . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 2.5 Validation of Fluid Approximation . . . . . . . . . . . . . . . . . . . . . . . 25 3 Setting Wave Release Times in the Presence of Uncertainty . . . . . . . . . . . . 29 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 3.2 The Case of a Single Wave . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 3.3 Adjusting Multiple Wave Release Times . . . . . . . . . . . . . . . . . . . . 31 3.4 An Empirical Procedure to Adjust Release Times . . . . . . . . . . . . . . . 34 3.5 Setting the Cuto Time . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 iv 3.6 Adjusting The Number of Waves . . . . . . . . . . . . . . . . . . . . . . . . 36 4 Wave Release Strategies for Systems with Multiple Order Classes . . . . . . . . 43 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 4.2 Single Wave Systems With Multiple Order Classes . . . . . . . . . . . . . . . 45 4.3 Multiple-Class, Multiple-Wave Systems . . . . . . . . . . . . . . . . . . . . . 52 4.3.1 A Two-Class, Two-Wave Order Release System . . . . . . . . . . . . 53 4.3.2 Systems with More Classes and Waves . . . . . . . . . . . . . . . . . 63 4.4 Implications for Practice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 5 Conclusions & Future Research . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 Appendices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 A Feasibility and Stability Conditions . . . . . . . . . . . . . . . . . . . . . . . . . 85 B Comparison of Search Methods and Details of Di erential Evolution . . . . . . . 87 v List of Figures 1.1 Orders arriving between consecutive cuto times are due on the next deadline. . 5 2.1 A single class, three wave system. di indicates the deadline on day i; wj is the release time for wave j. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.2 Di erent realizations of optimal w1. . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.3 Optimal wave release times: = 0:5: . . . . . . . . . . . . . . . . . . . . . . . . 21 2.4 Optimal wave release times: = 0:75: . . . . . . . . . . . . . . . . . . . . . . . 21 2.5 Optimal wave release times: = 0:95: . . . . . . . . . . . . . . . . . . . . . . . 21 2.6 Two perspectives on maximum possible NSD for di erent numbers of waves. . . 22 2.7 Change in NSD with respect to the number of waves when there is a xed time component associated with each release. . . . . . . . . . . . . . . . . . . . . . . 25 2.8 Accuracy of the uid approximation improves as orders become more indivisible. 26 3.1 Service performance versus utilization , for three wave release times. . . . . . . 30 3.2 Each curve shows how NSD changes for di erent levels of observed utilization with respect to planned utilization. . . . . . . . . . . . . . . . . . . . . . . . . . 33 3.3 Maximum expected NSD achieved when release times are adjusted to earlier times. 33 3.4 Ten di erent levels of planned utilization for the case study (left) and correspond- ing average NSD (right). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 vi 3.5 Performance results when the cuto time is adjusted. On the left, NSD for each of the 80 days in the data, sorted from high to low. On the right, Type 1 and Type 2 performance for di erent levels of cuto time setback. . . . . . . . . . . 36 3.6 Expected NSD considering only uncertainty (bottom curve) and expected NSD considering both uncertainty and the xed time component (top curve). . . . . 38 3.7 Loss associated with using N waves. . . . . . . . . . . . . . . . . . . . . . . . . 41 4.1 A three-class, three wave system. . . . . . . . . . . . . . . . . . . . . . . . . . . 43 4.2 Case 1: w1 d2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 4.5 Case 4: w1 >d1 and f1 d1 and f1 >d2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 4.7 Change in the system NSD in a two class, single wave system. . . . . . . . . . . 50 4.8 A wave is comprised of orders from their previous cycles (unworked inventory) and orders that arrive in the current day. . . . . . . . . . . . . . . . . . . . . . . 54 4.9 Solutions for a two class, two wave system ( 1 = 2). . . . . . . . . . . . . . . 58 4.10 Timing of waves and class mixtures when = 0:50. . . . . . . . . . . . . . . . . 60 4.11 Timing of waves and class mixtures when = 0:75. . . . . . . . . . . . . . . . . 60 4.12 Timing of waves and class mixtures when = 0:95. . . . . . . . . . . . . . . . . 60 vii 4.13 Admissible pairs of (w1;w2) indicate that the objective function surface is non- convex. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 4.14 Analytical models assume the server randomly selects orders. . . . . . . . . . . 69 4.15 Di erent sequencing rules result in di erent simulated average system NSD. . . 70 4.16 Model validation results with simulation when = 0:5. . . . . . . . . . . . . . . 71 4.17 Numerical solutions of analytical models. . . . . . . . . . . . . . . . . . . . . . . 73 B.1 Di erential Evolution Pseudo-code . . . . . . . . . . . . . . . . . . . . . . . . . 89 viii List of Tables 1.1 Operating characteristics of DCs (Source: van der Berg, 2010). . . . . . . . . . . 3 1.2 The top ten largest U.S. internet retailers? shipping policies (Source: http:// www.internetretailer.com/). . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.1 Validation of uid model via simulation. . . . . . . . . . . . . . . . . . . . . . 28 3.1 Adjusting release times in the presence of workload uncertainty improves ex- pected NSD. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 4.1 Cuto times for di erent deadlines d1 in a two class, two wave system ( 1 = 2). 62 4.2 Summary of performance metrics. . . . . . . . . . . . . . . . . . . . . . . . . . . 72 4.3 Computational times in seconds. . . . . . . . . . . . . . . . . . . . . . . . . . . 74 B.1 Comparison of objective function values in the preliminary run. . . . . . . . . . 88 ix Chapter 1 Introduction Within a supply chain, products need to be physically moved from the point of origin to the point of consumption. During this process, order ful llment centers receive product from suppliers, store it for a certain period of time, and ful ll customer orders. In this context, many production and distribution systems can be de ned as order ful llment systems. Distribution centers (DCs) receive orders from geographically dispersed customers who have increased expectations of getting their products at a desired level of quality and promised time of delivery. With the growing success of e-commerce, distribution centers often receive a large number of small orders which have to be ful lled within very tight time windows. From 1993 to 2007, smaller size shipments increased by around 107% . The result of these trends has been more complex distribution centers and more tightly controlled order ful llment systems. Today, the world?s newest and most successful companies, such as Amazon, Zappos, and WalMart have increased competitiveness by o ering better service promises. For example, Amazon customers are familiar with this deadline-driven o er: \Want it delivered Monday, July 11, 2013? Order it in the next 3 hours and 42 minutes, and choose One-Day Shipping at checkout." Amazon is taking its service one step further by investing billions to make next-day delivery standard and same-day delivery an option (Manjoo, 2012). But, how should DCs align their functions properly in order to ful ll orders in these deadline oriented environments? 2007 Commodity Flow Survey. 1 In our research, we address the problem of controlling order ow within a DC which operates in a deadline oriented environment. Our main purpose is to develop optimal ow control policies for DCs to increase their service levels. 1.1 Order Flow in Distribution Centers The basic ow of an order in a distribution center starts with receiving and put away. Receiving is the activity of gathering products from suppliers. The put away process moves products from receiving to designated storage locations. Because the ow is from suppliers to the DC, these activities are usually called inbound operations. Once the products are located in storage locations, they are available for outbound activities. These activities mainly include order picking, accumulation, packing, and shipping. Order picking, which involves gathering customer orders from storage locations, is the major activity in most DCs. When multiple orders are picked in groups (batches), they must be sorted based on destination. Because multiple customers may be assigned to the same destination, an accumulation process is necessary. Accumulated orders are transferred to packing and then to the shipping docks. Among all ful llment activities, order picking constitutes a signi cant portion of the total warehouse operating expense (de Koster et al., 2006). Although automated picking systems are also common (e.g. automated storage and retrieval systems, robotic picking), the majority of DCs employ picker-to-parts systems (de Koster, 2004). In picker-to-parts systems, pickers either pick orders individually (discrete picking) or in batches. The term sort-while-pick is used when pickers pick multiple products simultaneously and immediately sort them. The term pick-and-sort is used if sortation is after picking. From an operational point of view, pickers are segregated into zones and in each zone, pickers may either pick orders in a synchronized (pickers in di erent zones pick the same batch of orders simultaneously) or progressive way (where an order is completed after sequential processing in all zones). 2 1.2 Wave Processing Within a DC, picking orders in di erent ways has long been standard practice. One-size and multi-size batch picking are common in practice, especially when the number of orders (workload) is small. As daily workload increases, DCs require more complex processes. Daily workload can also uctuate signi cantly in complex ful llment environments. In highly complex systems, for example, maximum workload can even be more than twice the average daily workload (van der Berg, 2010). As a consequence, operating characteristics of DCs change dramatically (Table 1.1). Table 1.1: Operating characteristics of DCs (Source: van der Berg, 2010). Complexity Low Medium High Number of operators 15 15-45 45 Warehouse area (m2) 5;000 5,000-15,000 15;000 Orders per day 1;000 1,000-1,000 5;000 Warehouse area 5;000 5,000-15,000 15;000 Separate picking areas 1 2-3 4 Shipping pattern 1 wave 2-3 waves Individual Fluctuating workload 50% 150% 250% Orders can be picked in a single zone, representing the entire picking region, or in multiple zones. DCs with medium or high complexity require multiple zones. When multiple pickers in di erent zones pick for the same pool of orders simultaneously, order release is controlled with waves. A wave is basically a set of orders grouped by some criteria and which is released to the oor for processing at the same time. Attributes for grouping might be a mode of transportation, a group of stores in retail, high priority orders, orders requiring a speci c type of value-added service, or even an individual customer if it orders in large enough quantities. A basic wave planning and control process involves creation, release, and monitoring. In practice, wave planning is performed manually. Wave planners decide the criteria that determine which orders to include in the wave. Creating waves of orders that are all shipping via UPS, for instance, can ensure that those orders are ready when the trailers are scheduled 3 for pick-up (Franco, 2006). Wave planners are also responsible for determining the size of waves and releasing them to the oor. If the size (or length) of the wave is too large, it combines tasks and increases pickers productivity; however, it may not be possible to complete by the scheduled truck departure. On the other hand, if the workload is too small, productivity of the pickers will be low, but more orders will be ready for shipment by their deadline. Most of the world?s leading warehouse management systems (WMS) providers, such as SAP and Oracle, o er automatic wave management systems to alleviate ine ciencies of human-based wave management (SAP Business Solutions, 2012; Oracle Warehouse Manage- ment User?s Guide, 2012). Although these systems are able to monitor the progress of waves dynamically, they lack the ability to prioritize orders when the wave length is longer than two hours. As a consequence, wave planners usually release small waves manually, which creates underutilized resources and does not guarantee completion time of a wave by its deadline. Therefore, there is a need to develop policies that determine the number and size of the waves, their release times, and their contents. In practice, many DCs release waves sequentially. That is, a new wave does not be- gin until all the orders in the current wave are completed. Sequential processing of non- overlapping waves is sometimes referred as xed wave systems (Bozer et al., 1988). When orders are grouped in waves based on di erent attributes, xed waves create an easily con- trolled ow of work. For example, when destinations are sorted by distance from the DC and placed into waves (more distant customers are in the rst wave, to allow more time for transportation; nearer customers are assigned to later waves), there is no need for additional sortation. Because xed waves require pickers (and packers) to wait until all others complete their picks, it creates idleness. To prevent the potential idleness of workers, waves can over- lap (called dynamic waves). Once the orders are dynamically sent downstream, they have to be sorted out. A sortation system usually requires high investment, but more importantly, it adds additional processing time. 4 Cuto Deadline Cuto Deadline Orders arriving in this cycle . . . are due by this deadline Day 0 Day 1 Figure 1.1: Orders arriving between consecutive cuto times are due on the next deadline. 1.3 Next Scheduled Deadline Metric Whether a DC employs discrete, batch, or wave picking, ful llment processes should be designed in such a way that customers receive their orders at the promised time. When balancing internal operations and shipping deadlines, continuous processing versus a cyclical transportation environment should be considered. Doerr and Gue (2011) show that it is more important to make orders ready when the truck is on the shipping dock than it is to minimize ow time. The authors propose a metric called Next Scheduled Deadline (NSD), which measures the fraction of orders arriving during a speci c 24-hour period that are processed before a speci c truck departure (Figure 1.1). The 24-hour time window is de ned by consecutive \cuto times." We de ne cycle to specify the time interval between successive cuto times. To manage customer expectations, DCs publish \cuto times" and promise to ship orders by a deadline when they are ordered before that time. If a customer places his order after the cuto time (i.e., after the current cycle), he has no expectation that his order will ship that day (the order will be scheduled to the next day?s deadline). If the order is placed before the cuto time (i.e., before the current cycle), the customer expects his order to be shipped the same day. The cuto time clari es expectations and provides a measurable, customer-focused goal for DCs. When the cuto time is published, the ful llment system?s objective is to ship all orders that arrive between two consecutive cuto times by the deadline (which is usually before a speci c truck departure). While an early cuto time implies loss in revenues, a late cuto time 5 means too many disappointed customers. Consequently, the cuto times should be set as late as possible, while at the same time, minimizing compensations due to late deliveries. Table 1.2 shows the cuto times and the fastest delivery promises by the ten largest internet retailers (Internet Retailer, 2013). Table 1.2: The top ten largest U.S. internet retailers? shipping policies (Source: http://www.internetretailer.com/). Companies Shipping frequency Cuto time Compensation 1. Amazon.comInc. Overnight 4:30p.m. EST- Prime* Refund 6:45p.m. EST- Non-prime Refund 2. StaplesInc. Overnight 5 p.m. EST N/A 3. AppleInc. Overnight 3 p.m. local N/A 4. Walmart.com Nextbusiness day N/A N/A 5. DellInc. Nextbusinessday 5 p.m. local N/A 6. O ceDepot Inc. Sameday 10 a.m. local N/A 7. QVC (Liberty Corp.) Overnight 12 p.m. EST N/A 8. SearsHoldingCorp. Overnight 7 p.m. EST N/A 9. Net ixInc. Nextbusinessday N/A N/A 10. CDW Corp. Nextbusinessmorning 5 p.m. EST N/A * To increasecustomerloyalty, Amazon o ers two-day free shippingto its primecustomersand prioritizetheir orders. In addition to top ten internet retailers, many other companies publish their cuto times to attract customers and gain a competitive advantage. Amazon has taken its service one step further by o ering same day delivery with full refunds for late deliveries. In its recently introduced local express delivery policy, Amazon states \When you order using Local Express Delivery before certain cuto times, your items will be delivered the same day. If the order is placed after the deadline, your order will be delivered the next business day. For same day delivery, you can usually order as late as the times below" (Amazon.com, 2013). Because NSD counts the number of orders that arrive between two consecutive cuto times and how many of them are actually shipped, a 100% service performance can be achieved with an optimal wave release policy which maximizes the NSD and properly determines the cuto times. 6 1.4 Literature Review The relevant literature can be classi ed under three main categories: (1) Research which focuses on gaining competitive advantage by increasing customers? experiences, (2) Studies addressing the design and control of order ful llment systems, and (3) Studies relevant to uid approximations to queueing networks. The rst group of studies addresses time-based competition of ful llment systems in which the quality of service (Shang and Liu, 2011) is driven by how timely and accurately their orders are ful lled. These works, in general, address delivery time quoting, pricing, and capacity decisions with the goal of maximizing expected pro t. Most of the studies in this group focus on a uniform delivery time|delivery within a certain time windowy, for all customers. Because demand is driven by the lead time and the price, those studies focus on nding pro t maximizing strategies in time competitive environments. The ful llment system (generally a make-to-order company) selects a uniform delivery guarantee. When the published delivery lead time is short, more customers are attracted which may result in late delivery penalties. To reduce the risk of late completions, companies may publish longer quotes which may cause customers to seek other competitors. Chatterjee et al. (2002) present a two stage decision model to maximize the pro t for a company that quotes uniform delivery dates with incomplete information on real processing time (which makes it hard to estimate future capacity levels). In the rst decision step, the capacity is not known and the delivery time is quoted by the company (e.g. the marketing department) which a ects the operating characteristics (i.e. the delays in the system) in the second stage. The model considers an environment in which processing times and arrival rates are stochastic and consequently the decisions should re ect the unpredictability of new order arrivals and capacity. So and Song (1998) discuss the relation between delivery time and capacity expansion decisions. The authors developed a mathematical model to determine the optimal delivery time quote under xed and expandable capacities. yExamples include Pizza Hut?s 30 minutes delivery guarantee and FedEx?s 10:30 a.m. delivery guarantee 7 Rao et al. (2005) developed a model for integrating demand and production planning, which determines an optimal planning cycle and a corresponding guaranteed maximum pro- duction lead time. Similar to the previous studies, the model optimizes the expected pro t by quoting a uniform guaranteed delivery time to all customers, however with a particular focus on updating the production schedule in every xed point in time. Di erent than above studies, Duenyas and Hopp (1995) present a non-uniform delivery quote model in which the only controllable variable is accepting or rejecting a customer order. The authors address both nite and in nite capacitated single server queueing systems and demonstrate a control policy with optimal sequencing rules. The second group of studies focuses on ful llment operations which mainly include picking, packing, and shipping processes. de Koster et al. (2007) gives a survey of the design and control of order picking systems. In practice, manual picker-to-part systems include both discrete picking, in which orders are assembled one-by-one, and batch picking, in which multiple orders are picked together. They de ne wave picking as a manual picker-to-part method in which orders for a common destination are released simultaneously for picking. When synchronous picking is performed within small regions of the DC (called zone picking), release for zones is usually controlled with waves. Items in the orders are then consolidated in an automated sortation area (Johnson and Meller, 2002). Speaker (1975), Hu man (1988) and Frazelle and Apple (1994) de ne wave picking as a special case of batch-zone picking, where pickers pick very large batches based not on the number of items, but rather on a length of time. A comparison of di erent order picking policies, including wave picking, is given by Petersen (2000). Franco (2006) distinguishes between batch and wave picking and discusses drawbacks of wave picking from an industrial perspective. When the DC releases waves, it may be faced with the problem of sorting. To alleviate the sorting problem, DCs can use pick list generation algorithms (Owyong and Yih, 2006). 8 Some recent studies focus on an alternative policy called waveless picking (Bradley, 2007), in which orders are sent directly to picking upon arrival. Although applications di er in details (Perry, 2007; Morris, 2008), the objective is to decrease the non-productive walking of order pickers. Gilmore (2006) discusses and compares wave and waveless policies. Gallien and Weber (2010) develop a mathematical model to control waveless operations in order to maximize throughput. As discussed in Chapter 2, our modeling approach is based on uid approximation of an order ful llment system in which orders are released into waves. (Fluid models take the advantage of heavy-tra c conditions in queueing models where it is reasonable to replace discontinuities with continuos functions.) Exact queueing models are not appropriate for a number of reasons: (1) DC environment is too complex to analyze with exact models. Because wave policies imply queues and delays during or after the release times (i.e. rush hour), exact queueing analysis of such cases is di cult even for the simplest assumptions. (2) We are interested in systems with stationary and non-stationary arrivals. Arbitrary, non-stationary processes are typically beyond the reach of exact queueing models (Gupta et al., 2006). (3) The system dynamics and the performance measure we focus on cannot be analyzed by pure queueing models. The use of uid approximations to queueing models has been investigated by many researchers. The majority of studies have analyzed the behavior of queueing systems in heavy-tra c conditions. We present only studies relevant to our work. The rst model was developed by Newell (1973) to approximate the stochastic behavior of n-server service systems with a large n. Borovkov (1964, 1965) discussed limit theorems for mass service as well as for large values of waiting times. Ridley et al. (2004) implemented a uid approximation to a priority call center in which the system receives time varying arrivals. The authors investigate the e ectiveness of the approximation by comparing it with discrete event simulation. 9 Liu and Whitt (2010) introduced a deterministic uid model that serves as an approxi- mation for the Gt=GI=st+GI many-server queueing model. The system has a single class of time varying arrivals, general service time distribution (in parallel), time dependent number of servers (st) and abandonment time distribution. The uid model is intended to serve as an approximation for the queueing model when both the number of servers and the arrival rate are large, and then the system experiences occasional periods of signi cant overload- ing. Such approximations are usually justi ed by many-server heavy-tra c (MSHT) limits in which arrival rate and number of servers are increased with the same scaling factor (see Pang and Whitt, 2008 for details). Fluid models allow us to approximate the dynamics of a queueing system especially when there is a uctuating (non-stationary) arrival stream. A very detailed study on uctuating load was given by Gupta et al. (2006), in which system performance was studied under high loads. The authors addressed the performance experienced by customers that arrived to the system when it is in high or low load periods. Perry and Whitt (2011) analyzed multi-class and pools in a series of papers. The authors followed the same scaling methodology | using the standard many-server heavy- tra c scaling; such that both arrival rate and the number of servers are scaled up by a factor of n. Ward and Bambos (2003) identi ed the stability of queueing networks with deadlines. Each arriving job has a deadline, and a single server processes stationary arrivals under rst-come- rst served discipline. The job abandons the system if it is not completed by its deadline. The authors established stability conditions for this system. More studies can be found in Dieker (2006) for uid approximations and extremes in queueing models, Kella and Whitt (1996) for structural properties of storage networks, and Kella and Whitt (1998) and Whitt and Liu (2011) for very general queueing network approximations with uid models. 10 This research contributes to the deadline oriented order ful llment literature in two ways: (1) Our research has an emphasis on organizing internal operations to maximize a service objective. With the exception of Doerr and Gue (2011), previous research has focused on traditional measures such as throughput, ow time, or work-in-process inventory. (2) We know of no scienti c studies that have investigated systematic wave planning in the presence of deadline-oriented operations. Although wave picking is common in practice, there is still a gap between practice and academic research. Most of the previous academic studies in order picking have focused on e cient picking policies. 1.5 Problem Statement Our research considers ways in which service performance is maximized by releasing the ow of work optimally to follow on processes. To improve service performance, we propose to explore wave release policies in a DC that operates against a daily deadline. We believe this research is the rst comprehensive study of wave planning in a distribution environment, and certainly the rst that addresses deadline-oriented operations. We begin with a distribution environment in which the cuto time is common for all customers and ask, Problem 1. What is the optimal timing and number of waves in a DC that operates against a single daily deadline? Consider a DC that operates against a daily deadline (e.g. an internet retailer that ships all orders overnight). An important operational decision for this system is when to release orders to pickers. How should managers establish these times? Should they be equally distributed through the day, or does another pattern provide better service? Does the timing really matter from a customer service perspective? Our purpose is to answer these and related questions, all with the goal of maximizing service against a daily deadline. In our rst research question, we address ful llment systems in which both the arrival process and the server capacity are known. When one of these quantities (or both) is 11 uncertain, the utilization of the server becomes a random variable. What should be the release times if the server?s utilization is unknown? What would be the expected NSD? What is the risk of planning waves against low (or high) utilization? That is, Problem 2. How should the timing and number of waves be determined when utilization of the server is uncertain? In the rst two research questions, we address systems with a single order class. In a more complex ful llment environment, there are multiple order classes, corresponding to customer groups receiving service at di erent frequencies. This can be interpreted in our context as orders having di erent next scheduled deadlines, but the order ful llment system may choose to work any orders it likes in a wave. Problem 3. What should be the timing and the content of waves to maximize service per- formance across multiple classes of orders? This problem is considerably more complex because we must address not only the timing and the number of waves, but also the contents of waves. When a wave is released, it will include orders for the most imminent departure, but may also include orders for future deadlines. The advantage of doing so is obvious, but there is also a disadvantage. An order being worked early consumes capacity that might otherwise be used for an upcoming more urgent order. How to release orders into waves and when to release those waves must be considered for an optimal solution to this problem. 1.6 Organization of the Dissertation The remainder of this dissertation is organized as follows. In Chapter 2, we introduce optimal wave release strategies for order ful llment systems in which there is a single class of orders. We present the optimal release times for single and multiple wave release systems, and show how 100% service can be achieved by o ering optimal cuto times. We determine the optimal number of waves when there is an associated xed time component to waves. 12 We also present the underlying assumptions behind our models and validation of the uid model via a discrete event simulation model in this chapter. The models in this chapter constitute a foundation for analyzing more complex systems. In Chapter 3, we discuss how to set release times when the workload is uncertain. We show how to adjust single and multiple waves for a given density function of utilization. We also present a procedure to adjust the wave release times and the cuto time when the density function is only known empirically. In Chapter 4, we introduce systems of multiple order classes in which there are di erent deadlines every day. We rst address a basic system with multiple order classes and a single wave.We extend our discussion to multiple class, multiple wave systems and present our solution approach. We illustrate the use of the models with numerical examples and discuss how to implement those models in real applications. We o er conclusions and directions for future work in Chapter 5. 13 Chapter 2 Optimal Release Strategies for Order Ful llment Systems with Deadlines 2.1 Introduction Consider a DC that receives orders 24 hours a day, 7 days per week, serving a large enough area to require multiple days of transport to distant customers. Each customer receives service at the same frequency. We assume there is a single deadline such that each arriving order is assigned to a constant next scheduled deadline (NSD) time within a day. We assess the service performance of the DC with NSD. The DC has a xed capacity which depends on the workforce and which is su cient to process all arriving orders. Arriving orders accumulate in a bu er until the next wave is released, at which time the quantity of orders in that wave decreases at rate until the wave is complete. Upon completion of a wave, the server may go idle, or another wave can be released. The release of a wave before completion of the current wave (overlapping) is disallowed. While the server is working a wave, orders in the next wave accumulate, and the cycle continues. Figure 2.1 illustrates the inventory graphs for this system when there are three waves. Without loss of generality, we assume a day starts (and ends) at a deadline. That is, the start (and end) time of a day is equal to 0 (and 1). For now, we assume the cuto times are equal to deadlines. The workload in a wave is determined by the time since the most recent wave was released. For example, if orders arrive at constant rate of and there are three waves at times 0.3, 0.6, and 0.8 in a day of length 1, the workloads for each wave are (1 0:8) +(0:3 0) = 0:5 , (0:6 0:3) = 0:3 , and (0:8 0:6) = 0:2 . The rst wave consists of orders arriving between the last wave from yesterday until midnight and orders arriving from midnight until 14 di 1 diw1 w2 w3 Figure 2.1: A single class, three wave system. di indicates the deadline on day i; wj is the release time for wave j. the release of the rst wave today. The objective is to determine wave release times that maximize the fraction of orders completed by the deadline. 2.2 Single Wave Systems Consider a single wave system with wave release time w1. Day 0 is e ectively a \warm up day." Denote the release time of the wave by w1 and constant arrival rate by orders per day. Then, w1 orders arrive before the release on Day 0. The server nishes processing in w1 = time. Waves on all following days contain orders, so the server works = = of the time. NSD in a single wave system is a function of release time w1, arrival rate , and server capacity . By de nition, NSD = # orders worked that arrived in the current cycle and completed before di# orders that arrived in the current cycle : Because cycle de nes unit length of day between successive cuto times, the denominator is . The numerator is more complicated. Each wave consists of orders that arrived the previous day (but after the wave release, and therefore did not make the deadline) and orders that arrived during the current day. (Recall that we assume cuto times equal to deadlines). Only the latter count toward NSD, but we assume the former must be worked 15 before the current day?s work (i.e., we assume rst-come, rst-served discipline). Therefore, the number of orders released today has two components: (1 w1) orders from yesterday and w1 orders from today. If the server nishes the wave before the deadline, then NSD = w1 = w1: (2.1) If the server nishes the wave after the deadline, only (1 w1) orders are processed before the deadline, of which (1 w1) do not count toward NSD because they arrived the previous day. In this case, NSD = (1 w1) (1 w1) = (1 w1)( ) = (1 w1) 1 1 ; (2.2) where = = . The server nishes exactly at the deadline when w1 = 1 ; therefore, NSD = 8 >< >: w1; 1 w1 (1 w1)(1= 1); 1 w1: (2.3) and the maximum possible NSD occurs when w1 = 1 , for any value of . Figure 2.2 illustrates how NSD changes with the release time w1 for three values of utilization. We have just shown that Proposition 2.1. For a system with a single wave, the optimal release time w 1 = NSD = 1 : The proposition con rms two points of intuition. First, the server should begin work as late as possible in order to allow as many orders as possible to make it into the wave. Second, the wave should nish exactly at the deadline. 16 NSD w10 1 10.75 0.75 (a) = 0:75: NSD w10 1 10.5 0.5 (b) = 0:5: NSD w10 1 10.25 0.25 (c) = 0:25: Figure 2.2: Di erent realizations of optimal w1. Recall that NSD is de ned as the fraction of orders arriving between two consecutive cuto times that nish by the next deadline. As Doerr and Gue (2011) note, NSD is merely an accounting measure that can be manipulated by shifting the cuto time; that is, it can be increased by making the cuto time earlier, or decreased by making it later. In neither case do customers receive their orders any sooner or later. Nevertheless, the cuto time can serve at least two purposes; if it is published, the cuto time establishes expectations for customers (e.g., Amazon?s overnight delivery guarantee); if not, it can be used as an internal metric to motivate workers (Doerr and Gue, 2011). For our purposes, an optimal cuto time is the latest possible time for which NSD is 100 percent. Proposition 2.1, then, suggests Proposition 2.2. In a deterministic order ful llment system with a single wave, the optimal cuto time equals the optimal wave release time w1 = 1 . 2.3 Multiple Wave Systems In practice, distribution centers typically use multiple waves per day|usually 2{6, de- pending on the workload and number of destinations that must be accommodated. Our model of a multi-wave system is built on the simple insight that the number of orders, or load, in a wave is the product of the arrival rate and the time since the previous wave 17 was released. (We require that waves consist of all orders available to be released.) Before developing equations for a multi-wave system, we require two results. Proposition 2.3. In an optimal solution, work on the nal wave ends at the deadline. Proof. Let wN be the last wave release time in an N-wave system. When the last wave ends at the deadline, (1 wN) orders do not make the deadline and wN do. Therefore, NSD = wN = wN: Now assume to the contrary of the proposition that there exists a feasible set of release times such that wN does not nish at the deadline and the corresponding NSD wN. There are two cases: Case 1: Assume there is an optimal solution in which the nal wave wN nishes t before the deadline (0 < t < 1 wN), with NSD = x. Now, shift every wave release t to a later time. The system is still feasible, and now NSD = x+ t; a contradiction. Case 2: Assume there is an optimal solution in which the nal wave wN nishes after the deadline, such that NSD wN. Because the nal wave nishes after the deadline, the number of orders arriving today that make today?s deadline is the number worked before the deadline (1 w1) minus the number that arrived after the last wave in the previous day (1 wN) . Therefore, NSD = (1 w1) (1 wN) ; which can be rewritten as NSD = wN + 1 w 1 1: (2.4) Because the server begins working at w1, the number of orders completed before the deadline is no greater than (1 w1) , which must be less than total arrivals in a day . That 18 is, (1 w1) < (1 w1) < 1 1 w1 < 1: Substituting into Equation 2.4 implies NSD< >: j 1 N ; for = 1 j N+1 1 N ; for < : (2.6) Figures 2.3{2.5 show the optimal wave release times for di erent levels of utilization. Each horizontal line corresponds to a di erent system, with the indicated number of waves. Dots on the line correspond to optimal release times. Because there is no idle time between waves (Proposition 2.4), the time of the rst wave w1 does not change, but later wave times adjust as the number of waves increases. In fact, Equations 2.5 suggest a relationship between 20 d0 d1 Waves NSD 1 5 2 6 3 7 4 8 0.83 0.99 0.93 0.99 0.97 0.99 0.98 0.50 Figure 2.3: Optimal wave release times: = 0:5: d0 d1 Waves NSD 1 5 2 6 3 7 4 8 0.25 0.92 0.68 0.95 0.82 0.96 0.88 0.97 Figure 2.4: Optimal wave release times: = 0:75: d0 d1 Waves NSD 1 5 2 6 3 7 4 8 0.05 0.83 0.54 0.86 0.70 0.88 0.78 0.90 Figure 2.5: Optimal wave release times: = 0:95: 21 consecutive release time intervals: wn+1 wn wn wn 1 = ; which means that successively earlier wave lengths are 1= larger. This can be seen in Figures 2.3 and 2.4 especially. Recall that only orders arriving after the nal wave will not make the deadline, and therefore NSD = wN. In the expression above, substituting N for j and simplifying gives the optimal NSD for a system with N waves and utilization , NSD = wN = 1 N(1 ) 1 N : (2.7) As expected, wN = 1 when N = 1, and wN (and NSD) converges to 1 as the number of waves N!1. Figure 2.6 illustrates how NSD varies for systems with one to ve waves. As utilization increases, the maximum possible NSD decreases, converging eventually to (N 1)=N (Equation 2.6). The plot on the right shows how many hours before the deadline the cuto time should be set in order to achieve 100 percent NSD. 0.2 0.4 0.6 0.8 1.0 0.2 0.4 0.6 0.8 1.0 1.0 0.8 0.6 0.4 0.2 20 15 10 5 NSD Hours bef or e deadline 5 4 3 2 1 Utilization Utilization 5 4 3 2 1 # w aves Figure 2.6: Two perspectives on maximum possible NSD for di erent numbers of waves. The plot illustrates that an aggressive service promise via a late cuto time can be kept only by increasing the number of waves, or by adding capacity (decreasing utilization), or 22 both. The plot also shows that there is little marginal bene t to increasing the number of waves beyond four or ve. 2.4 Optimal Number of Waves We have shown that the last wave release time wN and therefore NSD converges to 1 as the number of waves N !1. Figures 2.3{2.5 show that there is little marginal bene t to increasing the number of waves beyond four or ve. From Equation 2.5, for j > 2, we have (wj wj 1) = (wj+1 wj): Rewriting this expression gives us the wave-size (or workload) ratio of consecutive waves: Lj+1 Lj = = ; where Lj denotes the workload of the j -th wave. More generally, the workload ratio of any two waves in a system with multiple waves can be written as: Lj+n Lj = n: Wave size ratios depend on the assumption that processing rate is constant without regard to the size of the wave. However, for many order picking systems, there is a xed time component to a wave. For example, in a manual picking system, workers must walk a tour to gather items in the wave. Thus, the processing time is comprised of a xed time T to walk the tour and a variable processing time Lp that depends on the number of picks in the wave L and the time p to make a single pick. The processing rate can be calculated by dividing the total number of picks with the total time required to traverse the picking area 23 and pick orders. When the number of orders in a wave is equal to L, the process rate is (L) = LT +Lp: Modifying Equation 2.5 with (L), T +p (1 wm 1N +wm1 ) = wm2 wm1 ; ... T +p (wmn wmn 1) = wmn+1 wmn; (2.8) ... T +p (wmN wmN 1) = 1 wmN: We use the RSolve function in Mathematica to solve for the optimal release times for day m and determine the steady state release times for an N-wave system. The right-hand side of the above equations refers to the allowed time of a wave in order to complete the work-content p (wmn wmn 1) including the xed time T (left-hand side). This observation leads to: Proposition 2.5. When processing rate is a function of walk time, pick time and wave size, the optimal number of waves N = 1 p T : Proof. For a system with N-waves, the rst wave is released at w1 = 1 NT p. Because w1 0, 1 NT p 0, and therefore the maximum number of waves is N = 1 p T : 24 Proposition 2.5 leads us to the plots in Figure 2.7, where we illustrate di erent levels of picking volumes when T = 0:05 (i.e., 1.2 hours of total xed time including operations in picking, packing, and shipping). 2 4 6 8 10 12 14 Waves 0.6 0.7 0.8 0.9 NSD ? p = 0 .4 ? p = 0 .6 ? p = 0 .8 Figure 2.7: Change in NSD with respect to the number of waves when there is a xed time component associated with each release. As p increases, the optimal number of waves decreases. For example, when p = 0:4, the optimal number of waves is N = 14 and maximum NSD=0.93. (However, there is little marginal return beyond 5 waves). Suppose p = 0:8, then the optimal number of waves is only 4 and the maximum NSD is 0.77. In this case, workers either require more time to pick items or must to pick more items, both resulting in longer processing times and therefore fewer waves. 2.5 Validation of Fluid Approximation In the previous sections, we assume a uid model, which assumes a continuous stream of work that ows into the system. In a typical order ful llment setting, however, orders arrive at discrete times. Therefore, we should determine when the uid approximation is valid. Fluid models as approximations to queueing systems have an extensive literature (see especially Bernd S. et al., 2011; Dai, 1995; Dai and Jennings, 2003). 25 Recall that, the accuracy of the uid approximation improves as the arrival and service rate or the number of servers becomes large (Section 1.4). The typical method of validating a uid approximation is to increase the number of servers and the corresponding arrival rate while keeping the utilization constant (Pang and Whitt, 2008). The arrival rate of a system with k servers is k = k . We use a similar but di erent approach. Because our interest is not the number of servers, but the length of the job (i.e., the service rate), we scale up the service rate parameter and keep the number of servers constant, thereby maintaining a constant utilization. dwN (a) 5 orders each of which takes two hours of processing. dwN (b) 20 orders each of which takes half hour of processing. dwN (c) In nitely many number of or- ders completed instantaneously. Figure 2.8: Accuracy of the uid approximation improves as orders become more indivisible. Figure 2.8 illustrates the approach. In the leftmost example, the nal wave of the day begins two hours before the deadline, in a system with ve workers. The nal wave contains ve orders, each taking two hours to complete. The work content, then, is 10 order-hours. The middle gure shows the same amount of work content spread among 20 orders, each taking 30 minutes to complete. When the wave is released, 15 orders wait in the queue until a server is free. The uid model is a limiting case in which the work content is divided into in nitely small portions, which are then completed in a continuous stream of output. If the time to process an order is stochastic, we should expect that some orders will nish after the deadline, and that the expected number of orders missing the deadline will vary based on how the work content is modeled|or said another way, based on how long 26 the expected processing time is with respect to the length of the nal wave. By contrast, the uid model assumes work is deterministic and that all orders in the nal wave make the deadline. Is this a problem? To answer this question, we simulated a three-stage order ful llment system, corre- sponding to the picking, packing, and shipping functions in a DC. To illustrate the point, consider a ful llment system in which each stage has twenty servers, representing workers in those functions. We consider three levels of utilization: 0.5, 0.75, and 0.95. For a four wave system, the lengths of the nal waves are equal to 24 (1 w4) = 56 minutes, 2.78 and 5.27 hours for = 0:5; 0:75; and 0:95, respectively. The goal is to determine at what ratio of expected processing time to length of nal wave the uid model breaks down, and whether common practice is typically less than that ratio. We begin with a very long expected processing time of 1,024 minutes (17 hours), then half the processing time in successive experiments until expected processing time is just 1 minute. Processing time is divided equally among the three stages, and processing times in each stage are exponentially distributed. We adjust the (exponential) arrival rate to maintain the appropriate utilization. Runs last 30 simulated days, with 3 days of warm-up and 25 replications. For each run, we compare the average simulated NSD with the uid model approximations, which were 96.7%, 88.6% and 78.1% for = 0:5;0:75;0:95. We show how to nd wave release times and resulting NSD levels in Section 2.3. Table 2.1 shows the results. As expected, the model performs very poorly when expected processing time is longer than the length of the nal wave (ratio greater than 1). When the ratio is less than about 0.3, the uid model approximation is within one percent of the simulated value. The uid model always overestimates NSD because it assumes full utilization of the (single) server and deterministic service time. In the simulated (and a real) system, the nal wave sometimes nishes early, but NSD is no higher because no more orders get processed before the deadline. However, sometimes the nal wave nishes late and orders fail to meet the deadline. NSD in this case is lower than predicted by the uid 27 Table 2.1: Validation of uid model via simulation. = 0:50 = 0:75 = 0:95 Processingtime E[NSD] error(%) ratio E[NSD] error(%) ratio E[NSD] error(%) ratio 3,072 1.8 94.8 53.9 2.6 85.9 18.4 3.2 74.9 9.7 1,536 8.4 88.2 26.9 13.2 75.3 9.2 12.8 65.9 4.9 768 29.5 67.1 13.5 37.4 51.2 4.6 37.2 40.9 2.4 384 60.9 35.8 6.7 63.0 25.6 2.3 62.7 15.4 1.2 192 79.4 17.2 3.4 79.4 9.2 1.2 74.1 3.9 0.6 96 88.8 7.9 1.7 86.6 1.9 0.6 78.0 0.1 0.3 48 93.8 2.9 0.8 88.3 0.2 0.3 78.1 0 0.2 24 96.2 0.5 0.4 88.4 0.2 0.1 78.1 0 0.1 12 96.6 0.1 0.2 88.5 0.1 0.1 78.1 0 0 6 96.7 0 0.1 88.5 0.1 0 78.1 0 0 3 96.7 0 0.1 88.5 0 0 78.1 0 0 model. In our experience, a ratio of average processing time to length of nal wave of about 1/3 is realistic for most distribution environments. 28 Chapter 3 Setting Wave Release Times in the Presence of Uncertainty 3.1 Introduction In the previous chapter, we discussed how to set release times when the workload is known. Most order ful llment systems, however, experience workload uctuations from day to day. Some uctuations are cyclical. Think of weekly patterns, such as peak workload on Sunday as experienced by online retailers (Bates, 2012). Because of workload uncertainty, planning operations based on an average daily workload results in poor performance. Work- load can also uctuate during the day as, for example, customers place orders in the morning or after work. Consider a DC that has been designed to accommodate a certain workload. There are two potential consequences: (1) If the DC has to process signi cantly higher volumes, then it will not nish all orders in time, resulting in lower NSD. (2) If the workload is lower, then the system completes orders before the deadline, and labor resources are underutilized. In addition to workload uncertainty, worker absenteeism and other causes of variable capacity lead to uncertainty in utilization. We have shown that wave release times (and the number of waves) can have a signi cant impact on NSD for a given utilization level. But how should the release times be set when utilization is uncertain? What is the risk of releasing waves too early or too late? 29 3.2 The Case of a Single Wave We start by describing how to adjust a single wave. For a xed release time w1, we showed that (Equation 2.3), NSD = 8 >< >: w1; 1 w1 (1 w1)(1= 1); 1 w1: Figure 3.1 shows NSD as a function of utilization for three speci c values of w1. NSD 0 1 10.25 0.75 0.50 0.25 0.50 0.75 w1 = 0:75 w1 = 0:50 w1 = 0:25 Figure 3.1: Service performance versus utilization , for three wave release times. When w1 = 0:5, for example, NSD is at its maximum 50% for all values of 0:5. In such cases, the server would nish at or before the deadline, but there is no increase in performance because all orders arriving after the release still miss the deadline. For > 0:5, performance on NSD drops o sharply because many orders in the wave nish after the deadline. The curves illustrate that, in the presence of uncertainty in workload or capacity (i.e., uncertain ), setting the release time too late (high w1) risks a signi cantly lower performance on NSD if utilization is actually high, but setting it too early (low w1) limits the highest possible NSD if utilization is actually low. Consider two systems|one with a xed utilization 0.5 and one with utilization dis- tributed Uniform[0,1]. The optimal release time for the rst system can be calculated from 30 Proposition 2.1: w1 = 1 = 0:5. To compute the maximum expected NSD in the second case, we use Equation 2.3 with density function ff( ) = 1; 0 1g. E[NSD] = Z 1 0 NSD(w1; )f( ) d (3.1) = Z 1 w1 =0 w1f( ) d + Z 1 =1 w1 (1 w1) 1 1 f( ) d = (w1 1) loge(1 w1): Taking the derivative and setting equal to zero, d dw1E[NSD] = loge(1 w1) + 1 = 0; which has solution w1 = 1 1e 0:6321; which is later than in the deterministic case. From Equation 2.3, maximum NSD for this release time is approximately 0.6321, but expected NSD is (1 1=e 1) loge(1 1 1=e) = 1=e 0:3679. Expected NSD for w1 = 0:5 would have been 0.3466. This simple example shows that releasing a single wave later than it otherwise would have been released improves expected NSD. 3.3 Adjusting Multiple Wave Release Times That uncertainty should lead to a later release time seems to violate the intuition that a later release time is \riskier." As we now show, it is not true that uncertainty of utilization leads to a later optimal release time for the more practical case of multiple wave releases per day. To see why, consider the optimal wave release times for many values of (call them ^ ). We can then compute expected NSD for each set of times, given a density of utilization f( ). Because the observed utilization may not equal the planned utilization ^ , we must 31 develop expressions for NSD for two possible cases. If < ^ , the last wave nishes early and there is idleness between consecutive waves. Only orders arriving after the last wave miss the deadline, so NSD = wN: If > ^ , the rst wave does not nish before the release of the second; the second does not nish before release of the third; and so on, and the last wave does not nish before the deadline. From Equation 2.4, NSD = wN + 1 w 1 1: From Proposition 2.3, we know w1 = 1 ^ . Therefore, ^ = 1 w1 and NSD = 8 >< >: wN; 1 w1 = ^ wN + 1 w 1 1; 1 w1 = ^ : (3.2) Figure 3.2 illustrates the relationship between planned utilization ^ and observed utiliza- tion . Planning for a high utilization means an earlier rst release time w1, and therefore a lower maximum NSD; planning for low utilization allows a higher NSD if observed utilization is less than planned, but increases the chance that observed utilization will exceed planned. Which curve provides the highest expected NSD depends on the density of utilization f( ). The curves in Figure 3.2 are easy to generate, and for a given f( ) we can easily integrate them to compute expected NSD. To illustrate the procedure, we assume utilization is distributed Uniform[0,1] and use Equation 3.2 to compute NSD for many values of planned utilization ^ . Each curve in Figure 3.3 is a piecewise-linear representation of expected NSD for 100 values of ^ . The maximum point for a single wave system is E[NSD] = 0:3679, which equals the optimal value of the single-wave model as expected. Release times for the best four-wave solution are w1 = 0:32; w2 = 0:597; w3 = 0:785; and w4 = 0:913; which produce E[NSD] = 0:855. Had we used the deterministic solution for = 0:5, release times would have been w1 = 0:5; w2 = 0:767; w3 = 0:9; w4 = 0:967, with E[NSD] = 0:813. In other words, adjusting 32 0.2 0.4 0.6 0.8 1.0 0.2 0.4 0.6 0.8 1.0 NSD Observed utilization Planned utilization 0.9 0.8 0.7 0.6 0.5 0.3 0.2 0.1 0.4 0.01 Figure 3.2: Each curve shows how NSD changes for di erent levels of observed utilization with respect to planned utilization. release times for the uncertain utilization increases expected NSD by 4.2 percent. Table 3.1 0.2 0.4 0.6 0.8 1 Planned utilization 0.2 0.4 0.6 0.8 NSD 8 . . . 3 2 1 # waves Figure 3.3: Maximum expected NSD achieved when release times are adjusted to earlier times. shows the percent improvement in expected NSD by adjusting the release times for other numbers of waves. 33 Table 3.1: Adjusting release times in the presence of workload uncertainty improves expected NSD. Number of waves 1 2 3 4 5 6 7 8 E[NSDj =0.5] (%) 34.7 67.9 77.5 81.3 83.0 83.9 84.3 84.5 E[NSD] (%) 36.8 68.4 79.8 85.5 88.9 91.2 92.7 93.9 Improvement (%) 2.1 0.4 2.3 4.3 5.9 7.3 8.4 9.4 0.2 0.4 0.6 0.8 0.95 Plannedutilization 0.2 0.4 0.6 0.8 1.0 NSD 0.2 0.4 0.6 0.8 1.0 Observed utilization 0.2 0.4 0.6 0.8 1.0 NSD Figure 3.4: Ten di erent levels of planned utilization for the case study (left) and corre- sponding average NSD (right). 3.4 An Empirical Procedure to Adjust Release Times In previous section, we showed how to adjust release times for a known density function of utilization. In practice, of course, the density is unknown, so variability of utilization must be handled in another way. Here we describe an empirical procedure to increase expected performance on NSD, given historical utilization data. (We used 80 day workload and shipment data of an industrial DC to analyze the variability of utilization). First, we generate performance curves for many levels of planned utilization, as in Figure 3.2. Each curve corresponds to a di erent set of release times and tells what NSD would be for any value of observed utilization. Next, for each value of planned utilization (and its associated curve) we simply record what NSD could have been for each day in the data set. The planned utilization giving the highest average NSD gives us the best set of release times. The left plot in Figure 3.4 shows curves corresponding to ten di erent levels of planned utilization. Data points on the bottom correspond to the observed utilization for each day 34 in the data set. The right plot shows average NSD for the ten levels of planned utilization. The highest value of NSD is 80.4%, which corresponds to a planned utilization of 80% and release times w1 = 0:200; w2 = 0:424; w3 = 0:690, and w4 = 0:841. As expected, the best value of planned utilization when accounting for variability of utilization is greater than the average utilization in the historical data. Higher planned utilization means earlier wave releases and a reduced chance that the nal wave will not nish before the deadline. 3.5 Setting the Cuto Time To this point we have assumed that the cuto time equals the deadline, which is hardly realistic if the cuto time is to be published to customers. How should a cuto time be established such that orders make the deadline \almost always?" Note that all orders arriving after the nal wave will miss the deadline, so the cuto time should be no later than wN, but setting it earlier than wN means the nal wave contains orders not due at the upcoming deadline. We conclude, then, that the cuto time should equal wN. Why not release the nal wave earlier to get a headstart on the nal orders arriving for the next deadline? De ne the cuto time setback to be the time between the cuto time and the deadline. The left plot in Figure 3.5 shows daily performance on NSD for the data set when release times are set according to ^ = 0:8 and the cuto time is equal to wN, for a setback of 3.8 hours. All customer requests were met on 53 of the 80 days, or 66.3% of the time, which is analogous to Type 1 service in inventory theory. The same policy yields a ll rate (Type 2 service) of 97.24%. The plot on the right shows Type 1 and Type 2 service levels for di erent levels of ^ . In each case, the cuto time is set to the new nal release time wN. 35 20 40 60 80 days 0.2 0.4 0.6 0.8 1.0 NSD SolidCircleSolidCircle SolidCircle SolidCircle SolidCircle SolidCircle SolidCircle SolidCircle SolidCircle SolidCircle SolidSquare SolidSquare SolidSquare SolidSquare SolidSquare SolidSquare SolidSquare SolidSquare SolidSquare SolidSquare 1 2 3 4 5 SetbackLParen1hrs.RParen1 0.2 0.4 0.6 0.8 1.0 Figure 3.5: Performance results when the cuto time is adjusted. On the left, NSD for each of the 80 days in the data, sorted from high to low. On the right, Type 1 and Type 2 performance for di erent levels of cuto time setback. 3.6 Adjusting The Number of Waves Proposition 2.5 shows how to determine the optimal number of waves for a given uti- lization . When is uncertain, N becomes uncertain and determining an optimal policy is more di cult. Consider a ful llment system in which planned utilization is less than observed utiliza- tion. High utilization implies releasing fewer waves; however, Equation 2.6 suggests more waves to achieve the same level of NSD. As a consequence, there is a risk of releasing more waves and missing some portion of the orders. Our objective is to determine the number of waves for a ful llment system in which the utilization is uncertain. We rst consider a system with a single wave release. As in Section 2.4, the xed time component is T. If the server completes the wave before the deadline, then the wave will be completed at w1 + = + T. Because the wave is completed before the deadline, w1 orders count toward NSD, and NSD = w1. If the server nishes the wave after the deadline, (1 w1 T) orders are processed, of which (1 w1) do not count toward NSD. In this 36 case, NSD = (1 w1 T) (1 w1) = (1 w1)( ) T = (1 w1) 1 1 T : Consequently, NSD = 8 >< >: w1; 1 w1 T (1 w1) 1 1 T ; 1 w1 T which is maximized when w 1 = 1 T. As expected, the optimal release time is now T earlier and as it increases, the wave should be released earlier, resulting in a lower NSD. When the workload is uncertain, the expected NSD can be calculated from E[NSD] = Z 1 0 NSD(w1; ;T)f( ) d : (3.3) To nd the optimal release time of a single wave, we take the derivative of the above expres- sion with respect to w1 and set it equal to zero. To illustrate the procedure, consider the following example. Assume that the server?s utilization is distributed Uniform[0,1]. What is the improvement in expected NSD when the release time of a single wave is adjusted based on the xed time wave component? Using the above equation, E[NSD] = Z 1 0 NSD(w1; ;T)f( ) d = Z 1 w1 T =0 w1 d + Z 1 =1 w1 T (1 w1) 1 1 T d = (w1 +T 1)log w 1 +T 1 T 1 and E[NSD] is maximized when 37 w 1 = 1 1e (1 T): When is uncertain, the adjusted release time is always later than what it would otherwise be: w1 1 ^ T. We nd the expected NSD by inserting w 1 into Equation 3.3: E[NSDjw 1] = 1 Te : The expected NSD is a linear function of T for the optimal value of w 1. Figure 3.6 shows that releasing a single wave later than it otherwise would have been released improves expected NSD. ELBracket1NSD wEqualLParen11Minus1Slash1eRParen1LParen11MinusTRParen1RBracket1 ELBracket1NSD wEqualLParen11Minus1Slash1eRParen1RBracket1 0.2 0.4 0.6 0.8 1.0 T 0.05 0.10 0.15 0.20 0.25 0.30 0.35 ELBracket1NSDRBracket1 Figure 3.6: Expected NSD considering only uncertainty (bottom curve) and expected NSD considering both uncertainty and the xed time component (top curve). For the above example, although considering the xed time component does not produce signi cant improvement (1.21%) when we account for multiple stages (e.g. packing and shipping), it produces more than 3.5%. Now, consider a multiple wave system. A lower server utilization suggests releasing more waves (lim !0N = 1=T). A higher utilization, on the other hand, implies fewer waves. What should be the optimal number of waves when the utilization is uncertain? 38 As before, when waves are completed before the deadline NSD = wN, otherwise the pickers process (1 w1 NT) many orders, of which (1 wN) of them arrive after the last release and thus do not count toward NSD. The resulting NSD is [ (1 w1 NT) (1 wN)] = wN + (1 w1 NT)= 1. Consequently, NSD = 8 >< >: wN; 1 w1 NT = ^ wN + 1 w1 NT 1; 1 w1 NT = ^ : (3.4) For given values of w1, wN, we can calculate NSD; however, for di erent values of , the number of waves N changes. To determine the optimal N that maximizes expected NSD, we follow a similar but slightly di erent procedure than the one in Section 3.3. The steps of the procedure are (1) For a given T and f( ), determine ^ = E[ ]. Determine N =b(1 ^ )=Tc. (2) Use the recursive system given in Equation 2.8 to determine w1, wN and NSDN. (3) Insert w1, wN into Equation 3.4 and use it in the following steps. (4) Calculate the loss ?( ) associated with releasing N waves. (4.1) Calculate down-side loss ?d( ): For ^ evaluate Equation 2.8 and determine resulting NSDd, ?d( ) = NSDd NSDN. (4.2) Calculate up-side loss ?u( ): For ^ , repeat (3) and determine resulting NSDu, ?u( ) = NSDN NSDu. (4.3) When the actual utilization is , total loss ?( ) = ?d( ) +?u( ). (5) Determine expected NSD and expected loss of using N waves. (5.1) Expected NSD: E[NSDjN] = Z 1 =0 NSD(w1;wN; )f( )d : 39 (5.2) Expected loss of using N waves: E[?jN ] = Z 1 =0 ?( )d : (6) Calculate expected NSD after expected loss E[NSDjN] E[?jN]. (7) While N > 0, decrement N by one and go to step 2. (8) Select the N that produces the maximum E[NSD]. We illustrate the procedure with the following example. Assume that the utilization is distributed Uniform[0,1] and T = 0:1. Proposition 2.5 suggests N = b(1 0:5)=0:1c = 5. When = 0:5, the rst and the last release times are w1 = 0; w5 = 0:8 resulting in NSD=0.8. Suppose that observed utilization is lower than planned, to say = 0:25. Although the maximum number of waves should beN =b(1 0:25)=0:1c= 7, because the number of waves were decided based on the planned utilization, the number of waves would still be ve. A lower observed utilization, however, suggests di erent release times: w1 = 0:25; w5 = 0:8664 (NSD=86.6%). Consequently, there will be a 6.64% loss in NSD due to releasing ve waves, but not adjusting the release times (downside loss). Now suppose that observed utilization = 0:75. The maximum number of waves should be N = b(1 0:75)=0:1c = 2. The resulting NSD would be 0.536. Because the maximum number of waves is two, releasing ve waves must produce a lower NSD. When there are ve waves, total xed time is 50% of the day, resulting in late completion of the waves (because remaining capacity will not be enough to complete the waves before the deadline). Releasing ve waves will produce an NSD equal to 0.397. The loss in NSD due to planning for low utilization is 0.536-0.397=0.139 (upside loss). If the planned utilization is equal to observed utilization (^ = = 0:5), then there would be no loss in NSD. For all other values of , there is either downside or upside loss which needs to be considered. We calculate the loss for each value of with increments of 0.1 (Figure 3.7). 40 N =1 N =2 N =3 N =4 N =5 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 10.0 0.2 0.4 0.6 0.8Loss in NSD Observed Utilization Figure 3.7: Loss associated with using N waves. In Figure 3.7, for data points where ^ , there would be a potential increase if release times were adjusted. In this case, releasing more waves decreases the loss (more orders will be completed). When observed utilization is higher than planned, there are too many waves and therefore the picking density is low (due to excessive xed time in operations). As a result, the last wave (in this example, the fth wave) will be completed after the deadline and NSD will be lower than a system with fewer waves. Taking the expectation of ?( ), we calculate the expected loss of releasing a certain number of waves. Subtracting each loss from its expected NSD value, we determine the optimal number of waves that maximizes NSD. For our example maxfE[NSDjN = j] E[?jN = j]g= maxf0;0:228;0:267;0:279;0:454g= 0:454 when N = 5. Recall that we have assumed T = 0:1 which means 2.4 hours per wave. In a more realistic ful llment system, the total xed time in all processes is no longer than two hours. If we assume T = 2 hours, then Proposition 2.5 suggests releasing (1 0:50)=0:083 = 6 waves. The solution for a six wave system suggests w1 = 0 and w6 = 0:837. In this case, the expected NSD E[NSDjN = 6;T = 2 hours] = 0:433. Following the procedure given above, the expected NSD is maximized when N = 2. Release times for two wave system are w1 = 0 and w2 = 0:502 which produce E[NSD] = 0:621 (18.8% higher than a six wave system). 41 In the above example, we assume a known density function for the server?s utilization and demonstrate the use of our procedure to adjust the number of waves. In practice, the density function is not known, therefore, historical utilization distribution should be used to determine the optimal number of waves. The procedure, in this case, will be slightly di erent than the one we present above (summation is required instead of integration). Nevertheless, we can use the approach given in Section 3.4 to determine the expected NSD and insert it into the Step 5 in the above procedure. 42 Chapter 4 Wave Release Strategies for Systems with Multiple Order Classes 4.1 Introduction Our rst and second research questions address systems with a single class of orders. That is, the DC ships all customer orders at the same time. In many order ful llment environments, there are multiple deadlines per day. Di erent deadlines may be the result of multiple delivery modes. For example, a retailer may have a deadline for FedEx, UPS, and even one or more LTL carriers. Some retailers, on the other hand, may control their own shipment operations. Customers farther from the DC may have earlier deadlines (to allow for timely deliveries); customers located closer to the DC may be assigned to later deadlines. Consider a DC that processes three di erent order classes with three waves (Figure 4.1). d1 d2 d3 3 1 2 w1 3 1 2 w2 3 2 1 w3 Figure 4.1: A three-class, three wave system. 43 Each class has a di erent deadline di;i2f1;2;3g. When a wave is released, it includes the most imminent orders, but may also include orders for future deadlines. (For waves 1 and 2, the most imminent orders are class 1 orders, having a deadline d1.) Should waves release all classes of waiting orders, or release only the most imminent ones? When should the waves be released, and what should be the content of the waves? In this chapter, our objective is to determine the optimal timing and mixture of waves in multiple order class, multiple wave systems. Before specifying an initial model, we need to de ne a \system NSD," which accounts for all order classes. Because there are multiple classes of orders (deadlines), each class has its particular NSD: NSDi = # class i orders that arrived in the current cycle and completed before di# class i orders that arrived in the current cycle ; i2f1;:::;Kg. Because the time between successive deadlines is one day, the denominator is i. Therefore, the total workload in a K class system is L = PKi=1 i. The numerator depends on multiple factors: the number of order classes K; arrival rates i;i2f1;:::;Kg, deadlines di;i 2 f1;:::;Kg and the number of waves N. When those parameters and variables are de ned, the system NSD can be written as a weighted combination of each NSDi: NSD = 1L KX i=1 iNSDi : (4.1) We use Equation 4.1 throughout this chapter to calculate system NSD. (Hereafter, NSD denotes the system NSD.) The server?s utilization is equal to the total workload divided by the total processing capacity : = KP i=1 i = L : 44 For now, we assume that the cuto time of each order class is equal to its deadline and that completed orders are instantaneously ready for shipment. In a system with K order classes, waves complete some classes before their deadlines and some after their deadlines. Releasing multiple order classes in waves re ects realistic operations. In practice, pickers pick orders from all classes in their picking tours. Consequently, waves include multiple order classes that are processed together. 4.2 Single Wave Systems With Multiple Order Classes Our approach is similar to the one given in Section 2.2. Because there is a single wave, all orders accumulated between the release times of two consecutive days constitute the total workload. That is, the server processes orders of all classes in a single wave. The simplest problem is to determine a single wave release time when there are two order classes with deadlines d1 and d2. The rst deadline d1 represents the earliest deadline in a day, followed by the second deadline d2. Without loss of generality, we assume d2 = 1. When should the single wave be released to maximize system NSD? We denote the completion time of a single wave by f1. When the wave is released, there are i many class i orders of which i(1 w1) arrived in the previous day and iw1 arrived in the current day. The total number of orders at time w1 is equal to the sum of orders from all classes: PKi=1 i = L. The server starts working at a rate of and completes the wave at f1 = w1 + L = w1 + : (4.2) What is the system NSD if the wave is completed after d1, or d2? The answer depends on its release time. For a two class, single wave system, there are ve possible cases. Case 1. The wave is released before the rst deadline and completes all work before the rst deadline. 45 Figure 4.2 illustrates this case. All orders from both classes are released at w1. The total number of orders released is equal to 1 + 2. With a rate , the server starts processing d2 d2d1w1 f1 1 + 2 Figure 4.2: Case 1: w1 d2. The system NSD is: NSD = 1 1 + 2 1 + 2 1 (d1 w1) + 2 1 + 2 1 + 2 1 (d2 w1) = 1(d1 w1) + 2(d2 w1) ( 1 + 2)2 : (4.9) Case 4. The wave is released after the rst deadline and the server completes the work before the second deadline (Figure 4.5). d2 d2d1 w1 f1 1 + 2 Figure 4.5: Case 4: w1 >d1 and f1 d1 and f1 >d2. Using Equation 4.8, the system NSD is: NSD = 2 1 + 2 NSD2 = 2 1 + 2 1 + 2 1 (d2 w1): (4.11) To demonstrate how each class NSD changes with respect to the single wave release time, consider the following example. Suppose there are two order classes having arrival rates 1 = 2. The rst and the second class of orders have deadlines d1 = 0:5, d2 = 1. The processing capacity of the server = ( 1 + 2)= , where = 0:5. Using NSD functions for each class, we obtain the piecewise functions of NSD1 and NSD2. Figure 4.7 shows each class NSD (and the resulting system NSD) with respect to wave release time w1. 49 0.2 0.4 0.6 0.8 1.0 w1 0.1 0.2 0.3 0.4 0.5 NSD NSD NSD 1 NSD 2 Figure 4.7: Change in the system NSD in a two class, single wave system. Figure 4.7 shows the change in system NSD when = 0:5. As expected, the function of NSD2 is the same as in single class, single wave systems. When w1 < 1 , the server completes fewer number of class 2 orders by d2. On the other hand, when w1 > 1 , the server completes the wave after d2, which produces lower NSD2. Consequently, NSD2 is maximized when w1 = 1 . Although releasing the wave at w1 = 0:5 maximizes NSD2, all class 1 orders would be missed by d1, which results in NSD1 = 0. What should be the single wave release time to maximize system NSD? Figure 4.7 shows that the system NSD is maximized when w1 1 . Although the system NSD does not change when w1 0:5, each class NSD is a ected by the release time. For example, when w1 = 0:3, NSD1 < NSD2 and when w1 = 0:2, NSD1 > NSD2. Consequently, the optimal release time depends on the importance of high NSD for each class. If the classes are equally important, then the optimal single wave release time w1 = 0:25. 50 Combining equations 4.5{4.11, we de ne the system NSD as a piecewise function of release time w1: NSD = 8> >>>> >>> >>>> >< >>>> >>> >>>> >>: (1 +w1) 1d1 + 2d2 1 + 2 ; if f1 d1 1 1 + 2 1 + 2 1 (d1 w1) + 2 1 + 2 1 (d 2 w1) ; if w 1 d1 d1 2 1 + 2 1 + 2 1 (d2 w1); if f1 >d2; w1 >d1 (4.12) Equation 4.12 is composed of continuous linear pieces forming a polygonal curve. Each piece forms a linear segment, and the system NSD is maximized at one of the extreme points of the segments. As a consequence, we can decompose the problem into distinct linear programming (LP) sub-problems. The maximum system NSD can be found by solving each sub-problem and then choosing the w1 that maximizes the system NSD. When there are more than two classes of orders, the number of subproblems increases. Fortunately, the structure of the problem allows us to generalize results of a single class system (Section 2.2) and develop the system NSD function for single wave systems with more classes. In a system with K order classes, for example, when the server completes the wave before the deadline di, only orders that arrive after the release do not count toward its NSD, and NSDi = 1 (di w1): When the server completes the single wave after the deadline di, we know that NSDi = (di w1)( i=L) i(di w1) i = i ( =L)(d i w1) (di w1) i = (di w1) 1 1 : 51 If the server starts processing the wave after di, non of the class i orders count toward its NSD, thus NSDi = 0. As a consequence, NSDi = 8> >>> < >>> >: 1 (di w1); if the wave completes before di (di w1) 1 1 ; if the wave starts before di and completes after di 0; otherwise and the system NSD for a single wave system is equal to PKi=1 iNSDi=L. 4.3 Multiple-Class, Multiple-Wave Systems When there is a single wave and multiple classes, the wave releases orders from all classes. Therefore a single wave is a mixture of all classes and each class comprises some proportion of the wave based on its arrival rate. In practice, distribution centers typically have multiple customer classes whose orders must be processed together (to achieve economies of scale in operations) within multiple waves. The server may process single or many classes in a wave. That is, an individual wave may be dedicated to a single class or to a composition of multiple classes. There are two limiting situations: (1) Each wave is dedicated to a single class |a class exclusive wave policy, in which orders of a single class comprise the total workload of a wave. (2) Compositions of waves are based purely on the arrival rates of classes | a pure mix wave policy. The number of orders from class i in a wave is equal to the workload of the wave times the proportion of total arrivals of class i to the total daily load. These two policies can be classi ed under a more general policy | a mixture policy, in which compositions of classes in waves are decision variables. We are interested in optimal mixture policies. This problem involves multiple order classes, each having a di erent deadline, and multiple waves per day. Because wave releases are common across order classes, they must be established to maximize multiple objectives. Further, each wave must have an optimal composition of classes which maximizes the system NSD. What should be the timing and the content of 52 multiple waves to maximize the system NSD? This problem is considerably more complex because we must address not only the timing of waves, but also which available orders to release in each wave. 4.3.1 A Two-Class, Two-Wave Order Release System We start by de ning a simple version of the multiple-class multiple-wave system in which there are two classes and two waves (K = 2;N = 2). Arrival rates and deadlines of the classes are 1; 2 and d1;d2 such that d1 d2. Without loss of generality, we assume d2 = 1. Two quantities are of interest: (1) the unworked inventory from the previous cycle, and (2) the orders that arrive within the current cycle. Recall that the cycle of a class is determined by the length of a unit time between two consecutive deadlines. Class 1 orders that arrive after the last release and before d1 in the previous cycle form the unworked inventory for this class. (Similarly, class 2 orders that arrive after the last release and before d2 form class 2?s unworked inventory.) Figure 4.8 illustrates the unworked inventory from previous cycles and the total number of orders that arrive in the current cycle before the rst release of the day. In Figure 4.8, the amount of unworked class 1 inventory depends on the release time of the last wave. Because w2 is the last wave before d1, the unworked inventory for class 1 is equal to 1(d1 w2). For class 2 orders, the unworked inventory by d2 is equal to 2(d2 w2). The total number of orders that arrive within the current cycle before the rst release is a function of the rst release time w1. Together with the unworked inventory, orders that arrive within the current cycle de ne the work content of a wave. Although releasing all unworked inventory rst is a guarantee for satisfying absolute business requirements (e.g., if an order misses a deadline, the DC should ship the order by the second deadline), because workers are indi erent to which order they process, there may be orders which will not be completed even after a couple or more deadlines. (We discuss 53 d1 d2 d1 d2w2 w2w1 Class 1 and 2 orders from their previous cycles Class 1 and 2 orders arrived in the current day Current cycle for class 1 orders Previous cycle for class 1 orders Current cycle for class 2 orders Previous cycle for class 1 orders Day 0 Day 1 Day 2 Figure 4.8: A wave is comprised of orders from their previous cycles (unworked inventory) and orders that arrive in the current day. the implications of di erent sequencing rules in Section 4.3.2). Consequently, the unworked inventory can be distributed to the waves to ship more orders that arrive in their current cycle. De ne yij as the fraction of unworked inventory of class i orders released in wave wj. If w2 d1, then the amount of unworked class 1 inventory released in wave j is 1(d1 w2)y1j; otherwise, 1(d1 w1)y1j. We write similar expressions for class 2 unworked inventory: unworked inventory for class 2 in the second wave is equal to 2(d2 w2)y21 and 2(d2 w2)y22. The number of class 1 orders that arrive in their current cycle is equal to 1(1 d1+w1). The number of class 2 orders that arrive in their current cycle is analogous: 2(1 d2+w1) = 2w1. We are now ready to write the necessary conditions for a two-class, two-wave system. Let xij be the number of class i orders that arrive in their current cycle and are released in 54 wave j (0 xij). Because the total number of class i orders is i(1 di +w1), we have the following condition for the rst wave: xi1 i(1 di + w1); i = 1;2. For the second wave, we have xi2 i(w2 w1); i = 1;2. The second wave must start after the server completes the total work in the rst wave, so w2 f1. Because we do not know which wave is the last one before a speci c deadline, we denote the release time of a wave before d1 by w . Then, completion time of the rst wave f1 = 8 >< >: w1 + x11 + 1(d1 w2)y11 +x21 + 2(d2 w2)y21 ; if w = w2 w1 + x11 + 1(d1 w1)y11 +x21 + 2(d2 w2)y21 ; otherwise. We allow the second wave to be completed after d2. We also require the condition to release all unworked inventory in waves, or P2j=1yij = 1;i = 1;2. Finally, wave release times should satisfy 0 wj 1; j = 1;2. The necessary conditions for a feasible solution are relatively easier to de ne, but the objective function is more challenging. By de nition, if the rst wave is completed before d1, then all orders of class 1 that arrive in the current cycle count toward its NSD. This quantity is, then, equal to x11. If f1 >d1, then only some proportion of class 1 orders will be completed by the deadline. The rst wave includes class 1 orders both from the previous cycle and the current cycle: x11 and 1y11(d1 w1), respectively. Whenf1 >d1, the number of class 1 orders completed by d1 is equal to (d1 w1) 1=( 1+ 2). Assuming that the server randomly selects which order to process (i.e. the server is indi erent to orders of unworked inventory and orders that arrive in the current cycle), the number of class 1 orders that arrive in the current cycle and are completed by d1 is equal to x11 (d1 w1) 1=( 1 + 2) = x11 + 1y11(d1 w1) . For convenience, we denote the number of class i orders that account for their NSD in wave j by nij. Simplifying terms, the number of class 1 orders completed by d1 when f1 >d1 is n11 = x11(d1 w1) [x 11 +y11(d1 w1)] : 55 Similarly, if f2 d1, the number of class 1 orders completed in wave 2 is x12; otherwise (w = w2) n12 = x12(d1 w2) [x 12 +y12(d1 w2)] : Let us now consider class 2 orders. If the rst wave completes before the second deadline (f1 d2, but not necessarily w1 d1), then all class 2 orders count toward its NSD which is equal to x21; otherwise n21 = x21(d2 w1) [x 21 +y21(d2 w1)] many class 2 orders count toward its NSD. Finally suppose that the server completes the second wave by d2. Because all class 2 orders that arrive in the current day and are ready by their deadline, x22 many of them count toward its NSD. On the other hand, if f2 > d2, then n22 = x22(d2 w2) [x 22 +y22(d2 w2)] many class 2 orders count toward its NSD. The completion time of wave 2 is analogous to f1: f2 = 8 >< >: w2 + x12 + 1(d1 w2)y12 +x22 + 2(d2 w1)y22 ; if w = w2 w2 + x12 + 1(d1 w1)y12 +x22 + 2(d2 w2)y22 ; otherwise. Then, we have the following expression for the number of class i orders that are processed in wave j which count toward the system NSD: nij = 8 >< >: xij; if fj di xij(di w ) [xij +yij(di w )]; otherwise (4.13) for i = 1;2 and j = 1;2. Note that variables nij are conditional expressions and we need to introduce new binary variables to the objective function in order to handle these expressions. Finally, we require the last wave to be completed before the release of the next day?s rst 56 wave. That is f2 1 +w1. Consequently, we formulate the problem as follows: max NSD = 2P i=1 2P j=1 nij 2P i=1 i s:t: x11 1(1 d1 +w1) x21 2(1 d2 +w1) x12 1(w2 w1) x22 2(w2 w1) y11 +y12 = 1 y21 +y22 = 1 f1 w2 f2 1 +w1 w2 1 0 wj j = 1;2 0 xij i = 1;2;j = 1;2 0 yij i = 1;2;j = 1;2 Example. Suppose that there are two classes of orders with stationary arrival rates 1 = 2. Assume the last deadline is at d2 = 1. For di erent deadlines of d1 and utilization levels of , what is the optimal system NSD for these systems? The variables of interest are fw1;w2;x11;x12;x21;x22;y11;y12;y21;y22g which can be de- termined by solving the above problem. Using Mathematica?s NMaximize function, we solve for the variables that maximize the system NSD for a two-class, two-wave system. A detailed discussion on the optimization methodology is given in Section 4.3.2. Each point in Figure 4.9 illustrates a solution for the system NSD with given d1 and values. 57 SolidCircle SolidCircle SolidCircle SolidCircle SolidCircle SolidCircle SolidCircle SolidCircle SolidCircle SolidCircle SolidSquare SolidSquare SolidSquare SolidSquare SolidSquare SolidSquare SolidSquare SolidSquare SolidSquare SolidSquare MedSolidDiamond MedSolidDiamond MedSolidDiamond MedSolidDiamond MedSolidDiamond MedSolidDiamond MedSolidDiamond MedSolidDiamond MedSolidDiamond MedSolidDiamond ? = 0.5 ? = 0.95 ? = 0.75 0.2 0.4 0.6 0.8 1.0 d1 0.2 0.4 0.6 0.8 NSD Figure 4.9: Solutions for a two class, two wave system ( 1 = 2). As expected, system NSD improves as the utilization decreases. When deadline d1 is early, the server cannot process many class 1 orders on time, and system NSD is low. This observation suggests that Observation 1. Setting a common deadline for all customer classes improves system NSD. What would be the potential improvement in system NSD if the DC were able to change the deadlines? For example, if class 1 orders have a deadline at noon (d1 = 0:5) and the DC is able to set the deadlines, the system NSD will be improved by 16.6% and 7.2% when there is a common deadline (for = 0:5 and = 0:75). Although improvements are signi cant, Observation 2. There is a little marginal bene t to having a common deadline when uti- lization is high. This can be seen especially when = 0:95. Because the system is heavily loaded, the server is able to complete fewer class 1 and class 2 orders by d1 and d2. When d1 is early, releasing class 2 orders in the rst wave results in fewer completed class 1 orders by d1, which produces a lower system NSD. Consequently, the rst wave must be class exclusive. What 58 should be the content of the second wave? Should it be class exclusive, or should it include some class 1 orders? Denote the percentage of class i orders in wave j by cij. Then, cij = xij + i(di w )yij2P i=1 [xij + i(di w )yij] : (4.14) Figures 4.10{4.12 illustrate the optimal timing of waves and mixtures of classes for di erent values of . The horizontal axis indicates di erent values of d1. Green and blue bars show the proportions of class 1 and class 2 orders, respectively; and dark (and light) bars indicate the proportions in wave 1 (and wave 2). The rst pairs of bars in the gures show the solutions when d1 = 0:1. When = 0:5, the solution allocates only class 1 orders to the rst wave and only class 2 orders to the second wave (a class exclusive wave policy). However, as utilization increases, the solution also allocates class 1 orders to the second wave (the rst pairs of bars in Figures 4.11 and 4.12). Why does the solution allocate more class 1 orders to the second wave? An increased utilization implies a later wave completion, consequently the solution allocates more class 1 orders to the second wave in order to complete more class 1 orders in the rst wave by d1. For = 0:5, the solution produces the same release times and the contents of waves for all d1 0:3; however, the resulting system NSD improves as the rst deadline is set to a later time. When d1 > 0:3, the number of class 2 orders in the rst wave and the number of class 1 orders in the second wave start to increase|suggesting more homogeneous waves. We also observe that the solution produces more homogeneous waves as the deadlines get closer for higher values of utilization. When d1 = d2, the solution proposes a pure mix wave policy for = 0:5. Consequently, Observation 3. When d1 is early, the solution suggests class exclusive waves and as d1 approaches d2, the solution suggests more homogeneous waves. 59 0 1 (%) d1 100 75 50 25 0:1 0:2 0:3 0:4 0:5 0:6 0:7 0:8 0:9 w2 = 0:83 w1 = 0:5 w2 = 0:82 w1 = 0:5 w2 = 0:61 w1 = 0:31 w2 = 0:45 w1 = 0:2 w2 = 0:33 w1 = 0:1 w2 = 0:2 w1 = 0 w2 = 0:2 w1 = 0 w2 = 0:2 w1 = 0 w2 = 0:2 w1 = 0 w2 = 0:1 w1 = 0 Figure 4.10: Timing of waves and class mixtures when = 0:50. 0 1 (%) d1 100 75 50 25 0:1 0:2 0:3 0:4 0:5 0:6 0:7 0:8 0:9 w2 = 0:7 w1 = 0:27 w2 = 0:64 w1 = 0:25 w2 = 0:66 w1 = 0:26 w2 = 0:69 w1 = 0:29 w2 = 0:6 w1 = 0:28 w2 = 0:5 w1 = 0:18 w2 = 0:4 w1 = 0:1 w2 = 0:3 w1 = 0:02 w2 = 0:2 w1 = 0 w2 = 0:1 w1 = 0 Figure 4.11: Timing of waves and class mixtures when = 0:75. 0 1 (%) d1 100 75 50 25 0:1 0:2 0:3 0:4 0:5 0:6 0:7 0:8 0:9 w2 = 0:56 w1 = 0:12 w2 = 0:79 w1 = 0:38 w2 = 0:76 w1 = 0:3 w2 = 0:63 w1 = 0:21 w2 = 0:6 w1 = 0:14 w2 = 0:46 w1 = 0:1 w2 = 0:4 w1 = 0:06 w2 = 0:3 w1 = 0 w2 = 0:2 w1 = 0 w2 = 0:07 w1 = 0 Figure 4.12: Timing of waves and class mixtures when = 0:95. 60 For other values of , the waves have nearly the same number of class 1 and class 2 orders when d1 = d2. This is because the rst wave contains only orders that arrive in the current cycle and all unworked inventory from the previous cycles is released in the second wave (i.e., the solution produced y11 = y21 = 0). As a result, the rst wave release time is later than the one in single class systems (e.g. the release time of the rst wave in a single class system is equal to 0.25 when = 0:75, but the solution found the rst release time as w1 = 0:27). Setting Cuto Times We have one remaining question corresponding to the service level: How many hours before the deadline(s) should the cuto time(s) be set in order to guarantee all customer orders are shipped by their deadlines? In single class systems with deterministic arrival and processing rates, we rst assumed the cuto time is equal to the deadline. Analytical solutions produce maximum NSD, and by setting cuto time equal to the last wave release time, we show that all customers who place their orders before the cuto time receive their orders by the next deadline. That is, service level (which is analogous Type 1 service in inventory theory) is 100% when cuto time is equal to the last release time. Because NSD counts the number of orders that arrived in the current cycle and shipped by its deadline, NSD is also 100% when cuto time is set optimally. Therefore, maximizing NSD is the same as achieving the latest possible cuto time for single class systems. Consider multiple-class systems. Denote cuto time of class i by ci. Recall that NSDi is de ned as the fraction of class i orders that arrived between successive cis and completed before di . To this point, we have assumed that ci = di. Numerical solution of the analytical model determines wave release times and their contents which maximize NSDi. Does setting cuto times equal to the last wave release times produce 100% service level? 61 Suppose there are two-class, two-wave release system in which 1 = 2;d1 = 0:8, d2 = 1. When = 0:75%, release times for this system are w1 = 0:26, w2 = 0:66, which are completed at f1 = 0:66, and f2 = 1. The resulting system NSD is equal to 66%. If cuto times are equal to deadlines (c1 = d1;c2 = d2), then service level is also equal to 66%. If there is a common cuto time for both classes which is set to w2, then all class 2 customers receive their orders by d2, thus service level for class 2 customers is equal to 100%. If class 1 orders would have the same cuto time, then all class 1 orders in the rst wave would be completed by d1; however, too few class 2 orders in the second wave would be completed by their deadline. That is, class 1 customers would have a service level which is less than 100%. If the cuto time were equal to w1, then all class 1 and class 2 customers who place their orders before c = w1 would receive their orders on time; however, because the cuto time for class 2 customers is early, they must wait longer for their orders. Now, suppose there is a unique cuto time for each class. When cuto times are equal to the last wave release times completed before the deadlines, all class 1 and class 2 customers receive their orders by d1 and d2, therefore service level is equal to 100%. Then, for our example, c1 = 0:26, and c2 = 0:66. Table 4.1 shows the sets of cuto times that guarantee 100% service level for di erent values of d1. Table 4.1: Cuto times for di erent deadlines d1 in a two class, two wave system ( 1 = 2). = 0:50 = 0:75 = 0:95 d1 c1 c2 c1 c2 c1 c2 0:1 0 0.10 0 0.10 0 0.07 0:2 0 0.20 0 0.20 0 0.20 0:3 0 0.20 0.02 0.30 0 0.30 0:4 0 0.20 0.10 0.40 0.06 0.40 0:5 0 0.20 0.18 0.50 0.14 0.14 0:6 0.10 0.33 0.26 0.60 0.14 0.14 0:7 0.20 0.45 0.29 0.29 0.21 0.21 0:8 0.31 0.61 0.25 0.66 0.3 0.3 0:9 0.50 0.825 0.25 0.25 0.38 0.38 1 0.83 0.83 0.70 0.70 0.56 0.56 62 Because solutions of the analytical model determine wave release times (and their con- tents) without regard to the cuto times (i.e. cuto times are not decision variables in the model), they do not produce optimal cuto times. That is, maximizing system NSD is not same as achieving the latest possible cuto times when there are multiple order classes. Cuto times in Table 4.1 can be used to establish lower bounds for optimal cuto times; however, a model that produces optimal cuto times is signi cantly more complicated and will be the subject of future research. 4.3.2 Systems with More Classes and Waves In this section, we extend two-class, two-wave systems to more general systems in which there are K di erent types of customers whose orders are released in N waves. Because the number of waves in a typical application is usually less than six and we are interested in systems with a reasonable number of deadlines (i.e. transportation modes), we focus on systems with up to six classes and six waves. Our objective is to nd a vector of variablesfw1;w2;:::;wN;c11;:::;cij;:::;cKNgthat maximizes system NSD. Recall that the number of class i orders in wave j that count toward the system NSD is de ned by nij in Equation 4.13. Because variables nij are conditional expressions, we introduce a new binary variable zij: zij = 8 >< >: 1; if fj di 0; otherwise for alli2f1;:::;Kgandj2f1;:::;Ng. That is, variableszij indicate if wavej is completed before or after deadline i. Incorporating zij into Equation 4.13, the system NSD function 63 may be written as follows: NSD = KP i=1 NP j=1 xij zij + (1 zij)(di w ) (x ij +yij(di w )) KP i=1 i : (4.15) Generalizing our ndings in a two class, two wave system, we determine the following necessary conditions. First, the number of class i orders that arrive in its current cycle and released in wave j should not be greater than the number of orders accumulated until the rst wave release time w1. That is, the solution decides how many class i orders to release and therefore the wave is selective. This condition should also hold for following waves. Consequently, we have xi1 i(dK di +w1) i2f1;:::;Kg; (4.16) xij 2(wj wj 1) 8i2f1;:::;Kg; 8j2f2;:::;Ng: (4.17) A second set of constraints is required to allocate unworked inventories (from the previous day) to the waves. De ning continuous variables yij, we have NX j=1 yij = 1 8i2f1;:::;Kg: (4.18) In practice, Equation 4.18 implies absolute business requirements such as shipment guar- antee by the \second scheduled deadline;" however, because pickers pick orders in a mixed sequence, all orders may still not be ready by their deadline. (We discuss the implications of order sequencing below.) Recall that the wave contents are determined by the unworked inventory from the previ- ous cycles and the number of orders that arrive in the current cycle. For all j2f2;:::;Ng, 64 the completion times of wave j fj = wj + KP i=1 x ij + i(wj wj 1)yij : The expression for f1 is more complicated. As we showed in a two-class, two-wave system, the rst wave contains yi1 of class i orders from its previous cycle. We do not know which is the last wave before di, so we denote the last wave release time before a speci c deadline by w . The total number of orders released in the rst wave is then equal toPKi=1 xi1+ i(di w )yi1 and the rst wave?s completion time f1 = w1 + KP i=1 x i1 + i(di w )yi1 : The optimal solution requires that each wave must be released after the completion of the previous one (non-overlapping waves). Therefore, fj wj+1 8j2f1;:::;N 1g; (4.19) fN 1 +w1 j = N; (4.20) We allow the server to complete the last wave after the last deadline. (Recall that we have assumed that the last deadline dK = 1, which de nes the current day.) The last wave of the current day should then be released before 1. Then, wN 1: (4.21) Together with the non-negativity of variables, necessary conditions 4.16{4.21 lead to the following optimization problem for a K class, N wave system: 65 maximize NSD = KP i=1 NP j=1 xij zij + (1 zij)(di w ) (x ij +yij(di w )) KP i=1 i subject to xi1 i(dK di +w1) 8i2f1;:::;Kg; j = 1 xij i(wj wj 1) 8i2f1;:::;Kg; 8j2f2;:::;Ng NP j=1 yij = 1 8i2f1;:::;Kg fj wj+1 8j2f1;:::;N 1g fN 1 +w1 j = N wN 1 zij2f0;1g 8i2f1;:::;Kg; 8j2f1;:::;Ng 0 wj j2f1;:::;Ng 0 xij i2f1;:::;Kg;j2f1;:::;Ng 0 yij i2f1;:::;Kg;j2f1;:::;Ng The model has linear and nonlinear constraints and a nonlinear objective function including binary variables and non-convexities in the continuous variables (we discuss this issue below); therefore it falls into a class of nonlinear optimization problems with binary variables. Although writing expressions for the constraints is relatively easy, because we do not know which wave is the last wave before a certain deadline, it is more challenging to write the objective function. We develop an algorithm in Visual Basic for Applications that evaluates all possible cases and generates the whole problem for a given number of classes and waves. Because nonlinearities imply many local maxima, solving problems optimally is di cult (if not impossible). Consequently, we choose a meta-heuristic solution algorithm to search for the optimal values of variables that maximize the system NSD. We present the solution approach and validation of the model below. 66 Optimization Methodology Before presenting the proposed optimization method, it is useful to explore the structure of K class, N wave problems. These problems have multiple variables which form a polytope in 3(NK + 1) dimensions. In addition to nonlinearities of variables, due to dependency between variables, we expect the objective function to be nonlinear and non-convex. Note that the objective function surface is not known in advance; however we can characterize it by generating many random solutions. For example, Figure 4.13 illustrates the surface of the objective function in a two class, two wave system in which 1 = 2 and d1 = 0:5;d2 = 1. 0.0 0.1 0.2 0.3w1 0.3 0.4 0.5 w2 0.2 0.4 0.6 System NSD Figure 4.13: Admissible pairs of (w1;w2) indicate that the objective function surface is non-convex. To expose the structure of the objective function, we generate more than 2,300 admis- sible pairs of release times (w1;w2) and plot them with the corresponding objective value. For this problem, we observe that the objective function is non-convex. Although this obser- vation does not imply that any other problem has the same non-convexity, due to an added complexity of the mixture of classes in waves, we expect the objective functions of more complicated problems to be both nonlinear and non-convex. 67 We use a meta-heuristic called di erential evolution (DE), which was rst introduced by Storn and Price (1995) to optimize problems with a real-valued, non-di erentiable, multi- modal objective functions with nonlinear constraints. Because the method searches on a search space, it does not guarantee optimality; however, our test results showed that di er- ential evolution is an appropriate method for multiple-class, multiple-wave problems. (We discuss the details of the method in Appendix B.) Model Validation The analytical models assume a uid model in which discrete arrivals are represented as a continuous stream of work. Consequently, we must validate the solutions before addressing systems with more classes and waves. We represent the ful llment system as a three-stage queueing network, corresponding to the picking, packing, and shipping processes. We assume each stage has twenty servers. For the validation experiment, we simulate systems with two classes and two waves at a utilization level of 0.5. Orders arrive to the system according to a Poisson process with rate 1 = 2 = 1 order per minute and processing times of workers are distributed exponentially with a rate of ve minutes. We choose the relatively simple case of two waves and two classes because it is much easier to interpret the results than it would be for more complex cases. An arrival to the system triggers a number of logical decisions to determine if it will be released in the next wave or not. Suppose that an order arrives in the current day before the rst release time. It is assigned to the queue of the rst wave. If the arrival occurs after the rst and before the second wave release time, it is assigned to the second wave queue. Otherwise, the order waits in the release queue until the rst release in the following day. Because each class has a di erent cycle regulated by its deadline, the amount of unworked inventory must be handled carefully. For class 2 orders that arrive after the last wave, y21 percent of them must wait until the rst release and y22 percent of them must wait until the second release in the following day. For convenience, denote the arrival time of an order in 68 the current day by tA. Delay time of a class 2 order is equal to w1 +d2 tA, if it arrived after w2 and assigned to the rst wave in the following day. (The delay time is w2 +d2 tA if the order is assigned to the second wave in the next day.) When a class 1 order arrives after w2, another decision is required. If tA d1, then it is assigned to the rst wave on the next day; otherwise it is assigned to the next day?s rst or second wave release queues. Once orders are released in waves, they are sequentially processed in picking, packing and shipping. The analytical models assume the server randomly selects which order to process. That is, the server is indi erent to the arrival times of orders. Figure 4.14 demonstrates the realization of this assumption. w1 d1 Current cycle inventory Unworked inventory Figure 4.14: Analytical models assume the server randomly selects orders. In Figure 4.14, the server completes the rst wave after d1 and some unworked class 1 inventory will not be processed by its (second) deadline. In practice, because pickers pick orders both from the previous and the current cycles in the same picking tour, there may be remaining unworked inventory at d1. Further, the amount of unworked and current cycle inventories are only known at the time of release. In simulation, each order arrival changes the number of orders in the release queue and because we do not know which cycle (and wave) the order should be assigned to preserve the class weights, it is not possible to maintain the class weights in waves. Consequently, we use rst-come, rst served (FCFS) and last-come, rst served (LCFS) sequencing policies to set bounds on the expected system NSD. In the FCFS policy, the server processes the orders as they appear in the queue. On the other hand, the orders that arrive in the current cycle are worked rst in the LCFS policy. In 69 both policies, the server has the same processing capacity; however, each rule may result in di erent system NSD. (Because the wave content is divided into two parts, the starting time of the second part can be thought as a sub-wave.) Figure 4.15 illustrates the implications of these two policies. w1 d1 Current cycle inventory Unworked inventory (a) FCFS policy results in lower system NSD. w1 d1 Current cycle inventory Unworked inventory (b) LCFS policy results in higher system NSD. Figure 4.15: Di erent sequencing rules result in di erent simulated average system NSD. When the FCFS policy is used, the server can only start the orders that arrive in the current cycle after completing all unworked inventory. All unworked inventory is completed by d1 (100% of orders are guaranteed for delivery by the second deadline); however, because the server completes the current cycle inventory after d1, there are too many class 1 orders that will not be ready by the deadline. Consequently, the average system NSD will be lower than the expectation. When orders are processed based on the LCFS policy, the server completes more current cycle orders, therefore the average system NSD will be higher than the expectation. (In this case, the DC should not promise a guarantee for the second day delivery.) In Figures 4.14{4.15, we ignored the presence of cuto times. Because cuto times correspond to the wave release times which are completed by the deadlines, when they are set optimally, FCFS and LCFS policies should produce the same service level (i.e., the plots above will be irrelevant). In each simulation run, we gradually increased the rst deadline d1 and kept the second deadline d2 = 1. We insert the output of the analytical model (release times and yij values) as an input to the simulation model. Runs last seven simulated days, with 30 replications. 70 MinusMinus MinusMinus MinusMinus MinusMinus MinusMinus MinusMinus MinusMinus MinusMinus MinusMinus MinusMinus MinusMinus MinusMinus MinusMinus MinusMinus MinusMinus MinusMinus MinusMinus MinusMinus MinusMinus MinusMinus 0.2 0.4 0.6 0.8 1.0 d1 0.2 0.4 0.6 0.8 Simulated average system NSD MinusMinus NSD FCFS ELBracket1NSD RBracket1 MinusMinus NSD LCFS Figure 4.16: Model validation results with simulation when = 0:5. We veri ed the simulation model by rst comparing the average utilization level of each stage with our expectation, and then inserting arbitrary release times to observe how system NSD changes. The average utilization of the workers was observed close to our expectation. As expected, di erent release times resulted in lower average daily system NSD. To validate the model, we compare the simulated average daily system NSD with the results of numerical solutions in each run. In Figure 4.16, the circular points correspond to the system NSD found by the solution of the analytical model for each value of d1. Lower and upper bounds (found by applying FCFS and LCFS rules) are shown with horizontal notches. We observe that the expected system NSD is within its bounds, suggesting that the model valid. Numerical Examples To justify the solution method and explore the underlying characteristics of more com- plex systems, we conduct a numerical analysis. Because the number of order classes and wave releases in most applications is at most six, we limit the experiment to the systems with fewer than six classes and six waves. 71 We rst test the performance of the optimization method by transforming multiple class systems into single class systems. That is, we set all deadlines equal to one while all other parameters remain unchanged. For each problem, we investigated three di erent utilization levels. There are 108 problems total. Closed form solutions for single class, multiple wave systems are available (Section 2.3), and we use Equation 2.6 to assess the performance of the search method for each set of problems. Table 4.2 summarizes the performance of the solution method. Table 4.2: Summary of performance metrics. Utilizationlevel =0:5 =0:75 =0:95 Numof %ofSolns Max Avg %ofSolns Max Avg %ofSolns Max Avg Waves Optimal Gap Gap Optimal Gap Gap Optimal Gap Gap 2 100% 0% 0% 40% 5.76% 3.74% 40% 2.98% 2.74% 3 100% 0% 0% 40% 0.70% 0.62% 40% 1.01% 0.71% 4 80% 1.43% 1.43% 40% 2.09% 0.75% 40% 1.13% 0.67% 5 80% 2.35% 1.52% 60% 1.74% 1.20% 40% 0.35% 0.31% 6 60% 2.15% 1.24% 60% 2.02% 1.67% 60% 0.49% 0.41% In Table 4.2, we solve six problems for a given utilization level and number of classes. Columns labeled with \% of Solns Optimal" represent the percentage of the solutions in which the method found the optimal solution. The method performs better when utilization is equal to 0.5; however, the percentage of problems for which the method produced optimal solutions appears to be only 13.3%. Although the method produced an optimal solution only for a small percentage of the test problems, the average gap between the numerical and the analytical methods is about 1.14%. Example. Suppose each class has the same arrival rate i = i+1;i 2f1;:::;5g, and the total processing capacity is equal to = P6i=1 i= where = 0:5;0:75;0:95. We assume that deadlines are equally distributed through the unit length of time. That is, di = i=6;i2f1;:::;6g. Figure 4.17 shows the numerical results of the analytical models. In Figure 4.17, each data point illustrates a system NSD found by solving the corre- sponding problem with K classes and N waves. The system NSD for each value of is shown 72 2 4 6 Num. of classes 2 4 6 Num. of waves 0.0 0.5 1.0 System NSD NSD nullnull 0.5 NSD nullnull 0.75 NSD nullnull 0.95 Figure 4.17: Numerical solutions of analytical models. with a di erent color. We measure the computational time in seconds. The method solves a problem in less than a minute on average (41.8 seconds); however, the computational time increases as the number of classes and waves increases. (We observe that the computational time does not signi cantly change for di erent values of utilization.) When K = 6;N = 6, the method obtained the solution in more than four minutes. Table 4.3 shows the details of computational times for each problem. Figure 4.17 led us to a number of observations: (1) When utilization increases, it impacts the system NSD. The impact of high utilization on system NSD is most signi cant when there is a single class of orders or when there is a single wave. (2) System NSD improves as the number of waves increases. As in single class systems, we observe that system NSD is improved as the number of waves increases, without regard to the number of classes. The maximum system NSD is observed when there is a single class of orders and the number of waves is six. On the other hand, the system NSD is minimum when there are six classes and 73 Table 4.3: Computational times in seconds. Number of waves Num. of classes 2 3 4 5 6 2 7.9 15.7 19.9 21.2 25.9 0.5 8.0 14.3 18.5 19.8 28.2 0.75 9.9 16.1 17.7 17.3 23.8 0.95 3 11.0 20.0 30.9 57.5 74.1 0.5 11.0 22.9 32.1 64.6 81.1 0.75 12.6 17.1 27.8 56.1 83.7 0.95 4 16.5 35.3 69.2 96.8 124.3 0.5 20.6 26.8 67.4 132.2 123.6 0.75 17.1 30.3 61.9 120.9 93.2 0.95 5 55.6 48.6 76.8 140.2 221.8 0.5 47.4 48.9 74.6 134.5 242.5 0.75 42.1 47.2 79.4 111.2 247.5 0.95 6 48.5 72.8 106.3 143.4 209.8 0.5 45.5 93.1 113.4 184.7 223.5 0.75 43.6 78.4 95.6 169.9 274.5 0.95 a single wave, which suggests that, (3) Systems with multiple deadline should compensate with more wave releases. 4.4 Implications for Practice The results of multiple-class, multiple-wave problems suggest implications for practice. When there are multiple classes, many factors a ect service performance, including the timing of waves and their contents. Increasing the number of waves seem to improve system NSD. The solutions of single class systems suggest that the rst wave should be released at w1 = 1 . Although the server is busy of the time, the rst wave release time should be set to a di erent time in multiple-class systems. This observation implies possible idleness of workers between completion and release times of waves. That is, the workers may go idle in order to improve the system NSD. Managers would have to think carefully how to use this time productively, perhaps with replenishment or other necessary operations. 74 In an order ful llment system that must ful ll orders of multiple classes, deadlines may signi cantly impact service performance. We have assumed deadlines that are evenly distributed throughout the day. What would the system NSD be if the deadlines were close to midnight? What should be the wave release times and contents? Should we expect later cuto times in this case? Earlier deadlines imply an earlier rst wave release, and when deadlines are closer to midnight, waves should be released later. We observe that system NSD improves when deadlines are close together and set as late as possible. Late wave release times and deadlines close to midnight also imply later cuto times. Consequently, DCs may o er late cuto times to their customers while maintaining an acceptable level of system NSD. The numerical results of the analytical models show that a wave should contain more orders of a particular class if that wave is the closest wave to the deadline of that class. When deadlines are closer, the waves should contain the same amount of orders from classes. The sequence of picks in waves also a ects the system NSD. Releasing late orders rst guarantees second day delivery, but it results in lower system performance. When pickers pick the most imminent orders rst, it will produce a higher next day ful llment rate; however, if an order is not delivered by its rst promised deadline, then it may also not be ready for the shipment by the second deadline. Managers should consider the tradeo between absolute business requirements and marketing promises. 75 Chapter 5 Conclusions & Future Research To attract more customers and increase the market share, order ful llment systems o er aggressive service promises. That kind of o er can put pressure on DCs in their ful llment operations. The purpose of this dissertation is to increase customer service by improving or- der releases with optimal timing of waves. In this dissertation, we present the rst systematic investigation of wave planning in deadline-oriented order ful llment systems. Contribution 1. We propose the rst uid model for a single-class, multiple-wave system that operates against a daily deadline. Systems with a single order class have correspondence in industry, especially in internet retailing where DCs promise overnight deliveries. We showed that properly timed waves should be smaller as the deadline approaches, rather than of uniform size, as intuition might suggest. The service performance can be increased by adding more waves, but there is a marginal bene t to increasing the number of waves beyond four or ve. The optimal number of waves is based on both the pick density and the xed time component of waves, and there should be fewer waves when the xed time component is long. The optimal policy provides a more stable workload as well as more tractable ow of work. Maximum service performance can be achieved by setting the cuto time equal to the latest release time. Contribution 2. We determined the optimal number and timing of waves in the presence of workload uncertainty. We addressed daily workload uctuations which may be the result of natural demand variability or the result of internal workforce variability. We discussed how to adjust wave 76 release times and the number of waves to minimize the risk of workload exceeding capacity. In the presence of workload uncertainty, we showed that sacri cing a little from the maximum achievable service level hedges against a drastic drop in performance if utilization is high. While the release time of a single wave is later than what it would otherwise be, multiple waves should be released earlier to maximize expected NSD. We discussed the implications of wave release times on the cuto times. We showed that the cuto time should be no later than the last release time. Setting the cuto time as late as possible is important for improving the service promise. Contribution 3. We showed how to establish the timing and contents of waves in order ful llment systems with multiple order classes. We extended our single class systems to multiple class systems in which customers are grouped into ordering clusters, and each cluster receives service at a di erent frequency. We modeled the problem as a nonlinear optimization problem and used a heuristic method to search for optimal solutions. We justi ed the mathematical models with discrete event sys- tem simulation. We presented numerical examples of multiple-class, multiple-waves systems and discussed the implications of the results. The methods in this dissertation rely on certain assumptions. To determine maximum system NSD, we have assumed that cuto times are equal to deadlines. Although this assumption maximizes system NSD, resulting cuto times are not optimal. Future research should consider cuto times as decision variables and determine optimal cuto times. We disallowed release of a wave if the server is busy. Although this assumption alleviates the need for additional sortation systems and makes the control of ow easier, it requires pickers to wait until all others complete their processing. To increase the productivity of workers, DCs might release waves dynamically, in which case waves are allowed to overlap and there is a requirement for sortation. Future research should address dynamic wave releases so that 77 worker productivity will be as high as possible and ful llment operations are designed to maximize service performance. The validation experiments of multiple-class, multiple-wave systems showed that sys- tem NSD could be improved by implementing a more general policy. Future research should develop more powerful models which consider the tradeo between absolute business require- ments and marketing promises. Future research should also consider uctuating work-force levels in order ful llment systems. An important issue is to determine an optimal work-force level. New models are needed to determine the minimum required (hourly) work-force levels that guarantee an acceptable level of service for systems that release orders in waves. Our research should be of interest to practitioners in at least two ways. First, it helps the workers to clearly articulate the service to be delivered. This is likely to help ful llment systems control the gap between customer expectations and actual performance. Second, it allows DCs to utilize guaranteed delivery as a powerful marketing strategy. Although wave planning is o ered by some commercial WMS providers, they do not determine the optimal number of waves nor the optimal release times. Our research could be used by WMS vendors or consultants for better operational planning. 78 Bibliography Amazon.com (2013). Local Express Delivery Ordering Deadlines. http://www.amazon.com/ gp/help/customer/display.html/?nodeId=201117750. Bates, D. (2012). The Secret to Online Shopping Revealed: Buy on a Tuesday, And Don?t Wait For the January Sales. http://goo.gl/aTo6o. Bernd S., Fabian, W., Dashkovskiy, S., Schonlein, M., Makuschewitz, T., and Kosmykov, M. (2011). Some Remarks on Stability and Robustness of Production Networks Based on Fluid Models. Dynamics in Logistics: 2. International Conference, LDIC, Bremen, Germany, 1:27{35. Borovkov, A. A. (1964). Some Limit Theorems in the Theory of Mass Service. SIAM, 9(4):550{565. Borovkov, A. A. (1965). Some Limit Theorems in the Theory of Mass Service, II Multiple Channels Systems. SIAM, 10(3):375{400. Bozer, Y., Quiroz, M. A., and Sharp, G. P. (1988). An Evaluation of Alternative Control Strategies and Design Issue for Automated Order Accumulation and Sortation Systems. Material Flow, 4:265{282. Bradley, P. (2007). Smoothing The Waves. DC Velocity. Chatterjee, S., Slotnick, S. A., and Sobel, M. J. (2002). Delivery Guarantee and the In- terdependence of Marketing and Operations. Production And Operations Management, 11:393{410. 79 Dai, J. (1995). On Positive Harris Recurrence of Multiclass Queueing Networks: A Uni ed Approach via Fluid Limit Models. Ann. Appl. Probab., 5:49{77. Dai, J. (1999). Stability of Fluid and Stochastic Processing Networks. Dai, J. and Jennings, O. (2003). Stability of General Processing Networks, in: Stochastic Modeling and Optimization. Springer. de Koster, R. (2004). How to Assess a Warehouse Operations in a Single Tour. Technical report, RSM Erasmus University. de Koster, R., Le-Duc, T., and Roodbergen, K. J. (2006). Design and Control of Ware- house Order Picking: a literature review. Technical report, Erasmus Research Institute of Management. de Koster, R., Le-Duc, T., and Roodbergen, K. J. (2007). Design and Control of Warehouse Order Picking: A Literature Review. European Journal of Operational Research, 182:481{ 501. Dieker, A. (2006). Extremes and uid queues. Ph.d. Thesis. Universteit van Amsterdam - Netherlands. Doerr, K. and Gue, K. R. (2011). A Performance Metric and Goal Setting Procedure for Order Ful llment Operations. Forthcoming in Production and Operations Management. Duenyas, I. and Hopp, W. J. (1995). Quoting Customer Lead Times. Management Science, 41(1):43{57. Franco, M. D. (2006). Batch vs. Wave Picking. Operations + Ful llment, pages 57{58. Frazelle, E. H. and Apple, J. M. (1994). Warehouse Operations. McGraw-Hill, NY. 80 Gallien, J. and Weber, T. (2010). To Wave or Not to Wave? Order Release Policies for Ware- houses with an Automated Sorter. Manufacturing Service Oper. Management, 12(4):642{ 662. Gilmore, D. (2006). Warehouse Management: To Wave or not to Wave - Part 2. Supply Chain Digest. Gupta, V., Harchol-Balter, M., Wolf, A. S., and Yechiali, U. (2006). Fundamental Character- istics of Queues with Fluctuating Load. SIGMETRICS Perform. Eval. Rev., 34:203{215. Hu man, J. R. (1988). Order Picking Systems. McGraw-Hill, NY. Internet Retailer (2013). Internet Retailer Top 500 Guide. http://www.internetretailer. com/. Johnson, M. E. and Meller, R. D. (2002). Performance Analysis of Split-Case Sorting Sys- tems. Manufacturing Service Oper. Management, 4:258{274. Kella, O. and Whitt, W. (1996). Stability and Structural Properties of Stochastic Storage Networks. J. Appl. Prob, 33:1169{1180. Kella, O. and Whitt, W. (1998). Linear Stochastic Fluid Networks. J. Appl. Prob, 36:244{ 260. Liu, J. and Lampinen, J. (2005). A fuzzy adaptive di erential evolution algorithm. Soft Computing, 9(6):448{462. Liu, Y. and Whitt, W. (2010). A Fluid Approximation for Large-scale Service Systems. SIGMETRICS Perform. Eval. Rev., 38:27{29. Manjoo, F. (2012). How Amazons Ambitious New Push for Same-day Delivery Will Destroy Local Retail. http://goo.gl/pfXUS. Morris, J. (2008). American Eagle Out tters on Cutting Labor Expenses. Logistics World. 81 Newell, G. (1973). Approximate Stochastic Behavior of n-Server Service Systems with Large n. Springer. Oracle Warehouse Management User?s Guide (2012). Wave Planning. http://goo.gl/ S3hgO. Owyong, M. and Yih, Y. (2006). Picklist Generation Algorithm with Order-Consolidation for Split-Case Module Based Ful llment Centres. International Journal of Production Research, 44-21:4529{4550. Pang, G. and Whitt, W. (2008). Heavy-tra c Limits for Many-Server Queues with Service Interruptions. Perry, D. (2007). Continuous Processing Using a Sorter. Technical report, Vargo Adaptive Software LLC. Perry, O. and Whitt, W. (2011). A Fluid Limit for an Overloaded X Model via An Averaging principle. Technical report, e. Petersen, C. (2000). An Evaluation of Order Picking Policies For Mail Order Companies. Production and Operations Management, 9:319{335. Rao, U. S., Swaminathan, M. J., and J., Z. (2005). Demand and Production Management with Uniform Guaranteed Lead Time. Production And Operations Management, 14:400{ 412. Ridley, A. D., Massey, W., and Fu, M. (2004). Fluid Approximation of a Priority Call Center with Time-varying Arrivals. Simulation Conference, 2003. Proceedings of the 2003 Winter, pages 69{77. SAP Business Solutions (2012). Wave Picks. http://goo.gl/dlkm9. Shang, W. and Liu, L. (2011). Promised Delivery Time and Capacity Games in Time-Based Competition. Management Science, 57(3):599{610. 82 So, K. C. and Song, J.-S. (1998). Price, Delivery Time Guarantees and Capacity Selection. European Journal of Operational Research, 111(1):28 { 49. Speaker, R. (1975). Bulk Order Picking. Industrial Engineering, 7-12:14{18. Storn, R. (1996). On The Usage of Di erential Evolution for Function Optimization. Tech- nical report, North American Fuzzy Information Processing Society. Storn, R. and Price, K. (1995). Di erential Evolution - A Simple and E cient Adaptive Scheme For Global Optimization Over Continuous Spaces. Technical report, International Computer Science Institute. Storn, R. and Price, K. (1997). Di erential evolution{a simple and e cient heuristic for global optimization over continuous spaces. Journal of global optimization, 11(4):341{359. Mathematica Tutorial (2013). Mathematica Tutorial: Numerical Nonlinear Global Opti- mization. http://goo.gl/D8cLJ. Ursem, R. K. and Vadstrup, P. (2003). Parameter Identi cation of Induction Motors using Di erential Evolution. CEC-2003. van der Berg, J. P. (2010). Integral Warehouse Management: The Next Generation in Transparency, Collaboration and Warehouse Management Systems. Management Outlook Publications. Ward, A. R. and Bambos, N. (2003). On the Stability of Queueing Networks with Job Deadlines. Journal of Applied Probability, 40:293{304. Whitt, W. and Liu, Y. (2011). Nearly Periodic Behavior in the Overloaded G/D/s+GI Queue. Stochastic systems, 1:1{71. 83 Appendices 84 Appendix A Feasibility and Stability Conditions A set of wave releases is feasible if each wave can be completed before the next one begins. Let N be the number of waves in a day, and let wi denote the time of the i-th wave release time everyday. The rst wave on day 1 is at w11. The wave is composed of all orders accumulated until that time, which is equal to (w11 0). The rst wave is completed at w11 + (w11 0)= . The completion time of the rst wave should be less than w12 to be feasible. The second wave which releases all orders accumulated within [w11;w12], is feasible if the rst one completes before the second wave. The feasibility condition for wave two is then w11 + (w12 w11)= w12. Similar expressions can be written for the rest of waves on day 1. On day two, the rst wave does not only release orders that arrive that day, but also releases the unworked orders from the rst day. The workload of the rst wave on day 2 should satisfy w21 + (1 w1N +w21)= w21 (the length of a day is represented as the interval of [0,1]). On day m, we have the following expressions for feasibility: wm1 + (1 w m 1 N +w m 1 ) w m 2 ; i = 1; wmi + (w m i w m i 1) w m i+1; 1 0. Otherwise there will be zero work left from the interval [wN 1;wN]. Then work-in-process inventory at end of day 1 is equal to I1 = (1 wN) + maxf (wN wN 1) (1 wN);0g: For the second day, work-in-process inventory is equal to I(2) =I(1) + maxf (wN wN 1) (1 wN);0g = (1 wN) + 2 maxf (wN wN 1) (1 wN);0g : 85 We de ne WIP for day m as: I(m) =Im+1 + maxf (wN wN 1) (1 wN);0g = (1 wN) +m maxf (wN wN 1) (1 wN);0g : Note that, in a steady state limm!1I(m) = (1 wN) only if (wN wN 1) = (1 wN), otherwise the system is feasible but not stable. Now, we establish the stability condition when I(0) 6= 0. Because the rst wave has w1+I(0) orders to process before the rst release of the next day, this quantity should not exceed the available capacity: w1 + w1 +I(0) 1; I(0) w1( + ): 86 Appendix B Comparison of Search Methods and Details of Di erential Evolution We solve multiple class, multiple wave system problems with Mathematica?s built in function NMaximize. The function provides a number of search algorithms including direct search methods such as Nelder-Mead, di erential evolution, simulated annealing, and random search (Mathematica Tutorial, 2013). Although Storn (1996) showed that di erential evolution is an appropriate method for solving real-valued, non-di erentiable, multi-modal objective functions with non-linear con- straints, we rst run preliminary tests to assess the performance of the method. The purpose of the preliminary tests is twofold: selecting the best method for our purpose, and assess the quality of the selected method. We test nine problems to compare the performances of di erent numerical methods: di erential evolution, simulated annealing, Nelder-Mead, and random search. Note that Mathematica?s built-in function NMaximize supports all these methods, therefore we test each one by specifying the method within the function. Based on the problem structure, NMaximize chooses which method to use automatically; however, it also lets users to select a speci c method (and its parameters) for their purposes. We start by running a number of preliminary tests to assess the performance of candidate optimization methods given in NMaximize. (Note that there are other possible heuristic methods such as particle swarm optimization which is not included in the function.) We compare the candidate methods as follows. We rst transformed multiple class problem into a single class system by setting a common deadline (d1 = d2 = ::: = dK) and then comparing the results with analytical solutions. During the preliminary tests, we assume 2{4 waves and test the performance of the methods at three levels of utilization: = 0:5; 0:75; 0:95. We limit computational time with 600 seconds. Because our rst task is to select an appropriate method, we use the default parameter sets of the methods during the preliminary tests. (Method speci c parameters can be found in Mathematica Tutorial, 2013.) The results of the preliminary test are given in Table B.1 in which best solutions are shown with asterisks. Table B.1 suggests that di erential evolution is an appropriate method for the problems we address. Simulated annealing (SA) performed very poorly when compared to the other methods (especially when N = 3 and as utilization increases). Although random search (RS) nds the optimal solutions for all values of utilization when N = 2, this method did not nd any solution within the speci ed time limit for N > 2. Di erential evolution attempts to nd global solution by executing three main routines: mutation, recombination, and selection. The rst step is to generate an initial solution which requires D dimensions to de ne the variables. (For a K class, N wave system, it requires N(3K + 1) dimensions.) In this step, randomly generated points between their bounds constitute an initial feasible solution S. The value of the objective function in the initialization step is denoted by f1. 87 Table B.1: Comparison of objective function values in the preliminary run. Waves Analytical (%) DE (%) SA(%) NM(%) RS(%) 2 83.3 83:3 83:3 83:3 83:3 0.5 3 92.9 92:9 32.2 91.4 - 4 96.7 96:7 96:7 96:7 - 2 67.9 67:9 55.2 67:9 67:9 0.75 3 81.8 81:3 31.6 79.7 - 4 88.4 88:4 87.6 88:4 - 2 53.7 53:6 40.0 52.6 53:7 0.95 3 69.9 69:5 35.3 69:5 - 4 78.1 78:0 77.8 77.8 - After specifying the control parameters (explained below) and generating the initial solution, the algorithm executes its routines to search the solution space. The algorithm iteratively repeats cycles of routines until it is terminated by stopping criteria. By default, these criteria are based on both the convergence of variables and objective function value in the successive iterations. We force the algorithm to terminate if the objective value di erence between two consecutive iterations is less than 10 6 or the computational time is greater than 600 seconds. Because di erential evolution uses sets of solutions in each generation, it requires a population size parameter Ps. In each generation, a population size of Ps is maintained. We denote the ith solution in generation G by SiG. A mutation routine is executed in the next step. Note that the set of solutions SiG;i2Ps contains vectors of variables xiG. For each target vector i, mutation operator rst selects three random points: xuG, xwG, and xzG. Then the operator generates trial vector xTG based on the selected points. This calculation is done by adding the scaled down (by a factor ) di erence of the two vectors to the rst one: xTG = xuG [xwG xzG] (operators \ " and \ " denote vectorial summation and scalar multiplication operators). Because mutation operation replaces all points (as a vector) in the solution by new points in each generation, it expands the solution space. Mutation routine generates the trial vector and with the crossover routine, the algorithm replaces some part of the target vector with this new vector. To crossover, a uniform generated random variable u is compared with the crossover constant pc. If u