School of Electrical Engineering

Marek Matusiak

Aalto-DD 62/2014

Optimizing Warehouse Order Batching when Routing Is Precedence Constrained and Pickers Have
Varying Skills

Aalto University
School of Electrical Engineering
School of Electrical Engineeringwww.aalto.fi

Marek Matusiak

BUSINESS +ECONOMYART +DESIGN +ARCHITECTURE
SCIENCE +TECHNOLOGYCROSSOVER
DOCTORALDISSERTATIONS

Aalto University

9HSTFMG*afgheg+

ISBN 978-952-60-5674-6
ISBN 978-952-60-5675-3 (pdf)ISSN-L 1799-4934
ISSN 1799-4934
ISSN 1799-4942 (pdf)

Optimizing Warehouse
Order Batching when
Routing Is Precedence
Constrained and Pickers
Have Varying Skills

DOCTORALDISSERTATIONSAalto University publication series
DOCTORAL DISSERTATIONS 62/2014

Optimizing Warehouse Order Batchingwhen Routing Is Precedence
Constrained and Pickers Have Varying
SkillsMarek Matusiak

A doctoral dissertation completed for the degree of Doctor of
Science (Technology) to be defended, with the permission of the
Aalto University School of Electrical Engineering, at a public
examination held at the lecture hall AS1 of the school on the 6th ofJune 2014 at 12.

Aalto University
School of Electrical Engineering
Department of Electrical Engineering and Automation
Generic Intelligent MachinesSupervising professor
Professor Emeritus Aarne Halme
Thesis advisors
D.Sc. (Tech) Jari Saarinen
Professor René de Koster
Preliminary examiners
Professor Wout Dullaert, Universiteit Antwerpen, Belgium
Assistant Professor Wilco van den Heuvel, Erasmus University
Rotterdam, the Netherlands
Opponent
Assistant Professor Adriana Gabor, Erasmus University Rotterdam,
the Netherlands

Aalto University publication series
DOCTORAL DISSERTATIONS 62/2014© Marek Matusiak
ISBN 978-952-60-5674-6
ISBN 978-952-60-5675-3 (pdf)ISSN-L 1799-4934
ISSN 1799-4934 (printed)
ISSN 1799-4942 (pdf)
http://urn.fi/URN:ISBN:978-952-60-5675-3
Unigrafia Oy
Helsinki 2014
FinlandAbstract
Aalto University, P.O. Box 11000, FI-00076 Aalto www.aalto.fi

AuthorMarek Matusiak
Name of the doctoral dissertation
Optimizing Warehouse Order Batching when Routing Is Precedence
Constrained and Pickers Have Varying Skills
Publisher School of Electrical Engineering
Unit Department of Electrical Engineering and Automation
Series Aalto University publication series DOCTORAL DISSERTATIONS 62/2014
Field of research Automation Technology
Manuscript submitted 28 April 2014
Permission to publish granted (date) 24 April 2014
Monograph

Date of the defence 6 June 2014
Language English

Article dissertation (summary + original articles)

Abstract
Warehouses are an important part of most supply chains. By batching customer orders androuting pickers effectively, warehouses aim to increase the efﬁciency of the order picking
process. The contributions of this thesis are related to two extensions to the optimizing of
picker-to-parts order picking operations in warehouses.
First, precedence constraints are introduced to picker routing, which pose new challenges toorder batching algorithms. This is due to the relative complexity of precedence-constrained
routing when compared to standard methods for routing pickers. The complexity is mitigatedby the use of an effective savings estimate in calculating the properties of large batches, which
reduces computation time signiﬁcantly. This savings estimate is used in a Large Neighborhood
Search algorithm, which outperforms heuristics from literature and compares well to optimal
solutions (1.2% mean error). An A* algorithm is used for precedence-constrained picker
routing. Compared to the performance of a reference warehouse, almost 16% (5000km) in total
travel distance can be saved during a three month period.
Second, forecasting models of order pickers' batch execution times are built with multilevel
modeling using a three-month set of operational data. It is shown that signiﬁcant differencesin picker performance exist. The forecasting models are used in an Adaptive Large
Neighborhood Search algorithm in ﬁnding the right picker for the job. Compared to a state-ofthe-art batching algorithm and to the current practice of assigning work in most warehouses,
9% is saved in total order picking time. These results show that differences among order
pickers should be taken into account when optimizing order picking operations in warehouses.
This is the ﬁrst work where forecasting models are used to predict the performance ofindividual order pickers, and where such models are exploited with a job assigning algorithm.

Keywords warehousing, routing, order picking, metaheuristics, batching, multilevel modeling,
worker modeling, business analytics, combinatorial optimization
ISBN (printed) 978-952-60-5674-6
ISBN (pdf) 978-952-60-5675-3
ISSN-L 1799-4934
ISSN (printed) 1799-4934
ISSN (pdf) 1799-4942
Location of publisher Helsinki
Pages 137

Location of printing Helsinki Year 2014urn http://urn.ﬁ/URN:ISBN:978-952-60-5675-3

Tiivistelmä
Aalto-yliopisto, PL 11000, 00076 Aalto www.aalto.fi

Tekijä
Marek Matusiak
Väitöskirjan nimi
Keräilyvaraston tilausten yhdistämisen optimointi keräilijöiden ominaisuudet huomioiden ja
reitityksen ollessa etusijarajoitettua
Julkaisija Sähkötekniikan korkeakoulu
Yksikkö Sähkötekniikan ja automaation laitos
Sarja Aalto University publication series DOCTORAL DISSERTATIONS 62/2014
Tutkimusala Automaatiotekniikka
Käsikirjoituksen pvm 28.04.2014
Julkaisuluvan myöntämispäivä 24.04.2014
Monografia

Väitöspäivä 06.06.2014
Kieli Englanti

Yhdistelmäväitöskirja (yhteenveto-osa + erillisartikkelit)

Tiivistelmä
Varastot ovat tärkeä osa teollisuuden ja vähittäismyynnin toimitusketjuja. Varastojen
keräilyprosessia on pyritty tehostamaan keruureittien optimoinnilla sekä yhdistämällä
tilauksia isommiksi kokonaisuuksiksi, jotka voidaan kerätä yhden keruureitin varrelta
lähtöpisteeseen välillä palaamatta. Tämän väitöskirjan tieteelliset kontribuutiot laajentavat
tieteen nykyisiä varastojen vähittäiskeräilyprosessin optimoinnin käytäntöjä.
Kerättävien tuotteiden välillä voi olla etusijarajoitteita, ts. jotkut tuotteet tulisi kerätä ennen
toisia. Perinteiset keräilijöiden reitittämiseen käytetyt algoritmit eivät ota tätä huomioon.
Etusijarajoitettu reitittäminen on laskennallisesti verraittain raskasta. Tilausten
yhdistämiseen käytetyt algoritmit kutsuvat reititysalgoritmia aina löytäessään uuden
tilausyhdistelmän. Reitityksen ollessa etusijarajoitettua, jatkuva reitittäminen tekee
kirjallisuudessa esitettyjen algoritmien käytöstä epäkäytännöllistä ratkaisuun kuluvan ajan
vuoksi. Tässä työssä esitetään tapa vähentää laskennallisesti raskasta reititystä käyttäen
säästettyyn matkaan perustuvia estimaatteja. Estimaatteja käytetään tilausten
yhdistämisalgoritmissa varsinaisen reitityksen sijaan. Tämä nopeuttaa ratkaisun löytämistä
laadusta tinkimättä: esitetty algoritmi löytää paremman ratkaisun tilausten
yhdistämisongelmaan nopeammin kuin mikään testatuista verrokkialgoritmeista, varsinkin
jos tilauksia on paljon. Optimaalisiin ratkaisuihin verrattuna, esitetty algoritmi pärjää hyvin,saavuttaen keskimäärin 1,2% epätarkemman tuloksen. Referenssivaraston alkuperäisten
keruureittien kokonaismatkaan verrattuna säästöä kertyi 16%, eli noin 5000 km kolmen
kuukauden ajanjaksona.
Kolmen kuukauden keruudatan pohjalta monitasomallinusmenetelmän avulla luodaan
keruutehtävän kokonaisaikaa ennustavat kerääjäkohtaiset mallit. Ennuste riippuu kerääjästä
ja työtehtävän parametreistä. Työntekijöiden välisten erojen tilastollinen merkittävyys
osoitetaan varianssianalyysia käyttäen. Ennustemalleja käytetään adaptiivisessa
hakualgoritmissa, joka etsii parhaita keräilijä/tilaus/tilausyhdistelmä-kombinaatioita isosta
hakuavaruudesta. Verrattuna tieteen nykytilaan, jossa tilaukset yhdistetään kokonaismatkaa
minimoiden, mutta ei oteta keräilijöiden taitoja huomioon, saavutetaan 9% säästöt
kokonaiskeruuajassa. Nämä tulokset osoittavat, että keräilijöiden taidot tulisi ottaa huomioon
keräilyvaraston toimintojen optimoinnissa.
Avainsanat varastointi, reititys, keräily, metaheuristiikka, tilausten yhdistäminen,
monitasomallintaminen, työntekijöiden mallintaminen, liiketoimintatiedon
hallinta, kombinatorinen optimointi
ISBN (painettu) 978-952-60-5674-6
ISBN (pdf) 978-952-60-5675-3
ISSN-L 1799-4934
ISSN (painettu) 1799-4934
ISSN (pdf) 1799-4942
Julkaisupaikka Helsinki
Sivumäärä 137

Painopaikka Helsinki
Vuosi 2014urn http://urn.ﬁ/URN:ISBN:978-952-60-5675-3

PrefaceThis research was conducted in the Centre of Excellence in Generic Intelligent Machines Research, partly funded by the Finnish Academy of
Sciences, in the Department of Automation and Systems Technology during the years 2010-2014. Funding was in part received from the Motti
project, which is a part of the Energy and Life Cycle Cost Efﬁcient’ Machines (EFFIMA) research program, managed by the Finnish Metals andEngineering Competence Cluster (FIMECC), and funded by the Finnish
Funding Agency for Technology and Innovation (TEKES), research institutes and companies. Their support is gratefully acknowledged.
I would like to thank Professor Aarne Halme for the opportunity to
pursue my Doctorate and for his role as my supervisor. Professor Renéde Koster of Erasmus University has been integral to the content of histhesis, you are a great mentor and an inspiration. My instructor Jari
Saarinen inspired me to pursue this subject: thank you for telling me
to leave the robots alone, and for the many ideas and discussions. The
pre-examiners of this thesis, Professor Wout Dullaert and Assistant Professor Wilco van den Heuvel have my gratitude for taking the time to readthrough the thesis and helping to make it better.
Second, I would like to thank the people who helped with work donein this thesis. My father, Professor Jerzy Matusiak helped me to organize the thesis and supported me the whole way. My roommate at work,
Eemeli Aro, with whom I had numerous discussions on my research;
Janne Paanajärvi for being a fountain of knowledge in all things scientiﬁc;Professor Leo Kroon for modelling the precedence constrained routing andbatching; my boss, Professor Ville Kyrki for his advice on things related tothe thesis; Antti Maula for his help on many coding issues; Johan Grönholm for his work on building the distance matrix for the warehouse; Jani
Mähönen of Rocla for presenting me with these warehousing problemsiPreface(and for the data!); Olli Haavisto for his regression analysis primer; andAssistant Professor Murat Tarakci for the idea of using multilevel analysis to model the pickers.
The research done in this thesis would not have been possible without
my periodic trips to the Erasmus University in Rotterdam to visit Professor René de Koster, funded by a generous and ﬂexible travel grant from
Teknillisen korkeakoulun tukisäätiö. Thank you.
Finally, I would like to thank my family: my wife Salla, who has beenboth amazing and amazingly supportive; my great kids for sleeping betterof late; and my mother for always being there for me.

Espoo, May 7, 2014,

Marek Matusiak

ii

Contents

PrefaceiContents

iii

List of Figuresvii

List of Tables

xi

List of Symbols

xiii

List of Acronyms

xviii

Glossaryxx

1. Introduction

1

1.1 Background and Motivation . . . . . . . . . . . . . . . . . . .

1

1.2 Objectives. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7

1.3 Contributions of the Thesis . . . . . . . . . . . . . . . . . . .

7

1.4 Related Work . . . . . . . . . . . . . . . . . . . . . . . . . . . .

91.4.1 Routing in Picker-to-Parts Warehouses . . . . . . . .

9

1.4.2 Routing with Precedence Constraints . . . . . . . . .

11

1.4.3 Order Batching . . . . . . . . . . . . . . . . . . . . . .

12

1.4.4 Heterogeneous Workers and Jobs . . . . . . . . . . . .

15

1.4.5 Very Large Scale Neighborhood Search . . . . . . . .

16

1.4.6 Multilevel Modeling . . . . . . . . . . . . . . . . . . . .

18

1.4.7 Forecasting Order Picking Time. . . . . . . . . . . .

18

1.5 Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

18

1.6 Author’s Contribution . . . . . . . . . . . . . . . . . . . . . . .

19

1.7 Assumptions . . . . . . . . . . . . . . . . . . . . . . . . . . . .

19

2. Relevant Mathematical Models and Algorithms

21

iii

Contents

2.1 Mathematical Models . . . . . . . . . . . . . . . . . . . . . . .

21

2.1.1 Combinatorial Optimization Problems . . . . . . . . .

21

2.1.2 Routing . . . . . . . . . . . . . . . . . . . . . . . . . . .

22

2.1.3 Order Batching . . . . . . . . . . . . . . . . . . . . . .

23

2.1.4 Generalized Assignment Problem . . . . . . . . . . . .

24

2.2 Algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

25

2.2.1 Dynamic Programming Algorithm for Optimal
Batching . . . . . . . . . . . . . . . . . . . . . . . . . .

25

2.2.2 Clarke and Wright Based Algorithms . . . . . . . . .

26

2.2.3 Variable Neighborhood Search . . . . . . . . . . . . .

27

2.2.4 Attribute Based Hill Climber . . . . . . . . . . . . . .

29

2.2.5 Large Neighborhood Search . . . . . . . . . . . . . . .

31

2.2.6 Adaptive Large Neighborhood Search . . . . . . . . .

31

2.3 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

32

3. Order Batching when Items Have to Be Picked in a Strict
Sequence35

3.1 Problem Description . . . . . . . . . . . . . . . . . . . . . . .

36

3.1.1 Notation . . . . . . . . . . . . . . . . . . . . . . . . . .

37

3.2 Solution Approach . . . . . . . . . . . . . . . . . . . . . . . . .

39

3.2.1 Estimating Savings from Batching More than Two
Customer Orders . . . . . . . . . . . . . . . . . . . . .

40

3.2.2 Pick Route Construction . . . . . . . . . . . . . . . . .

43

3.2.3 Batching . . . . . . . . . . . . . . . . . . . . . . . . . .

46

3.3 Validation and Results . . . . . . . . . . . . . . . . . . . . . .

46

3.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

55

4. Batching when Routing with a Generic TSP Algorithm

59

4.1 Problem Description . . . . . . . . . . . . . . . . . . . . . . .

59

4.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

61

4.3 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

62

5. Forecasting Batch Execution Time for Individual Pickers

65

5.1 Warehouse and Dataset Description . . . . . . . . . . . . . .

66

5.2 Data Cleaning . . . . . . . . . . . . . . . . . . . . . . . . . . .

66

5.3 Multilevel Modeling . . . . . . . . . . . . . . . . . . . . . . . .

68

5.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

77

6. Assigning Orders and Batches to Pickers of Varying Skill79

6.1 Joint Batching and Generalized Assignment Problem . . . .

79

iv

Contents

6.2 Note on Problem Complexity . . . . . . . . . . . . . . . . . . .

82

6.3 Solving the Joint Batching and Generalized Assignment
Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

82

6.3.1 Routing Heuristic . . . . . . . . . . . . . . . . . . . . .

82

6.3.2 Initial Assignment of Batches . . . . . . . . . . . . . .

83

6.3.3 Assignment of Customer Orders Based on Average
Picker Productivity . . . . . . . . . . . . . . . . . . . .

84

6.3.4 Adaptive Large Neighborhood Search Algorithm . . .

84

6.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

87

6.4.1 Parameter Calibration . . . . . . . . . . . . . . . . . .

87

6.4.2 Data Preparation . . . . . . . . . . . . . . . . . . . . .

87

6.4.3 Experimental Setup . . . . . . . . . . . . . . . . . . . .

88

6.4.4 Time Savings and Comparison to Original Solution .

89

6.4.5 Impact on Picker Productivity . . . . . . . . . . . . . .

91

6.4.6 Effect of Skill on Batch Assignment . . . . . . . . . .

94

6.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
7. Conclusions

105

Bibliography

107v

Contentsvi

List of Figures

1.1 The Supply Chain Process. Adapted from Beamon (1998). .

2

1.2 An order picker towing two empty roll containers. Image
courtesy of Rocla. . . . . . . . . . . . . . . . . . . . . . . . . .

3

1.3 Top view of a warehouse three cross-aisles (marked with arrows). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

10

1.4 S-shape (left) and largest-gap (right) routing (Roodbergenand De Koster, 2001a). . . . . . . . . . . . . . . . . . . . . . .

10

1.5 Aisle-by-aisle (left) and optimal (right) routing (Roodbergenand De Koster, 2001a). . . . . . . . . . . . . . . . . . . . . . .

11

2.1 VNS neighborhoods and search heuristics. Adapted from
Albareda-Sambola et al. (2009).

. . . . . . . . . . . . . . . .

28

3.1 Top view of the warehouse with seven drop-off locations indicated at the top (marked 1-7). The empty bin depot isdenoted with D. The 57 aisles are unidirectional, while thethree cross-aisles are bidirectional. The aisle-to-aisle distance is 5.5 m the slot-to-slot distance is 3.7 m, which isalso the width of the central cross-aisle. . . . . . . . . . . . .

36

3.2 Solution approach for solving the batching problem with estimated savings with PCES. . . . . . . . . . . . . . . . . . . .

40

3.3 A truck picking three customer orders into three separate
bins. Items in bold have already been picked; the underlined
one denotes the last pick and italicized ones are possibilities
for the next transition, representing state (2,2,1,3). . . . . .

44

vii

List of Figures

3.4 Left: a Hasse diagram of the state space formed by two customer orders with two items to be picked into each of them.
Inside the circles, the state name is at the top, the state vector is at the bottom. The ﬁrst member of the state vector
refers to the customer order where the last pick occurred,
the following two refer to the customer order location indices. The highest location index (3 in this case) represents
the delivery location for the given customer order. Right:
possible transitions and distances from the state m9 in the
Hasse diagram. . . . . . . . . . . . . . . . . . . . . . . . . . .

44

3.5 Error % in comparing batch ranks, ordered by the batchwith the highest total rank for 10,000 batches of three customer orders sampled from a set of 1,536 customer orders.
Only 4,900 are shown here so that small ranks can be distinguished. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

50

3.6 Comparison of PCES to other heuristics and an optimal solution for batch sizes of 3. The top ﬁgure shows the qualityof the solution when compared to the C&W(ii), while thebottom one shows the computation time for each wave sizefor solving an instance of 1,536 orders.

. . . . . . . . . . . .

53

3.7 Comparison of PCES to other heuristics and an optimal solution for batch sizes of 4. The top ﬁgure shows the qualityof the solution when compared to the C&W(i), while the bottom one shows the computation time for each wave size forsolving an instance of 1,536 orders.

. . . . . . . . . . . . . .

54

3.8 Percent of travel distance saved for optimal (solid lines) and
PCES based batching (dashed lines) as a function of the
wave size R in customer orders, for batches of size N = 2,
N = 3, and N = 4, in comparison to N = 1. The mean errorof PCES batching compared to optimal is 1.2% when N = 3.

55

4.1 Comparison of PCES to other heuristics and an optimal solution for batch sizes of 3. The top ﬁgure shows the qualityof the solution when compared to the C&W(i), while the bottom one shows the computational time for each wave size.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.1 Cross-validation residuals for ﬁnal full model vs.

batch

time, grouped per picker. . . . . . . . . . . . . . . . . . . . . .

viii

63

73

List of Figures

5.2 Dendogram of picker groups formed with hierachical clustering . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

75

6.1 Savings generated by different algorithms, compared to theoriginal batching with forecast execution times during a virtual day. 7.15h minimum time to qualify, 29 pickers. The
real batch execution times for the original batching are 3%
smaller than the forecast in this case. . . . . . . . . . . . . .

90

6.2 Savings generated by different algorithms, compared to theoriginal batching with forecast execution times during a virtual day. 7.4h minimum time to qualify, 20 pickers. The
real batch execution times for the original batching are 1%
smaller than the forecast in this case. . . . . . . . . . . . . .

90

6.3 Change in average productivity by picker category. 7.15h
minimum time to qualify, 29 pickers . . . . . . . . . . . . . .

92

6.4 Change in average productivity by picker category. 7.4h
minimum time to qualify, 20 pickers . . . . . . . . . . . . . .

93

6.5 Change in number of lines assigned by picker category.
7.15h minimum time to qualify, 29 pickers. . . . . . . . . . .

93

6.6 Change in number of lines assigned by picker category. 7.4h
minimum time to qualify, 20 pickers. . . . . . . . . . . . . . .

94

6.7 Picker categories after ﬁnal assignment. 7.15h minimum
time to qualify, 29 pickers on average. . . . . . . . . . . . . .

95

6.8 Picker categories after ﬁnal assignment. 7.4h minimum
time to qualify, 20 pickers on average. . . . . . . . . . . . . .

95

6.9 The average number of order lines in a batch assigned toa picker vs. the corresponding ln(Lines) forecasting model
coefﬁcient of that picker, for two different ways of solving
the BatchGAP. . . . . . . . . . . . . . . . . . . . . . . . . . . .

98

6.10 The travel distance to pick a batch assigned to a picker vs.
the corresponding ln(T ravel) forecasting model coefﬁcient ofthat picker, for two different ways of solving the BatchGAP.

99

6.11 The average mass of a batch assigned to a picker vs. thecorresponding ln(M ass) forecasting model coefﬁcient of thatpicker, for two different ways of solving the BatchGAP. . . . 100
6.12 The mean pick level of a batch assigned to a picker vs. thecorresponding ln(Level) forecasting model coefﬁcient of thatpicker, for two different ways of solving the BatchGAP. . . . 101

ix

List of Figures

6.13 The average volume of a batch assigned to a picker vs.
the corresponding V ol forecasting model coefﬁcient of thatpicker, for two different ways of solving the BatchGAP. . . . 102

xList of Tables

3.1 Results from rerouting and combined rerouting and rebatching of the original batches. . . . . . . . . . . . . . . . . .

55

5.1 Summary of data cleaning. The cleaning categories are in
the order they are applied — thus Speed is ﬁrst and Picker
out is last. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

67

5.2 Mean values and coefﬁcients of variation per batch before
and after data cleaning.

. . . . . . . . . . . . . . . . . . . . .

67

5.3 Multilevel model selection procedure . . . . . . . . . . . . . .

71

5.4 Model betas for 5LLLLX with scaled data, effects on ln(T ime). 74
5.5 Picker clusters, cluster sizes and mean coefﬁcients orderedby cluster height in Figure 5.2. . . . . . . . . . . . . . . . . .

76

xi

List of Tablesxii

List of Symbols

SymbolDescriptionarZero-one vector indicating which orders from O are included
in batch rao,rZero-one variable indicating whether order o is included inbatch rbBin (e.g., a roll container) to which a customer order is to bepicked intoci

The cost to pick order ici,jThe cost to pick a batch containing orders i and jcr

The cost to pick all items in batch r

dg,kDistance between items g, k ∈ Gd(mi , V)¯ i , V)d(mDistance to state mi for a set of orders V
Estimated distance from state mi to the goal for a set oforders V

dr

Travel distance for batch re(s)Destroy function for a large neighborhood search algorithm
for solution sf (s)Function which maps a solution s to its costg

An item (order line), contains information to which customer order it belongs to, the pick location and the numberof units to be pickedh(s)Neighborhood search heuristic (a move), which enables
searching the neighborhood of s for new solutionsi(s)Repair function for a large neighborhood search algorithm
for partial solution s

lr

Number of order lines (items) in batch r

mrState vector of the A∗ -algorithm in Section 3.2.2

oCustomer orderxiii

List of Symbolso(g)Function indicating the customer order to which item g belongs to

pg

Sequence number of item gqhProbability of choosing heuristic h with the roulette wheel
methodr

Batch (a combination of multiple customer orders)r(b)Function returning the batch to which bin b belongs tosSolution to a problems∗

Optimal solution to a problem

tr,w

Batch execution time for picker w to complete batch r

u0,w

Between group error term in a multilevel model for group(picker) w

v

Vector collecting adjustable parameters for an Adaptive
Large Neighborhood Search algorithmverEstimated batch rank for batch r

vvr

Real batch rank for batch rw

An order picker

xr

Parameter vector for batch r

BDistance matrix

C

Number of possible batch combinations in a solution (or partof it)

Cn+

Number of possible batch combinations in a solution forbatches of size larger than n

D

Depot

E(eκw )The Smearing Estimate (Duan, 1983) for picker w

I

Instance of a combinatorial optimization problem

J

Number of iterations of an algorithm

Lw

Productivity of picker w

Mmax

Maximum working time during a virtual day

Mmin

Minimum working time during a virtual day

Mw

Maximum working time for picker w

NMaximum capacity of a picking device, i.e., the batch size

K0

Starting temperature of a simulated annealing algorithm

P

Maximum number of items in a customer order

Q

Random number with a uniform distribution

R

Order picking wave size, i.e., the number of orders considered when running a batching algorithm

xivList of Symbols

R2

R-squared, the amount of variance explained by a regression model

Rc2

Conditional r-squared, the amount of variance in a multilevel model explained by ﬁxed and random effects

2
Rm

Marginal r-squared, the amount of variance in a multilevel
model explained by ﬁxed effects

Si,j

Savings resulting from batching two customer orders i andj

Smax

Maximum savings

Snew

New savings

Sorg

Original savings

Stot

Total savingsSP

Savings resulting from batching the customer orders in set
P

Sr

Savings of batch r

SS

Sum of squared residuals

S̄i,j,kEstimated savings resulting from batching customer ordersi, j and k

S̄P

Estimated savings resulting from batching the customer orders in set P

T

Temperature in a simulated annealing routine

T0

Initial temperature in a simulated annealing routine

Tn

Temperature in a simulated annealing routine at iteration
no. n

V

Time complexity of an algorithm (big-oh)

UgIndicates the relative position of item g in the tour thatpicks gXr

Zero-one decision variable whether batch r is included in asolution

Xg,k,rZero-one decision variable whether batch r goes directlyfrom item g to item k

Xw,rZero-one decision variable whether batch r assigned topicker w

Yo,b

Zero-one decision variable whether order o is picked to bin bA

Set of all aisles

B

Set of all binsG

Set of all items (order lines)

Gf

Set of ﬁrst items

xv

List of Symbols

Gl

Set of last items

Go

Set of all items included in order oH

Set of all neighborhood search heuristics of an algorithm

N (s)The neighborhood of solution s

Nh (s)The neighborhood of solution s when searched with heuristic h

Lr

Set of drop-off locations for batch r

O

Set of all customer ordersP

Set of batchesR

Set of all possible batchesR+

Set of all batches with positive savingsR+e

Set of all batches with positive estimated savingsR+v

Set of all batches with positive real savingsS

Set of feasible solutions for a problem

V

Subset of customer orders V ∈ O

W

Set of all order pickers

Yr

Set of all possible tours for batch r

Greek
βi,w

Forecasting model coefﬁcient i for picker wδNumber of iterations after which an update occurs

r

Proportional batch rank error for batch r

γ0Common intercept for a multilevel modelκrWithin-group error termλDecay parameterμMean valueφ

Cooling rate in geometric simulated annealing schedule

σ1

Update to the score of a heuristic if a new solution is theglobal bestσ2

Update to the score of a heuristic if a new solution is better
than a current one

σ3

Update to the score of a heuristic if a new solution is accepted

σ2Variance (within group)

τ0

Between-group variation in intercepts

τ1

Between-group variation in slopes

θh

The weight of a heuristic hΨhScore of a heuristicxvi

List of Acronyms

AcronymDescription

A*

An extension of Dijkstra’s algorithm, which uses a costto-goal heuristic to minimize exploration

ABHC

Attribute Based Hillclimber

AGV

Automated Guided Vehicle

ANOVA

Analysis of Variance

ATSP

Asymmetric Traveling Salesman Problem

AIC

Akaike information criterion

ALNS

Adaptive Large Neighborhood Search

AS/RS

Automated storage and retrieval system

CPLEX

Optimization software by IBMCV

Coefﬁcient of variation, i.e., standard deviation divided
by the mean.

C&W

Clarke and Wright savings algorithm

GAP

Generalized Assignment Problem

GLS

Generalized Least Squares

ICC

Intraclass Correlation Coefﬁcient

LKH

Lin-Kernighan-Helsgaun algorithm, capable of solving
TSP and ATSP problems

LNS

Large Neighborhood Search

OLS

Ordinary Least Squares

PCES

Precedence Constrained Estimated Savings based
batching

REMIX

A random destroy-random repair neighborhood searchheuristic

SA

Simulated Annealingxvii

List of Acronyms

SA-REMIXA Large Neighborhood Search algorithm, which uses
the REMIX neighborhood search heuristic and a simulated annealing search scheme

SKU

Stock Keeping Unit

SOMB

Self-organization Map batching algorithm

SOP

Sequential Ordering Problem

TSP

Travelling Salesman Problem

VIF

Variance Inﬂation Factor

VLNS

Variable Large Scale Neighborhood Search

VNS

Variable Neighborhood Search

WMS

Warehouse Management System

xviii

Glossary

batch execution time the total time to complete a batch of orders, starting at the depot having empty roll containers on board the picking
truck, and ending ending at the depot after the delivery of the orders.
batching the act of combining multiple customer orders that are to bepicked simultaneously during one pick tour, with the aim of saving
order picking time during a pick tour.
capability the ability of a worker to do a job (may be binary).
customer order consists of order lines, each line for a unique product orstock keeping unit (SKU), in a speciﬁed quantity (De Koster et al.,2007).
destroy a method for breaking apart some part of a solution.
distribution center an often large warehouse supplying one or multiple
retailers.
ﬁxed effect in multilevel modeling, the average model of all groups.
high-level a level of picking items that is reachable by pickers only withassistance (usually mechanical).
insertion see repair.
low-level a level of picking items that is reachable by pickers without
assistance (usually mechanical).
order short for customer order.

xix

Glossaryparts-to-picker an order picking method, where the picker remains stationary at a pick station and items are transported to him/her via
some, usually automated, method.
picker-to-parts an order picking method, where the picker travels between warehouse aisles to pick items.
precedence constraint a constraint specifying that a certain thing
must be dealt with before some other thing.
random effect in multilevel modeling, the deviation from the ﬁxed effectfor a group.
removal see destroy.
repair a method for reinserting a destroyed part of a solution to form a
new solution.
roulette wheel method a method of randomly choosing a neighborhood
search heuristic from a set of heuristics based on the heuristic’s performance scores.
skill the speed of at which a job can be carried out.
wave a set of orders that are processed together in an algorithm, see also
wave picking.
wave picking the act of picking a certain number of orders for a common
destination or a completion deadline, often started simultaneously
by multiple pickers.
zone picking an order picking method where a picker is responsible forpicks in a certain part, i.e., a "zone", in the warehouse.

xx1. Introduction

1.1

Background and Motivation

Supply chains are networks that are composed of various business entities(suppliers, manufacturers, distributors and retailers) who collaborate inan effort to: (i) acquire raw materials, (ii) convert the acquired materials
to speciﬁed products, and (iii) deliver the products to retailers (Beamon,
1998). Traditionally, materials ﬂow forward and information ﬂows backward in the supply chain.
A supply chain is composed of two basic processes: the Production Planning and Inventory Control Process and the Distribution and Logistic Process, see Figure 1.1 (Beamon, 1998). Production planning describes thedesign and planning of the entire manufacturing process. Inventory control describes design and management of storage policies and procedures
for raw materials and inventories of unﬁnished and ﬁnished products.
The Distribution and Logistic Process speciﬁes how products are transported from the warehouse to retailers. Products can be transported directly to retailers, or to a warehouse supplying the retailer (a distribution
center).
Entities involved in a supply chain can use anticipatory or responsive
business models, or a combination of both (Bowersox, 2011). Anticipatory
business models are based on forecasting demand, producing or acquiringand storing the goods ﬁrst and hoping that the forecasts hold. Responsivemodels sell ﬁrst and manufacture and/or procure goods after the sale, delivering them to the customers as the goods become available. Most retailers use anticipatory business models, as they try to offer goods basedon forecast demand (and sometimes the availability of goods).
Warehouses are an important part of most supply chains. More efﬁcient1

Introduction

Figure 1.1. The Supply Chain Process. Adapted from Beamon (1998).

supply chains have resulted in smaller warehouses being replaced withfewer large warehouses. These distribution centers provide timely andeconomical inventory replenishment for retailers (De Koster et al., 2007,
Bowersox, 2011). As warehouses become larger, companies are more reliant on their efﬁcient functionality: in many cases, a company’s whole
supply chain is dependent on the performance of a single central warehouse. Between the time an order arrives at a warehouse and its eventual
arrival at its destination, there is ample opportunity for errors and room
for improvement (De Koster et al., 2007).
Warehouses can employ humans as pickers and/or be automated to somedegree. When humans are employed, three picking methods can be distinguished: picker-to-parts, parts-to-picker, and put systems (De Koster et al.,2007). In picker-to-parts order picking, the picker drives or walks among
the aisles to pick items. Parts-to-picker systems employ automated storage and retrieval systems (AS/RS) to bring the goods to the picker who
works at the depot or a pick station. Finally, in put systems (order distribution systems), items for multiple orders are picked into a container
either in a picker-to-parts or parts-to-picker manner, offered to an order
picker for distribution over customer orders. Zone picking (Goetschalckx
and Ashayeri, 1989) is a scheme where the picker is responsible for all
picks in some area of the warehouse (such as a number of aisles) andplaces the picked goods on an automated guided vehicle (AGV) or a conveyor belt. Bartholdi et al. (2001) propose Bucket Brigade picking system,a type of work sharing zone picking system, where the pickers are not restricted to a single zone, but may operate in many, with the restrictionthat they may not pass each other. If the pickers are sequenced fromslowest to fastest, a balance of work will spontaneously emerge resulting
in high throughput.

2

Introduction

Figure 1.2. An order picker towing two empty roll containers. Image courtesy of Rocla.

This thesis focuses on optimizing processes in the type of warehousemost common in the world today: the low-level picker-to-parts warehouse,a type which, according to De Koster et al. (2007), forms a very large majority (approximately 80%) of all warehouses worldwide. The popularity ofsuch warehouses can be attributed to the relative ease of setting up such afacility and the low initial capital cost when compared to AS/RS systems.
In low-level picker-to-parts warehouses, the order picker picks requested
items from storage racks or bins while traveling along the storage aisles(see Figure 1.2). Conversely, in high-level picker-to-parts warehouses,
cranes or lifts mounted on order picking trucks are employed to lift the
picker up so that he or she can reach all of the stored items (De Koster
et al., 2007).

3

Introduction

Tasks in a Warehouse

There are many activities that are needed to getmaterials in and out of the warehouse. Listed roughly in the order of
execution from the arrival to the departure of goods, and according to
Tompkins et al. (2003), they are:
1. Receiving is a collection of the following activities: (i) the orderly receipt
of all materials coming to the warehouse; (ii) assuring the quantity andquality of the items; and (iii) disbursing the materials.
2. Inspection and quality control are extensions of the receiving process
and are done, e.g., when suppliers are inconsistent in quality or the
goods must be inspected for some other reason.
3. Repackaging is done when goods are received in bulk and are to be
stored in some other manner, e.g., in smaller units or with other products.
4. Putaway is the act placing materials in the warehouse.
5. Storage is the physical containment of the goods while it is awaiting
demand.
6. Order picking is the process of removing items from storage to meetdemand. It is the basic service around which most warehouse designs
and processes are based.
7. Postponement is an optional step done after the picking process, where
items are repackaged for more convenient use. Instead of doing this in
the repackaging step, performing it after picking has the advantage ofallowing more ﬂexibility in the use of on-hand inventory.
8. Sortation can be done by the picker by picking items directly to theright customer order, or can be done after picking, if items are not pickedby order, but by type, for example.
9. Packing and shipping can include the following tasks: (i) checking order completeness; (ii) packaging merchandise in containers appropriate for shipping; (iii) preparation of shipping documents; (iv) weighing4Introduction

orders to determine shipping charges; (v) grouping orders to await an
outbound carrier; and (vi) loading trucks.
10. Cross docking is the act of moving inbound goods directly to the shipping dock.
11. Replenishing involves moving items from reserve to primary picking
locations.

Order Picking

To remain competitive, warehouses need to continuously

improve the efﬁciency of their processes. According to Gu et al. (2007),
research on improving warehouse operation efﬁciency has focused on fourareas: order picking, receiving, storage, and shipping. Of these four areas,
order picking is the most important, and it is the most expensive process
in distribution centers. De Koster et al. (2007) deﬁne order picking as "theprocess of retrieving products from storage (or buffer areas) in response toa speciﬁc customer request". In the literature, of the total cost of operatinga distribution center, the cost of order picking is estimated to be as highas 55% (Drury, 1988) to 65% (Coyle et al., 1996). In spite of the advent ofhighly mechanized warehouses, orders are still picked manually in mostwarehouses, with pickers traveling in the warehouse to retrieve the items
for an order (De Koster et al., 2007). These manual picking processes are
increasingly supported by advanced technologies, such as pick-by-lightsystems, supporting mobile terminals, and voice-picking systems, whichhave made the order picking process more reliable and efﬁcient (Berger
and Ludwig, 2007, Weaver et al., 2010, Baumann, 2013). These computerassisted picking processes generate extensive data logs, including who
performed what pick at which location and at what point in time.
Travel time is the largest component of the total time of the order picking process (batch execution time), with a contribution of up to 50% (Tompkins et al., 2003). Much research has been devoted to methods reducing
order pickers’ travel times, including modifying the warehouse layout, using a different storage policy, batching orders, and routing pickers. For
this thesis, the most relevant literature is related to the batching of customer orders. Customer orders consist of order lines, each line for a uniqueproduct or stock-keeping unit (SKU), in a speciﬁed quantity (De Koster
et al., 2007). Order batching is deﬁned as combining multiple customer orders into a single assignment for an order picker with the goal of minimiz-5

Introduction

ing travel or batch execution time. When optimizing for travel time, travel
speed is generally assumed to be constant to justify travel distance based
optimizations. Other components of total batch execution time, such as
pick and setup time, are also usually assumed to be constants, and thus
easily left out of optimization models. Wave picking is used if orders share
a common destination or completion deadline. Customer orders belonging
to a wave are processed, and often batched, together.
Customer orders can have precedence constraints, i.e., some items of theorder must be picked to the relevant container before others. The problemof batching orders and routing pickers while respecting the precedence
constraints of products is common, particularly in retail organizations,but such restrictions may also play a role in other warehouses (Dekker
et al., 2004, Chan and Kumar, 2008). Precedence constraints may varyin nature. They may be due to weight restrictions (heavy products at thebottom of the roll container), fragility (light at the top), shape and size
(big boxes at the bottom), stackability, but also preferred unloading sequence due to family grouping on the customer’s shelves. Although orderbatching and picker routing have received attention in the literature, thecombined problem of order batching and picker routing while respecting
the precedence constraints of the products (and in this thesis, including
potential multiple drop-off points in a route) has not. In the sample warehouse of this thesis, at most three customer orders (or orders for short)
can be batched in a single pick tour, each possibly for a different customer(i.e., a store) and a different delivery location within the warehouse. Eachcustomer order is to be preferably picked in a ﬁxed sequence, due to familygrouping in the company’s stores.
Order pickers are people with different skill sets and varying physicalcapabilities. Some are experienced and others are new at their jobs. Interviews with warehouse managers and Warehouse Management System(WMS) suppliers reveal that in addition to the skills and physical capabilities of a picker, work motivation plays a key part in the picker’s productivity. Besides batching and routing, the second key question of this thesis
is how to make the most of the available workforce. Using regression analysis on detailed log data provides a way of modeling the pickers based on
past performance as the data contain timestamps of all the pick events,
travel distances, and the mass and volume of picked items. Using models
forecasting the pickers’ batch execution times based the on log data, thestandard way of just minimizing travel time can be extended to optimiz-6

Introduction

ing the total order picking time. Furthermore, the models allow for better
assignment of work among the pickers based on the strong areas of theirpast performances, and thus may motivate the pickers more.

1.2

Objectives

The general objectives of this thesis are to improve the efﬁciency of pickerto-parts order picking where the batching of customer orders plays animportant part. In more detail, the objectives are:
• To include precedence-constrained routing in joint picker routing andorder batching problems.
• Improve the efﬁciency of state-of-the-art order batching algorithms in
the presence of precedence-constrained routing.
• To study order picker skill models based on operational data and toshow that pickers have signiﬁcant differences in skills, which affect expected job completion times based on the job parameters.
• To develop a methodology for order and batch assignment based on
skills of individual pickers and to show that signiﬁcant savings areachieved when compared to the current state-of-the-art.

1.3

Contributions of the Thesis

In this thesis, the following extensions to the related work presented in
Section 1.4 are made.
This thesis introduces precedence constraints on picker routing.
Precedence-constrained routing is combined with batching customer orders to ﬁnd the combination of batches that maximizes the total travel
distance savings of pick waves. As new batches are explored, their travel
distance savings need to be calculated. This is usually done by routing
the batches and comparing batched travel distance to the travel distanceto pick individual orders in the batch. Items must be picked in a strict
sequence to the customer orders. However, the complexity of precedence-

7

Introduction

constrained routing poses new challenges to order batching. Thoroughexploration of combinatorial batch spaces with state-of-the-art batchingalgorithms becomes impractical, as these algorithms need to route eachtime a new batch is found. A new order batching algorithm is introduced
to counter the time complexity of optimal precedence-constrained routing:a fast Large Neighborhood Search (LNS) algorithm, which uses an innovative way of estimating batch savings for batches of three or more customer orders to speed up the processing time. It is shown that this novel
way of estimating batch savings is accurate for those batches that savethe most in travel distance. Incidentally, the high-savings batches are
the ones that generate the best solutions. For the precedence-constrained
routing, an optimal A* algorithm is used.
In the presence of precedence constrained routing for medium to large
wave sizes, the batching algorithm, Precedence Constrained Estimated
Savings based batching (PCES), provides solutions of better quality than
any tested state-of-the-art order batching algorithm. Only one of thetested algorithms, C&W(i) (Clarke and Wright, 1964, De Koster et al.,1999), is consistently faster than PCES, but it is much worse in solution
quality in all cases. For small wave sizes, the mean error from optimal
is 1.2%; and for medium and large wave sizes, PCES outperforms othermethods by 1% to 3%. When compared to original batching and routing
by the company, more than 15% in total travel distance is saved. It is also
shown that the savings estimates and the batching part of PCES, SAREMIX, can be applied to batching problems when the routing does notcontain precedence constraints, with a similar solution quality as with theprecedences.
Second, using actual picking data, multilevel modeling is used to build
forecasting models for the pickers’ batch execution times. It is shownthat a signiﬁcant amount of the total variance in batch execution time,
13%, results from differences in picker skill. To the author’s knowledge,
this is the ﬁrst such application of using multilevel modeling to exploit
operational data, at least in the warehousing context.
Third, the varying skills of the order pickers are exploited by extending the standard set-partitioning model (Gademann and Van de Velde,2005) for the batching of orders with a generalized assignment problem.
An Adaptive Large Neighborhood Search algorithm (ALNS) (Ropke andPisinger, 2006a) is used to search a very large combinatorial space forsolutions to assign orders and batches to the right pickers. Via the combi-

8

Introduction

nation of the forecasting models and the ALNS algorithm, it is shownthat almost 10% of total batch execution time can be saved by taking
the skills of the pickers into account when compared to a state-of-theart batching solution, but without considering picker skill. Compared toassigning work based on picker productivity (picking speed), taking skill
into account saves 6% of the total time. To the author’s knowledge, this is
the ﬁrst work where worker-speciﬁc forecasting models are used in conjunction with a work-assigning algorithm. Also, the proposed model andproblem setup are novel. Finally, ALNS has not been previously applied
to order picking problems.

1.4

Related Work

This section presents relevant related work. First, routing methods inpicker-to-parts warehouses are presented, followed by routing with precedence constraints. Second, an overview of order batching methods is presented. Third, research where differences among workers skills and factors affecting workers’ performances are discussed. Fourth, an overview
of the class of Very Large Scale Neighborhood Search algorithms is presented. Multilevel modeling is brieﬂy touched upon. Finally, a study
where regression models of an order picking process have been built is
discussed.

1.4.1

Routing in Picker-to-Parts Warehouses

A classic algorithm for order picking tour construction is the exact andpolynomial time algorithm ﬁrst presented in Ratliff and Rosenthal (1983)
and extended in De Koster and Van der Poort (1998) to include multiple
drop-offs. The algorithm presented in Ratliff and Rosenthal (1983) is further extended in Roodbergen and De Koster (2001b) to include a middle
aisle. This extended version assumes a parallel-aisle warehouse with amaximum of three cross-aisles (see Figure 1.3) and does not account forprecedence constraints or unidirectional travel.
Next to optimal routes, heuristics are often used. In S-shape routing(Randolph, 1993), also known as traversal routing, an order picker travels
the whole subaisle if he or she needs to pick an item in the subaisle visited.
An exception is made for the last aisle if the number of aisles is odd. Withlargest-gap routing (Randolph, 1993), each subaisle, except for the ﬁrst

9

Introduction

Figure 1.3. Top view of a warehouse three cross-aisles (marked with arrows).

Figure 1.4. S-shape (left) and largest-gap (right) routing (Roodbergen and De Koster,
2001a).

and the last one visited, is exited on the entry side. Figure 1.4 shows
examples of the S-shape and largest-gap routing methods. The ﬁrst andthe last visited aisles are traveled completely. A heuristic combining Sshape and largest-gap routing is presented in Roodbergen and De Koster(2001a). In aisle-by-aisle routing (Vaughan, 1999), each aisle is visited
only once, i.e., ﬁrst all items in the ﬁrst pick aisle are picked, then all
in the second, and so on. See Figure 1.5, with the optimal solution as areference.
Theys et al. (2010) consider the applicability of the Lin-KernighanHelsgaun (LKH) TSP algorithm (Helsgaun, 2000) to routing order pickers.
They compare LKH with several routing heuristics and obtain savings of
up to 47% in travel distance. Helsgaun (2000) ﬁnds optimal solutions
with the LKH algorithm for all previously solved TSP instances available
at that time, including a 13,509-city problem, which was the largest problem instance solved to optimality at the time. LKH is used in Chapter 4 totest the SA-REMIX batching algorithm with non-precedence-constrained
routing and as a reference algorithm in comparing the effectiveness of a

10

Introduction

Figure 1.5. Aisle-by-aisle (left) and optimal (right) routing (Roodbergen and De Koster,
2001a).

routing heuristic used in Chapter 6.

1.4.2

Routing with Precedence Constraints

A dynamic programming algorithm for solving the single-vehicle manyto-many immediate request dial-a-ride problem (multiple customer destinations, each with a possibly unique drop-off location), which is similar
to the A* routing algorithm presented in Section 3.2.2, is introduced in
Psaraftis (1980b, 1983). In Psaraftis (1980b) each customer has threepossible states: not picked up, picked up, and delivered. In Section 3.2.2,
a single customer order is analogous to the customer, but with many more
possible states than three. Furthermore, Psaraftis (1980b, 1983) does notuse an A* type heuristic to speed up the search. The time complexity ofPsaraftis’ algorithm is O(n2 3n ) (where n is the total number of customers),
and it solves problems for up to ten customers.
Kubo and Kasugai (1991) introduce the Precedence-Constrained Traveling Salesman Problem (PCTSP) and a branch-and-bound method for ﬁnding exact solutions for cases of up to 49 locations with acceptable computation times. The precedence-constrained path construction problem can bemodeled as a case of the Sequential Ordering Problem (SOP) (Escudero,1988). The SOP can be formulated as an Asymmetric Traveling Salesman
Problem (ATSP) with precedence constraints. In the SOP, paths usually
have a start and a ﬁnish position that differ from each other, while ATSP
paths ﬁnish where they started. In the general ATSP case, each of the11Introduction

non-visited cities can be the next target with each iteration, but in the
SOP, the set of the next possible cities is limited as deﬁned by a directed
graph formed from the problem. The routing problem presented in Section
3.2.2 is a special case of the one presented in Kubo and Kasugai (1991)
since items in each customer order must be picked in a strict sequence.

1.4.3 Order Batching
The savings algorithm by Clarke and Wright (1964), C&W(i), as well as
an extension of it, C&W(ii), are used for batching orders in De Koster
et al. (1999). They are compared to seed algorithms using two routing
strategies: S-shape and largest-gap. Seed algorithms consist of two distinct steps: seed order selection and order addition. A single order isselected as the seed order based on criteria, e.g., the largest number of
items, the longest travel time, or the farthest item. Additions can be doneusing different rules such as adding the order that minimizes the sum of
the distances of every item of the seed and the closest item in the order,or minimizing the additional number of aisles to be traveled. The authors ﬁnd that seed algorithms work best with S-shape routing and large
pick device capacity, while savings algorithms work best with largest-gap
and small pick capacity. C&W(ii) consistently outperforms C&W(i), but iscomputationally more expensive.
Albareda-Sambola et al. (2009) use a Variable Neighborhood Search(VNS) algorithm to batch orders. It uses six different local exchange
schemes incorporated into three different search neighborhoods of varying
sizes to ﬁnd good batches. Within each neighborhood, all moves belonging to that neighborhood are tried, and the one that results in the largest
savings is chosen. The larger neighborhoods are searched as needed. If acurrent one fails to produce a better solution, the next (larger) neighborhood is explored for a better solution until no improvement can be made.
Albareda-Sambola et al. (2009) compare VNS to the C&W(i), C&W(ii), andseed algorithms, and it consistently outperforms them. The authors ﬁndthat the best performing algorithm from the literature is C&W(ii), which
is on average 2% worse than VNS. Solution quality comes with added computational complexity. When compared to C&W(ii), a much larger part ofthe combinatorial batch space is explored. For the most complex instancerun (250 orders), VNS took almost six times as much time to reach its
solution. For routing batches, the authors use the computationally inexpensive S-shape, largest-gap, and combined heuristics.

12

Introduction

Henn and Wäscher (2012) introduce an Attribute Based Hill Climber(ABHC) order batching algorithm (Whittley and Smith, 2004).

This

ABHC uses a set of attributes to guide the search out of local minima
and a set of moves to explore the neighborhood of a current solution. Twomoves are used to search for new solutions: shifting an order from one
batch to another batch and swapping two orders between two batches.
As the attribute that gave the best results, the authors used a pair ofcustomer orders. Each attribute is given a value, which initially equals
inﬁnity. As a solution (a set of batches) is found, it will have a total travel
distance. That solution is accepted if at least one of the attributes has acurrent value that is higher than the current solution’s value (the total
travel distance), after which that attribute’s value is set to the current
value. S-shape and largest-gap routing are used as routing algorithms.
Comparisons are made between three local search algorithms, a Tabu
Search heuristic, and C&W(ii). ABHC compares favourably to these algorithms, with over 2% better results than those achieved with C&W(ii).
A method for batching orders is introduced in Gademann and Van de
Velde (2005). The order-batching problem is modeled as a set-partitioning
problem. A column-generation algorithm is used to solve the linear programming relaxation. They rely on the polynomial time algorithm presented in Ratliff and Rosenthal (1983) to calculate the route length. Theyﬁnd that the maximum batch size has the largest impact on the solution
time.
Large-scale order batching in parallel aisle warehouses is consideredby Hong et al. (2012b). The authors introduce a route-selecting orderbatching (RSB) formulation, a route bin packing formulation for calculating lower bounds, and a heuristic route-packing based order batchingprocedure (RBP). Instead of constructing routes for batches, they use theS-shape routing method to construct a set of all possible routes given awarehouse layout to which batches can be assigned. Given a batch, a bestﬁt route can be selected with RSB. They assume that each order picker canpick up to ten orders to ten bins, and an average order size is two lines. A
comparison of RBP is made to C&W(ii) and seed algorithms, with the RBP
dominating the other algorithms and C&W(ii) performing better than theseed heuristics.
In comparison to the batching methods in the literature, which use multiple deterministic neighborhood search heuristics, the method presented
in Chapter 3, PCES, uses a single ﬂexible stochastic heuristic based on

13

Introduction

random moves, and is not constrained to a set of deterministic moves.
Other similar methods (Gademann and Van de Velde, 2005, AlbaredaSambola et al., 2009, Henn and Wäscher, 2012) do explicit batch evaluation (and thus routing), which is computationally expensive, while PCES
estimates the best batches using information inherent in batches of one
and two customer orders. Furthermore, the problem in Chapter 3 also
includes precedence constraints in the routing. With computationallymore expensive routing algorithms, such as those used in Chapters 3 and4, solving the batching with any such deterministic batching algorithm
becomes intractable due to exponential growth in calls to routing algorithms. In combination, the ﬂexible heuristic along with the estimate
enables escaping local minima and exploring larger parts of the set of
possible batches more thoroughly for good options. In addition to goodsolution quality, PCES is also computationally fast. With precedenceconstrained routing, it is slower than C&W(i) but faster than C&W(ii)and provides good solutions to instances that C&W(ii) cannot solve in areasonable time.
A simulated annealing (SA) algorithm for order batching is introducedby Hong et al. (2012a). The authors deal with narrow-aisle systems withthe possibility of picker blocking. In addition to the batching problem, theauthors also deal with the sequencing problem, which is also NP-hard.
They use four neighborhood search heuristics in their SA-algorithm: theﬁrst one exchanges two orders from different batches; the second swaps
the sequence of two batches; the third and fourth impose an acceptancecondition on the order exchange to limit the size of neighborhoods to be
searched. Their method aims to minimize total retrieval time, which is thesum of pick and travel time and time due to delays. The algorithm uses
a geometric cooling schedule (Cohn and Fielding, 1999), while the criteria
for accepting a worse batch are the same as PCES. Hong et al. (2012a) use
computationally light traversal routing, and thus have no need to makeuse of savings estimates. Similarly, as in Chapter 6, the simulated annealing algorithm presented in Hong et al. (2012a) solves a combinationof two NP-hard problems. They consider assignments to a set of homogeneous pickers, while the main contribution of Chapter 6 is to extend theorder batching problem to include the (generalized) assignment of orders
and batches to order pickers of varying skills.
Hsieh and Huang (2011) introduce two batch-construction heuristics
based on data-clustering methods: K-means Batching based on the14Introduction

K-means algorithm by MacQueen (1967); and Self-organisation Map(SOMB) based on the Self-Organizing Map (Kohonen, 1990). The Kmeans algorithm can be classiﬁed as a seed algorithm. It uses a selection rule of selecting the order with the minimum travel distance and an
addition rule of adding an order which causes the minimum incrementin travel distance when combined with the existing batch. The process
is greedy and sequential and very similar to seed rules presented previously in De Koster et al. (1999). The SOMB algorithm groups orders by anorder-relativity measure using the SOM. The measure is a weighted sumof the number stock-keeping units (SKUs) in an order and the number ofthe same aisles covered between two orders.
Of these algorithms, C&W(i) and C&W(ii), the VNS of AlbaredaSambola et al. (2009), and the ALNS by Henn and Wäscher (2012) areused in comparison algorithms to the batching algorithm in Chapter 3.
C&W(i) is used in constructing initial batching solutions in chapters 3
and 6, and the VNS in Chapter 6. These algorithms are presented in
more detail in Chapter 2.

1.4.4

Heterogeneous Workers and Jobs

Assigning jobs to workers based on their capabilities has been studied
relatively well. A distinction is made between capability (the degree towhich a person is able to carry out a certain task, regardless of the time
it takes) and skills (the speed at which such a task can be carried out).
In this thesis, all pickers have the capability to pick the batches, albeit
at different speeds. Differences in skills may be caused by, the ability
to lift heavy objects, motivation, stacking ability, and knowledge of thewarehouse — many of which are not directly quantiﬁable.
Taking advantage of workers’ capabilities has been studied and modeled
as a generalized assignment problem (GAP). For extensive surveys on thesubject, see Cattrysse and Van Wassenhove (1992) and Pentico (2007).
Campbell and Diaby (2002) solve a GAP for cross-trained workers withan assignment heuristic. Each worker has a predetermined value fromzero to one, characterizing his or her ability to work in a department (todo a job). Chapter 6 extends the order batching problem with the GAP toinclude knowledge of the pickers’ skills in batching and assigning orders
to pickers.
Another stream of research investigates how worker speed depends on
endogenous and exogenous factors (Bendoly et al., 2006). Powell and15Introduction

Schultz (2004) show that workers increase speed in the presence of a visible backlog. Doerr and Arreola-Risa (2000) investigate how the variabilityof task completion times is affected by different tasks and workers. Theyﬁnd that the most signiﬁcant factors affecting the variability are: (1) the
worker who did the work and (2) the interaction effect between the worker
and the task, while the task type in question had no signiﬁcant inﬂuenceon its own. Juran and Schruben (2004) use personality and demographic
data to predict task execution times of a collaborative two-worker task.
They ﬁnd that including information on individual workers has signiﬁcant impact on simulation accuracy compared to assuming no differencesamong workers. Bartholdi and Eisenstein (1996) show that if workers
are sequenced from slowest to fastest in typical production line, a stable partition of work will spontaneously emerge, with the production rate
converging to the maximum possible value. These studies reinforce thehypothesis of picker skills playing a signiﬁcant role in total batch execution time, which is further studied in chapters 5 and 6.

1.4.5 Very Large Scale Neighborhood Search
Very Large Scale Neighborhood Search (VLSN) is a class of neighborhood
search algorithms, where the searchable neighborhood size grows exponentially with the instance size or the neighborhood is too large to be
searched explicitly (Pisinger and Ropke, 2010). Ahuja et al. (1998) classify
VLSN methods into three categories: variable-depth methods, networkﬂow based improvement methods, and methods based on subclasses solvable in polynomial time. Pisinger and Ropke (2010) extend these categories by including Large Neighborhood Search (LNS) algorithms.
The joint batching and generalized assignment problem dealt in Chapter 6 belongs to this class, as does the order batching problem when batch
size grows, especially when batch evaluation is computationally expensive, which imposes a practical limit on explicitly searching the whole
neighborhood (see Chapter 3). To solve the order batching part in Chapter 3, a LNS-type algorithm is applied. In Chapter 6, an Adaptive Large
Neighborhood Search (ALNS) algorithm (an extension of LNS) is used tosolve a joint batching and generalized assignment problem. The LNS andALNS algorithms are detailed more thoroughly in sections 2.2.5 and 2.2.6,
respectively.

16

Introduction

Large Neighborhood Search
According to Pisinger and Ropke (2010), neighborhoods in LNS (Shaw,
1997) methods are deﬁned by destroy and repair methods. They are also
called in removal and insertion methods in other sources, e.g., in Ropke
and Pisinger (2006a) and Gharehgozli et al. (2013). The destroy methodsof neighborhood search heuristics break down a part of a solution in a predeﬁned manner, while the repair methods try to ﬁnd good alternatives toa solution by constructing another feasible solution in the neighborhoodof the solution based on the repair method associated with the neighborhood search heuristic. These methods typically contain some element ofrandomness in choosing what parts of a solution to modify. The degree ofdestruction (Pisinger and Ropke, 2010) of the destroy method should be
carefully chosen. Destroying too small parts of a solution can result inmyopic exploration, while too large changes can result in suboptimality
and repeated re-optimization.

Adaptive Large Neighborhood Search
Adaptive Large Neighborhood Search, ﬁrst proposed in Ropke andPisinger (2006a) to solve the Pickup and Delivery Problem with Time
Windows, extends the LNS by adding a learning component to the search.
Typically, ALNS uses several heuristics to form the search neighborhood.
Details of each heuristic’s runtime performance are kept in memory. Eachheuristic has a score, which gets updated positively if good solutions arefound and/or the heuristic ﬁnds new interesting solutions. One heuristic
is selected with each iteration, with a weighted random selection criterion called the roulette wheel method (Pisinger and Ropke, 2010). As with
the LNS, a heuristic is generally composed of a combination of destroy
and repair methods. However, in Section 6.3.4, dedicated heuristics areused, which are not naturally divisible into both destroy- and repair-type
methods. The set of all heuristics H is formed by (not necessarily all)combinations of the destroy and repair methods and the possible dedicated heuristics. ALNS has been shown to work very effectively on a wide
variety of computationally hard problems, such as container scheduling(Gharehgozli et al., 2013), vehicle routing (Ropke and Pisinger, 2006a)
and scheduling technicians and tasks with varying skills (Cordeau et al.,2010).

17

Introduction

1.4.6 Multilevel ModelingMultilevel modeling is a regression analysis technique, that can be used
to distinguish group level effects in some population. A standard linear
multivariate regression can be extended to a linear multilevel model byconsidering varying intercepts and slopes for each group. Fixed effects
of a multilevel model contain the overall (mean) model encompassing information from all groups, quite similarly as a single multivariate model
would. In linear multilevel regression, the random effects or variations tothe ﬁxed effect intercept and slope are calculated for each group. In thismanner, forecasting models with group-speciﬁc parameters can be found.

1.4.7 Forecasting Order Picking Time
Larco Martinelli (2010) studies the problem of how to assign products
to different storage locations, with the objective to minimize both order picking cycle time and worker discomfort. By reducing discomfort
felt by the workers, the well-being of the workers increases, and may
also yield the long-term beneﬁts of increased productivity (Kuijt-Evers
et al., 2007). Larco Martinelli (2010) uses extensive warehouse log data
to build regression models for both objectives. Subsequently, the models
are used in a dual-objective assignment model to optimize product placement. Larco Martinelli (2010) ﬁnds that there is an efﬁcient frontier between the two objectives. If either cycle time is solely prioritized, worker
discomfort is high and vice versa.

1.5

Outline

Chapter 2 presents algorithm descriptions and mathematical models related to the subject. Chapters 3 to 6 contain the scientiﬁc contributions of
this thesis. Chapter 7 concludes the thesis.
The main contribution in Chapter 3 is the modeling and solving of a
joint precedence-constrained picker routing and order batching problem.
A journal article based on this chapter has been published: (Matusiak
et al., 2014). These results are extended in Chapter 4 by relaxing theprecedence constraints in the routing. In Chapter 5, the large order picking data-set obtained from the Finnish retailer is processed and used toform forecasting models for the execution times of batches for pickers. The

18

Introduction

forecasting models constructed in Chapter 5 are utilized in Chapter 6 tosolve a joint-routing, order batching and generalized assignment problem.

1.6

Author’s Contribution

Chapter 2 presents mathematical models and algorithms relevant to thisthesis. All the models in that chapter are from the literature. The algorithm in Section 2.2.1 has been developed by the author. The otheralgorithms presented in Chapter 2 have been implemented by the author,
as detailed in the relevant references. Chapters 3 to 7 hold the main contributions of the thesis. These chapters are the work of the author, apartfrom the exceptions noted below. The problem description and mathematical model presented in Section 3.1 and the CPLEX code used in a partof the experiments in Chapter 3.3 were formulated and coded by Professor Leo Kroon. The Lin-Kernighan-Helsgaun algorithm for solving theTraveling Salesman Problem, used in Chapter 4 and Chapter 6, has beendeveloped by Keld Helsgaun and is free for research use.

1.7

Assumptions

There are number of assumptions made in this thesis. In Chapter 3, thepickers are always assumed to travel with a constant speed, the time topick an item to the container is assumed constant, as is the setup time(getting on and off the truck). Congestion is not considered in the optimization, as are no other dynamic factors. Also, it is assumed that theorder pickers are alike. The mathematical model does not consider time
windows for the deliveries’ of the orders, as the log data does not contain
information on the truck departures.
In Chapter 5, the standard assumptions of non-multicollinearity, normality, and non-heteroscedasticity are made in the regression. However,the data are heteroscedastic, which is partly compensated by the logarithmic transformations applied to the data.
Some of the assumptions in Chapter 6 are the same as in 3: congestion, deadlines for orders and system dynamics are not considered. While
batches are assigned, they are not sequenced in any manner. Stochasticity in the picking process is not taken into account, i.e., the problem issolved in a deterministic manner not considering the uncertainties that19

Introduction

could be extracted from the regression model.

202. Relevant Mathematical Models andAlgorithms

This chapter is organized as follows. In the ﬁrst part of this chapter,
combinatorial optimization problems are formally deﬁned. The presented
nomenclature is used in problem deﬁnitions later. Three separate combinatorial optimization problems are described to illustrate subproblems,
each of which plays a part in the joint models presented in chapters 3,
4 and 6. These are the traveling salesman problem, the order batchingproblem, and the generalized assignment problem.
The second part of this chapter starts with descriptions and pseudocode
of some of the algorithms implemented for comparison purposes and initial solutions. Finally, Large Neighborhood Search and Adaptive Large
Neighborhood Search metaheuristics are detailed, which are later applied
to order picking problems in chapters 3, 4 and 6.

2.1

2.1.1

Mathematical Models

Combinatorial Optimization Problems

Following Pisinger and Ropke (2010), let I be an instance of a combinatorial optimization problem, where S is the set of feasible solutions forthe problem. Let f : S → R be a function which maps a solution toits cost. A globally optimal solution s∗ is such that cost is minimized:f (s∗ ) ≤ f (s),∀s, s∗ ∈ S. A neighborhood of a solution s ∈ S is deﬁned

as N (s) ∈ S. N is function which maps s to a set of solutions S  ⊆ S. If
f (s ) ≤ f (s), ∀s ∈ N (s), then s is said to be a locally optimal solution forneighborhood N (s).
A neighborhood search heuristic h(s) ∈ H, or a move, is a function whichenables searching of the neighborhood Nh (s) of a solution s ∈ S. H is theset of all neighborhood search heuristics of an algorithm.

21

Relevant Mathematical Models and Algorithms

A steepest descent algorithm is an algorithm explores N (s) for solution
s. In this case, N (s) is searched for the best solution s ∈ N (s), i.e., s =arg maxs ∈N (s ) f (s ). If f (s ) < f (s), an update s = s is performed.

2.1.2 Routing
The problem of routing order pickers in conventional multi-parallel-aisle
warehouses is often modeled as a Steiner Traveling Salesman Problem(Steiner-TSP) (Cornuéjols et al., 1985, De Koster et al., 2007, Theys et al.,2010). The Steiner-TSP includes Steiner and non-Steiner nodes: Steiner
nodes can be visited multiple times, and are typically located at the end of
aisles. Each non-Steiner node represents a Stock Keeping Unit’s (SKU) location from which items need to be picked during the current tour. Theyset al. (2010) propose that for warehouses with unconventional layouts, theSteiner nodes can be discarded. Let G  be the set of items to picked during
a tour. Among other things, an item contains information about locationof the products to be picked. Let D denote the depot in the warehouse,and let G = G  ∪ D. Also, let Xi,j be a binary decision variable, whichequals one if item i is picked directly before item j, and zero otherwise.
The problem can now be represented as a standard TSP:
min

 

ci,j Xi,j ,(2.1)

i∈G j∈G\isubject toXi,j = 1 ∀j ∈ G \ i,

(2.2)

i∈G



Xi,j = 1 ∀i ∈ G \ j,

(2.3)

j∈G

U1 = 1,2 ≤ Ui ≤ |G|(2.4)∀i = 1,

(2.5)

Ui − Uj + 1 ≤ (G − 1)(1 − Xi,j ),(2.6)

Xi,j ∈ {0, 1}(2.7)∀i, j ∈ G,

where (2.2) and (2.3) guarantee that each item location is entered andexited once, respectively. Equations (2.4), (2.5), and (2.6) are the subtour
elimination constraints of Miller et al. (1960). If for some i, j ∈ G, ci,j =cj,i , the problem is said to be asymmetric; otherwise it is symmetric.

22

Relevant Mathematical Models and Algorithms

2.1.3

Order Batching

A customer order o ∈ O consists of a subset Go of the set of all items (order lines) G to be picked during some period of time, where O is the setof all customer orders. The items in each order have to be picked to thesame container, or in the case of the sample retailer, to a roll container.
A combination of customer orders that can be picked together during apick tour by one order picker is a batch, and the set of all possible batches
is denoted by R. For each feasible batch r ∈ R, a version of a Traveling
Salesman Problem (TSP, see Section 2.1.2) can be deﬁned. These TSPs
are usually solved with heuristics (Ratliff and Rosenthal, 1983, Randolph,
1993, Theys et al., 2010), or in some special cases with optimal algorithms(Ratliff and Rosenthal, 1983, De Koster and Van der Poort, 1998, Roodbergen and De Koster, 2001b). In the commonly used order batchingmodel of Gademann and Van de Velde (2005), presented below, solving
the TSP is included implicitly only as the travel cost. This is not exactly
correct, but it does make the model much simpler when the focus is not
so much on the routing, but on the order batching algorithm. In Chapter
3, a joint precedence-constrained TSP and order batching problem is presented, while in the model of Chapter 6 the TSP is similarly absent fromthe model as below.
Following Gademann and Van de Velde (2005), the order batching problem, which is central for the research presented in this thesis, can be modelled as a set partitioning problem. Every batch r ∈ R is characterized bya zero-one vector ar = (a1,r , ..., a|O|,r ), where ao,r = 1 if and only if ordero ∈ O is included in batch r. Batch r is feasible only if o∈O ao,r ≤ N ,where N is the maximum customer order capacity of the picking device(the maximum batch size). Let cr be the cost to pick all items of batch r.
Furthermore, let Xr be a zero-one variable, such that

Xr =

⎧
⎪
⎨1

if batch r is selected⎪⎩0 otherwise.

(2.8)

In order to minimize the total cost to pick all batches, the model can be23Relevant Mathematical Models and Algorithms

described as follows:
min



(2.9)

c r Xr ,r∈Rsubject toao,r Xr∀o ∈ O,(2.10)∀r ∈ R,(2.11)r∈R

Xr ∈ {0, 1}

where (2.10) guarantees that each order is assigned to one batch and one
batch only. Gademann and Van de Velde (2005) show that for batch sizes
larger than two (N > 2), the order batching problem is NP-hard.

2.1.4 Generalized Assignment Problem
The Generalized Assignment Problem (GAP) is an NP-hard optimization
problem (Fisher et al., 1986). Let tw,r be the cost for agent w ∈ W tocomplete job r ∈ R. Let Xw,r be a binary decision variable, which is one
if agent w is assigned to job r, and zero otherwise. Furthermore, let aw,rbe the capacity absorption when assigning job r to agent w, and Mw is the
maximum capacity of agent w. The objective is to minimize the total costof assigning all jobs to agents. The problem can be formulated as follows:
 

min

tw,r Xw,r ,(2.12)w∈W r∈R

subject toaw,r Xw,r ≤ Mw

∀w ∈ W,(2.13)r∈R



Xw,r = 1 ∀r ∈ R,

(2.14)w∈W

Xw,r ∈ {0, 1}∀w ∈ W, ∀r ∈ R,

(2.15)

where (2.13) is to enforce the capacity constraints for each agent and(2.14) guarantees that all jobs are assigned each to one agent. In thecontext of assigning batches to pickers, each job is a batch, and an agent
is a picker. Furthermore, in this thesis, aw,r = tw,r , i.e., the capacity absorption for an agent to complete a job is the same as the cost to completeit.

24Relevant Mathematical Models and Algorithms

2.2

2.2.1

Algorithms

Dynamic Programming Algorithm for Optimal BatchingThis algorithm searches for an optimal assignment of orders to batches.
It is used to ﬁnd optimal batching solutions for instances of limited size insections 3.3 and 4.2. Let O be the set of orders to be processed, and N the
maximum batch size, with |O| being the total number of orders. Each nodein the algorithm is a pair composed of a state a, indicating which orders
are included in a solution, and the state’s cost ca . The algorithm operates
using a priority queue-memory structure. As new nodes are added, thepriority queue is sorted in ascending order by cost. Thus the lowest cost
node is always on top. As the state, a binary vector a = {a1 , ..., a|O| }, isused. If ai = 1, order i is included in a node’s solution, otherwise not. As
the algorithm runs, and a new state an+1 is reached, a map data structure m is used to keep track of the state of the entering state an so thatm(an+1 ) = an . To ﬁnd the (binary) batch r that was added to an to reachan+1 , the following formula is used: r = an ⊕ an+1 , where ⊕ is the binary
exclusive-or operator. The steps of the algorithm are detailed below.
Step 0. Add a node with zero cost and state a = {0, ..., 0} to the priority
queue.
Step 1. Read the top node p, and its state a.
Step 2. If p contains all orders, extract the batch combinations in thefound solution from m, return p and exit.
Step 3. Find the smallest index i in a that equals zero.
Step 4. Find the set of feasible batches Ri that include order no. i, but
are not yet included in a (up to batch size of N ).
Step 5. ∀r ∈ Ri generate new states ar .
Step 6. ∀r ∈ Ri ﬁnd the costs car of ar .
Step 7. Remove p from the priority queue.
Step 8. ∀r ∈ Ri add new nodes pr = {ar , car } to the priority queue, andset m(ar ) = a Go to Step 1.

25

Relevant Mathematical Models and Algorithms

2.2.2 Clarke and Wright Based Algorithms
The order batching variants of the Clarke and Wright (1964) routing algorithm, C&W(i) and C&W(ii) (De Koster et al., 1999), are presented below.
Both algorithms operate by updating a savings matrix, iteratively adding
feasible pairs of orders or order clusters to the current solution. The savings matrices indicate the savings if two orders or groups of orders arebatched together. Generally, the algorithms aim to add a pair of orders
or order clusters with the maximum savings value in the matrix to thecurrent solution, if feasible.

Basic variant: C&W(i)This basic variant of the C&W algorithm provides a fast solution for thebatching problem, which is very usable for initial solutions and subsequent local search. It requires the routing of all combinations of one andtwo customer orders from the set of all orders O, which is usually computationally manageable (see Section 3.2.1). The algorithm consists of thefollowing steps (De Koster et al., 1999):
Step 1. Calculate the savings for all possible customer order pairs{i, j} ∈ O, where O is the set of all customer orders.
Step 2. Sort the savings in descending order.
Step 3. Select the pair with the highest savings.
Step 4. Three cases can be distinguished:(i) Neither order has been previously included in the current
solution s. Add a batch consisting of orders i and j to solutions.
(ii) Either i or j is included in s, but not both. Maximum batch
capacity permitting, add the as-of-yet unincluded order to thebatch of the included order. Otherwise proceed to Step 5.
(iii) Both orders have already been included in batches: proceed to
Step 5.
Step 5. Select the order combination with the next highest savings, andproceed to Step 4.
If all customer orders o ∈ O are not yet present in solution s, add a single
order batch for each remaining order. Finally, to get the total routing cost,
the batches in the ﬁnal solution are routed.

26

Relevant Mathematical Models and Algorithms

Recalculation of the Distance Matrix: C&W(ii)The C&W(ii) algorithm considers order clusters as orders and replaces
steps 4 and 5 with the following:
Step 4. Combine order clusters to form a new cluster, if allowed by the
maximum batch size. If not, choose the next combination andrepeat Step 4.
Step 5. If all orders have not been included in the current solution,
proceed with Step 1, otherwise ﬁnish.
C&W(ii) is computationally much heavier than C&W(i), as it requires
evaluation (and thus routing) of batches consisting of more than two customer orders inside the procedure to recalculate the distance matrix, instead of just routing batches as the algorithm ﬁnishes.

2.2.3

Variable Neighborhood Search

The VNS algorithm presented in Albareda-Sambola et al. (2009) uses
three neighborhoods of different sizes to search a new solution in theneighborhood of solution s ∈ S. The ﬁrst neighborhood, N1 , is composedof a single heuristic move, which moves one order from a batch to another. The second neighborhood, N2 , is composed of N1 and two additional
moves: one that transfers two orders from one batch to another, and onethat transfers two orders to different batches. Neighborhood N3 includes
N2 in addition to three other moves: the ﬁrst swaps orders between two
batches, the second transfers two orders from two separate batches to one
other batch, and the third moves an order from batch a to batch b and anorder from b to batch c. For a given instance I of a combinatorial problem, the VNS always reaches the same solution; there is no randomness
involved in the search. See Figure 2.1 for an illustration: ra , rb , rc ∈ R represent batches involved in each iteration of a search heuristic, oi , oj ∈ O
the related orders, while the heuristics are denoted by hi ∈ H, wherei ∈ {1, ..., 6}.
When the algorithm runs, there must exist at least one empty batch sothat it is always possible to transfer an order from a batch. Batch capacity
constraints are always maintained, and the feasibility of the solution ispreserved throughout the search. The ﬂow of the algorithm is detailed
below in the following steps.
Step 0. Form an initial solution s0 by including one batch per order in27

Relevant Mathematical Models and Algorithms

Figure 2.1. VNS neighborhoods and search heuristics. Adapted from Albareda-Sambola
et al. (2009).

addition to an empty batch. Set s = s0 .
Step 1. Find the best neighboring solution s in N1 (s). If f (s ) < f (s), sets = s , repeat Step 1; otherwise go to Step 2.
Step 2. Compute the savings values for merging all batch combinations.
If a solution with positive savings values exists among all
combinations, ﬁnd the solution s where the merging of two batches
will result in the biggest savings, set s = s , and repeat Step 2.
Otherwise proceed to Step 3. This step is similar to applying
C&W(ii) to the batching problem.
Step 3. Find the best neighboring solution s in N1 (s). If f (s ) < f (s), sets = s , repeat Step 3; otherwise go to Step 4.
Step 4. Find the best neighboring solution s in N2 (s). If f (s ) < f (s), sets = s , repeat Step 4; otherwise go to Step 5.
Step 5. Find the best neighboring solution s in N3 (s). If f (s ) < f (s), sets = s , go to Step 4; otherwise exit.
To ﬁnd the total routing cost, the sum of all costs of the batches batches
present in s is calculated.

28

Relevant Mathematical Models and Algorithms

2.2.4

Attribute Based Hill Climber

The Attribute Based Hill Climber (Whittley and Smith, 2004) can be classiﬁed as a Tabu Search-type algorithm. Instead of using a classic Tabu-listof visited solutions which should be avoided, a memory of attributes andtheir corresponding values is maintained.
Let s ∈ S be the current solution, and s ∈ S be the new candidate solution. For each problem, a set of neighborhood search heuristics (moves) is
deﬁned. Two sets of attributes can be distinguished by applying a moveto reach solution s from solution s: the entering (s − s) and the leaving(s − s ). The entering attributes form the set of attributes that do not exist
in solution s but are present s , and vice versa for the leaving attributes.
Regardless of the type of combinatorial problem, the ABHC can be implemented with three algorithms, presented below (Whittley and Smith,2004). Algorithm 1 returns a boolean value for the acceptance of candidate solution s , Algorithm 2 contains the attribute memory update function, and Algorithm 3 is a "steepest ascent/mildest descent" local searchalgorithm.
In Algorithm 1, line 1, a candidate solution s is accepted if it contains an
attribute for which a bigger value has previously been stored in memory.
Otherwise if the new solution has a bigger cost (line 3), the new solution isrejected. Finally, if the cost of s is higher than that of the biggest (worst)
cost in s, it is rejected. Otherwise, if moving from solution s to s deletes
one of the worst attributes in solution s (lines 6-11), s is accepted.
Input: Solutions s, s , attribute memory amem
Output: Acceptance of solution s
1 if ∃a ∈ (s − s) s.t. f (s ) < amem[a] then2
return true;
3 else if f (s ) ≥ f (s) then4
return false;
5 else6for i = 1 to |s − s | + 1 do
7
a = worst(i, s);8
if f (s ) ≥ amem[a] then9
return false;
10
else if a ∈/ (s − s ) then11
return true;
12
end
13
end
14 endAlgorithm 1: ABHC acceptance
Algorithm 2 updates the attribute memory as needed only: if the cost of29

Relevant Mathematical Models and Algorithmss is higher than that of s, only the new (entering) attribute values need
to be updated (lines 1-2). Otherwise, all of the attributes are updated.
Input: Solutions s, s , attribute memory amem

1 if f (s ) ≥ f (s) thenamem[a] = min(amem[a], f (s )) ∀a ∈ (s − s);
3 else4amem[a] = min(amem[a], f (s )) ∀a ∈ s ;5 end
Algorithm 2: ABHC memory update2

A single local search iteration of the ABHC is described in Algorithm 3.
On line 1, the worst |s − s | + 1 attributes in s are identiﬁed and stored.
The cost of the current best solution (best) is initialized inﬁnity on line
2. The neighborhood of s is iterated through updating the best solution
when a better solution is found and accepted (lines 3-7). Finally, updates
to the attribute memory are done on line 9 if a new candidate solution isfound in N (s).
Input: Current solution s, attribute memory amem
Output: Updated solution s
1 identifyWorstAttributes(s)2 best = null {initialize f (best) = ∞}3 for all s ∈ N (s) do
4
if f (s ) < f (best) and accepts(s , s) then5best = s
end
6
7 end
8 if best = null then9
update(s, best)10s = best
11 endAlgorithm 3: ABHC main routineHenn and Wäscher (2012) use ABHC to solve the order batching problem. Two types of attributes are explored: the presence of a pair of orders{o1 , o2 } in any batch, and the assignment of an order o to a batch r (i.e.,
the pair {o, r}). They ﬁnd that using the pair of orders as the attributeleads to better performance in their algorithm. Two deterministic heuristics form the neighborhood of a solution s: a swap of two orders between
batches and a move of an order from one batch to another (h4 and h1 in
Figure 2.1, respectively).

30Relevant Mathematical Models and Algorithms

2.2.5

Large Neighborhood Search

According to Pisinger and Ropke (2010), neighborhoods in LNS (Shaw,
1997) methods are deﬁned by destroy and repair methods. They are also
called in removal and insertion methods in other sources, e.g., in Ropke
and Pisinger (2006a) and Gharehgozli et al. (2013). The destroy method ofa neighborhood search heuristic h ∈ H breaks down a part of a solution sin a predeﬁned manner, while the repair methods try to ﬁnd good alternatives to s by constructing another feasible solution s in the neighborhoodof s based on the repair method associated with h. These methods typically contain some element of randomness in choosing what parts of s tomodify.
Let i(s ) denote a repair and e(s) destroy function of a LNS algorithm
for a partial solution s ∈/ S and solution s ∈ S, respectively. Algorithm
4 shows the pseudocode implementation of a LNS algorithm. As input,
a feasible starting solution is needed. The algorithm iteratively applies
destroy and repair methods in succession, and it stores new best known
solutions until an exit criterion is met (such as a cap on the number of iterations). The acceptance function on line 4 can take the form of a greedy
heuristic or something else, e.g., a simulated annealing acceptance criterion (see Section 6.3.4).
Input: Feasible solution s
Output: Best found solution s
1 sbest = s
2 while Stop criterion is not met do3s = i(e(s))
4if accept(s, s ) then5s = s
end
6
7
if f (s ) < f (sbest ) then8sbest = s
9end
10 end11 return sbest
Algorithm 4: Large Neighborhood Search (Pisinger and Ropke, 2010)

2.2.6

Adaptive Large Neighborhood Search

The probability of choosing heuristic h ∈ H with the roulette wheel
method is
qh = θhh ∈H θh,

(2.16)

31

Relevant Mathematical Models and Algorithms

where θh ∈ R+ is the weight of heuristic h. As the algorithm runs, someheuristics will perform better than others. Good performance of a heuristic will increase its weight and thus the probability of it being selectedagain by (2.16). Each of the heuristics has a score, which is updated whenever the heuristic is selected based on the heuristic’s search success. Eachscore Ψh is updated by:
⎧
⎪
σ1
⎪
⎪
⎨
Ψ h = Ψ h + σ2
⎪
⎪
⎪
⎩
σ3

if the new solution is the global best.
if the new solution is better than the current one.
if the new solution is accepted.
(2.17)

An update to a heuristic’s weight θh is done every δ iterations byθh = λ/ζh θh + (1 − λ)Ψh ,(2.18)

where λ is the decay parameter, which controls the rate of change in theweights, and ζh ≥ 1 is a normalization factor which reﬂects the computational effort that h requires (Ropke and Pisinger, 2006b). At the beginningand after an update of the scores, i.e., every δ iterations, the values Ψh are
set to zero ∀h ∈ H. If a hill climbing scheme is not used in the algorithm
updates to the scores occur only when a new global best solution is found,i.e., Ψh can only get the value σ1 .
A pseudocode for the ALNS is shown in Algorithm 5 (Pisinger andRopke, 2010). The main differences between algorithms 5 and 4 are theinclusion and updates of the heuristics’ selection probabilities on lines 2
and 12, respectively. Additionally, instead of explicitly referring to eachpair of destroy and repair methods, their combinations are marked moregenerally with h.

2.3

Summary

Each of the presented combinatorial optimization problems is NP-hardby itself. Combining this with the realistic instance sizes dealt with in
this thesis, which in most cases result in very large search spaces and/orlong computational times, local non-optimal methods are in practice theonly option to solve the problems. The presented algorithms, apart fromthe C&W methods and the optimal batching algorithm, can be classiﬁed
as local search methods. All the other algorithms are deterministic in32

Relevant Mathematical Models and Algorithms

Input: Feasible solution s
Output: Best found solution s
1 sbest = s
2 ph = {1, ..., 1}
3 while Stop criterion is not met do
4select neighborhood search heuristic h ∈ H using ph5s = h(s)6
if accept(s, s ) then7s = s
end
8
9
if f (s ) < f (sbest ) then10
sbest = s
end
11
12
update ph13 end14 return sbest
Algorithm 5: Adaptive Large Neighborhood Search (Pisinger andRopke, 2010)
the sense that their exploration does not depend on any random routine.
For reasonably sized problems with no computationally expensive subprocedures that must be run when new solutions are found, they function
very well, as can be seen from the following chapters. However, whensearch spaces grow in size, such as when considering two of the presented
problems in unison, or if they have expensive subprocedures, deterministic search algorithms cannot explore the space of solutions in reasonable
time.

33

Relevant Mathematical Models and Algorithms

34

3. Order Batching when Items Have to
Be Picked in a Strict Sequence

Batching orders and routing pickers while respecting precedence constraints is complex.

Order batching is NP-hard for batches of size

three or higher (Gademann and Van de Velde, 2005). Routing pickers
while respecting precedence constraints has complexity of O(N 2 (P + 1)N )(Psaraftis, 1980a), where N is the number of customer orders in the batch
and P is the maximum number of items in an order. Finding exact solutions for large orders and large batches becomes rapidly intractable.
The company used as a reference in this thesis has 2,000 orders per day
and up to 50 products per order, which have to be batched in batches ofthree orders. Batching aims to reduce total travel distance by combining multiple similar orders (which have to visit overlapping parts of thewarehouse). As the problem has to be solved multiple times per day andbecause of large and exponentially growing search spaces and relatively
time consuming routing exact computation of the order batching is not
feasible and thus a heuristic method is used.
The problem is solved using a generic solution method. For the batching, a simulated annealing-based combinatorial search algorithm basedon maximizing total savings in travel distance is used. Selecting theorders to be batched is based on comparing all pairs of orders and estimating the savings due to combining a larger number of orders. This
reduces computation times substantially, while the error compared to optimal combinations is within 1.2%. The routing can be solved using anA∗ -type shortest path algorithm. The method is exact and general, notconstrained by the particular warehouse layout. The method is readily
applicable to any warehouse where batch picking is applied in waves andthe orders have strict internal precedence constraints.
The contribution of the research presented in this chapter is a structured approach for the joint order batching and picker routing problem

35Order Batching when Items Have to Be Picked in a Strict Sequence

Figure 3.1. Top view of the warehouse with seven drop-off locations indicated at the top(marked 1-7). The empty bin depot is denoted with D. The 57 aisles are
unidirectional, while the three cross-aisles are bidirectional. The aisle-toaisle distance is 5.5 m the slot-to-slot distance is 3.7 m, which is also thewidth of the central cross-aisle.

in warehouses with any layout or any strict sequence constraint using
an innovative way of estimating batch savings. The method is fast, andcompares well to optimal solutions and heuristics from literature. In thereference case, travel distance savings of nearly 16% could be achieved.

3.1

Problem Description

In this section the combined batching and routing problem as a variantof the precedence-constrained Traveling Salesman Problem is described.
The process is as follows. A given number of trucks pick the customer
orders. Each customer order represents a roll container (a bin in this context) to which items must picked to in a strict pre-speciﬁed sequence. The
items do not have to be stacked on top of each other in this order, if thereis space in the bin to do otherwise. Each truck starts completely empty
and then receives a set of orders to be picked. The number of received
orders per truck does not exceed N , the number of bins carried by a truck.
Before the items of an order can be picked, for each order an empty bin
has to be collected at the empty bin depot (see Figure 3.1). Subsequently,
each truck travels through the warehouse to pick the items of the orders
in the pre-speciﬁed sequence per order. Once all items of an order havebeen picked, the truck delivers the order at its drop-off location. Once allorders have been completed and the truck is empty again, it receives a
new batch of orders. A batch is the set of orders handled in a route by atruck between the moment the truck receives this set and the moment it
has dropped off all orders and is empty again.
Given the assignment of orders to batches, the sequence in which thebatches are carried out is not relevant for the total picking time.

36Order Batching when Items Have to Be Picked in a Strict Sequence

3.1.1

Notation

There is a set O of orders to be carried out. Each order consists of anumber of items that have to be picked at speciﬁc pick locations. All itemsof an order are picked in a roll container of one of the available trucks(see Figure 1.2). Once all items of an order have been picked, it must be
delivered at its pre-speciﬁed drop-off location (one of the locations marked
with 1-7 in Figure 3.1).
The set of all ordered items is indicated by G. Every item g ∈ G ischaracterized by three attributes: pick location; the customer order towhich it belongs to; and the number of the units that need to be picked.
The order to which item g belongs is indicated by o(g) ∈ O. The itemsof each order must be picked in a pre-speciﬁed sequence. The sequencenumber of item g ∈ G of order o(g) ∈ O is indicated by pg . The distance
between items g and k ∈ G is equivalent with the distance between thecorresponding pick locations. This distance is denoted by dg,k .
The ﬁrst item of an order corresponds with picking up an empty rollcontainer at the empty roll container depot. The last item of an ordercorresponds with the delivery of the order at its drop-off location. The
sets of ﬁrst and last items are denoted by Gf and Gl , respectively.
Note that the model described in this section is only meant to explain
the combined batching and precedence-constrained routing problems indetail. For computational purposes other models and techniques that aredescribed later on in this chapter are used.
Model

In order to model the batching and routing problem, (i) the orders

have to assigned to the batches, and (ii) the sequence of the items of eachorder needs to be kept track of. To that end, a set P of potential batches is
introduced. The number of potential batches should be sufﬁcient to carryout all orders. For example, |P| = |O| is certainly enough (but often much
too large), since this corresponds to the situation that each order is pickedon its own. Each batch can contain up to N bins. The set of bins is denotedby B. The batch to which bin b ∈ B belongs is denoted by r(b) ∈ P.
For each order o ∈ O and bin b ∈ B a binary decision variable Yo,b is introduced to indicate whether or not order o is picked in bin b. Furthermore,
let Xg,k,r be a binary decision variable indicating whether batch r ∈ P
goes directly from item g ∈ G to item k ∈ G. Finally, let Ug ∈ N indicate
the relative position of item g ∈ G in the tour that picks the order of itemg. The relative position is the integer position of g in the tour.

37Order Batching when Items Have to Be Picked in a Strict Sequence

Now the model can be described as follows:
 

min

(3.1)

Yo,b = 1∀o ∈ O(3.2)

Yo,b ≤ 1∀b ∈ B(3.3)

r∈Pg∈G k∈G\{g}subject toXg,k,r

dg,k

b∈B



o∈O



Yo(g),b∀g ∈ G \ Gl ∀r ∈ P(3.4)

Yo(k),b∀k ∈ G \ Gf ∀r ∈ P(3.5)∀g, k ∈ G : g = k

(3.6)∀g, k ∈ G : o(g) = o(k) ∧ pg > pk

(3.7)

b∈B:r(b)=rk∈G\{g}Xg,k,r =

Xg,k,r =g∈G\{k}b∈B:r(b)=r

Ug − Uk + |G| ·Xg,k,r ≤ |G| − 1r∈R

U g − U k ≥ pg − phYo,b ∈ {0, 1}∀o ∈ O ∀b ∈ B

Xg,k,r ∈ {0, 1}
Ug ∈ N∀g, k ∈ G : g = k ∀r ∈ P∀g ∈ G(3.8)(3.9)(3.10)

According to (3.2), each order is assigned to exactly one bin, and (3.3)
states that each bin covers at most one order. Thus the orders are correctly
assigned to the bins.
According to (3.4), each location related to an item (except a drop-off)is exited so that the corresponding order has been assigned to one of thebins in the batch. Similarly, constraints (3.5) state that each item, except
a ﬁrst location of an order, is entered so that the corresponding order isassigned to one of the bins of the batch. It follows that all items of anorder are picked by the same batch.
Furthermore, constraints (3.6) are based on the TSP subtour elimination constraints of (Miller et al., 1960). According to (3.6), the values Ug
represent the (relative) positions of the items in each batch: if a batch
goes directly from item g to item k, then (3.6) implies Uk ≥ Ug + 1. The
values Ug need not be consecutive integers within a batch. Anyway, for
two items g and k of the same batch, Ug < Uk stands if and only if item
g is picked before item k. It also follows that the sequence of items of abatch does not contain subtours.
Furthermore, constraints (3.7) state that two items belonging to the38Order Batching when Items Have to Be Picked in a Strict Sequencesame order are picked by the same batch in such a way that their (relative) positions in the batch are at least the a priori speciﬁed difference.
Thus Ug < Uk if and only if pg < pk . It follows that the precedence constraints of the items of each order are satisﬁed.
Constraints (3.8), (3.9) and (3.10) specify the binary and non-negative
character of the decision variables.
The objective function (3.1) represents the total routing costs of thebatches. Thus the optimal solution of (3.1)–(3.10) is the solution withminimum total routing costs.

3.2

Solution Approach

The algorithm presented in this section, Precedence-Constrained Estimated Savings based batching, or PCES for short, incorporates two distinct sub-algorithms:
1. an optimal A∗ -algorithm for the routing;
2. a simulated annealing-based combinatorial search algorithm, SAREMIX, for searching the space of possible batches.
Figure 3.2 shows the ﬂow diagram of PCES. The initialization is performed by inputting a selected set of customer orders O, in step 1, usuallymore than can be processed at once. In step 2, the next set of R customer
orders is selected for batching. Then batches with two customer orders areconstructed and routed in steps 3 and 4. In step 5, estimates of savings forbatches larger than two are calculated. To ﬁnd the best batch allocation,
a combinatorial search is performed in the space of batches with real andestimated savings in step 6. Estimated savings are calculated as neededin steps 5 and 6. In step 7, optimal routes are found for those batches thatmaximize the savings. Steps 5 and 6 are repeated for a ﬁxed number of
iterations to ﬁnd a good solution for the order-batching problem. If somecustomer orders are left in O, the algorithm returns to step 2, otherwise
it terminates with its solution in step 8.

39Order Batching when Items Have to Be Picked in a Strict Sequence

Figure 3.2. Solution approach for solving the batching problem with estimated savings
with PCES.

3.2.1 Estimating Savings from Batching More than Two
Customer Orders
The main reason for batching customer orders is to generate savings intravel distance compared to picking these customer orders independently.
If computationally feasible, all possible routes can explicitly and exhaustively be constructed from all possible customer order combinations anda combinatorial search performed in the exhaustive space, thus minimizing the travel distance. The number of possible combinations C up to
the maximum batch size N of all R customer orders can be calculated by40Order Batching when Items Have to Be Picked in a Strict Sequence

Equation (6.9), below.

C=

N


⎛
⎝

⎞
R

⎠

(3.11)ii=1

subject to R ≥ N and C, R, N ∈ N.
Especially with larger batch sizes (N ≥ 3), explicitly evaluating all possible batch combinations quickly becomes intractable due to the large
number of combinations and the complexity of the routing. This is countered by estimating the savings made for batches of more than two customer orders.
Instead of performing a combinatorial search in the explicitly evaluatedsearch space of all batch travel distances, only a part of the space consisting of batches each with real or estimated savings is searched. Realsavings are found by constructing routes for batches of size two and subtracting their travel distances from the respective single-order travel distances. Estimates for savings of batches of size three or larger are calculated on the basis of the real travel distances of batches of sizes one
and two. With PCES based batching, in addition to routing all batches ofsize two and one, only those batches with three or more customer orders
which are found using SA-REMIX need to be routed. Compared to routing
all batches with three or more orders, the number of calls to the routing
algorithm for batches of size three or more is reduced from

C3+ =

N

i=3

⎛
⎝

⎞
R

⎠

(3.12)

ito R/N (plus a possible remainder if N > 3).
In the following the way to estimate the travel savings for batches ofthree and more orders is explained. Let i, j and k be three orders that canbe batched together. If by batching order i with j, some travel distance issaved when compared to picking them independently, batching is beneﬁcial. The savings resulting from batching the two customer orders i and j
can be calculated by:
Si,j = ci + cj − ci,j(3.13)

where ci,j is the batched travel distance of customer orders i and j, ci
and cj the travel distances for picking customer orders i and j independently. The travel distances in Equation (3.13) are the real travel dis-41

Order Batching when Items Have to Be Picked in a Strict Sequencetances which have been evaluated with the routing algorithm presented
in Section 3.2.2.
Likewise, let batching i with k and j with k be beneﬁcial. This would
then imply that batching all three of these customer orders into a batchof size three would also be beneﬁcial. If these three customer orders arebatched, the savings can be calculated by:
Si,j,k = ci + cj + ck − ci,j,k(3.14)

Expanding from this, if i, j and k are batched, the savings can be estimatedby:S̄i,j,k = Si,j + Si,k + Sj,k(3.15)

Double counting can be eliminated by dividing (3.15) by two. However,as these three order batch savings estimates are used in the same search
space with the real two order batch savings, an added emphasis is given
to the larger batches by not eliminating the double counting. By combining (3.13) and (3.15), the estimated savings of a size three batch can becalculated:S̄i,j,k = 2ci + 2cj + 2ck − ci,j − ci,k − cj,k(3.16)Let V be a set of n arbitrary customer orders {1...n} and cn be the travel
distance of customer order n. Using this notation, if a batch is formed
from the customer orders in V, the savings can be found by:SV =

ci − cV(3.17)

i∈V

Using the same notation as above, the total savings resulting from batching N customer orders can be estimated:
S̄V =

N
−1


(

N


(ci + cj − ci,j ))(3.18)

i=1 j=i+1

Using (3.18), estimates can be calculated for all possible combinations up
to a given batch size. To use this equation, only the batch tour distances ofbatches of size one and two are needed. The real savings for batches of size
two are calculated, while for size three and more they are estimated. For
batches of size one, the savings value is set to zero. In practice, estimates
are calculated as they are needed, i.e., in steps 5 and 6 of Figure 3.2.

42

Order Batching when Items Have to Be Picked in a Strict Sequence

3.2.2

Pick Route Construction

The A∗ -algorithm introduced by Hart et al. (1968) is used to ﬁnd the optimal pick route. This algorithm’s states are organised in a lattice, as
states can be reached through multiple paths. The A∗ -algorithm’s stateis deﬁned by the vector m = (m0 , m1 , ..., mN ), similar to Psaraftis (1980b)
and Psaraftis (1983). m0 ∈ {1, ..., N } deﬁnes the customer order in which
the last product was placed. mi , i ∈ {1, ..., N }, is the number of locations
visited for customer order i (including the delivery location).
An example is shown in Figure 3.3. The algorithm’s state starts withempty bins, and the picker at the depot. The state in the starting location
is m1 = (m0 , 0, ..., 0), where m0 is an arbitrary customer order. To pick
an item to customer order i, mi is incremented by one and set m0 = i.
Keeping track of m0 is important, as the decision of which of the orders
to pick next depends on the location of the storage slot of the item thatwas last picked. Otherwise there would be up to N possible locations
where the picker could be when calculating distances to the next possible
pick locations instead of just one. Figure 3.4 gives an example of the state
space for two customer orders with two items to be picked in each of them.
The arcs between two adjacent states represent possible state transitions
by moving to the next location and picking the item required. The A∗ algorithm searches the lattice by using the exact travel distance to reacha state plus an estimate for the distance to reach the ﬁnal state, where
all customer orders have been picked and delivered and the picker is at
the depot. As long as the estimate underestimates the remaining travel
distance, the algorithm is guaranteed to ﬁnd the optimal solution as itprogresses through the lattice (Hart et al., 1968). Psaraftis (1980b) uses
a dynamic programming approach with no estimated distances to prune
the search space. He also considers additional constraints which are notused here, such as customer preferences between waiting and riding, andtime windows (Psaraftis, 1983).
The shortest total travel distance to reach a state m for a given set ofcustomer orders V is denoted by d(m, V) . For each state and the selected
set of customer orders, an estimate for the remaining travel distance from¯
state m to the depot and add it to the state cost is calculated by d(m,V).
The single customer order with the longest remaining travel distance to¯
the goal is found and d(m,V) is set to equal that distance. Let d¯i (m, Vi ) be
the travel distance to goal by only visiting the locations of customer order43Order Batching when Items Have to Be Picked in a Strict Sequence

Figure 3.3. A truck picking three customer orders into three separate bins. Items inbold have already been picked; the underlined one denotes the last pick anditalicized ones are possibilities for the next transition, representing state(2,2,1,3).

Figure 3.4. Left: a Hasse diagram of the state space formed by two customer orders
with two items to be picked into each of them. Inside the circles, the statename is at the top, the state vector is at the bottom. The ﬁrst member ofthe state vector refers to the customer order where the last pick occurred,
the following two refer to the customer order location indices. The highest
location index (3 in this case) represents the delivery location for the given
customer order. Right: possible transitions and distances from the state m9in the Hasse diagram.

44

Order Batching when Items Have to Be Picked in a Strict Sequencei, i ∈ {1, ..., N }. The estimate is calculated by:
¯
d(m,
V) =

max

i∈{1,...,N }d¯i (m, Vi ).

(3.19)The estimate of the total distance to the goal state is calculated from¯ē(m, V) = d(m, V) + d(m,V). Doing this prevents some of the unnecessary
exploration of states that are not part of the optimal tour.
As an example, if there are two customer orders with two products to bepicked in each and a delivery location for both orders, a lattice is formedby all the possible states mi , ∀i ∈ {1, ..., 26}. In Figure 3.4, all states
are enumerated with m1 being the starting state and m26 being the end
state. Arrows indicate the reachability of states. The number of statesin the lattice and the time complexity for the algorithm (Psaraftis, 1980a)
can be found from (3.20) and (3.21), respectively. Note that in (3.21) the
product is preceded by N 2 , and not by N as in (3.20). This is due to thefact that all states in the lattice have at most N incoming arcs (see Figure
3.4).
N

M = O(N|Vi + 1|),

(3.20)|Vi + 1|).

(3.21)i=1
N

V = O(N 2i=1where the sum term removes all states that are infeasible such as [1, 0, 1].
|Vi | is incremented with one as the delivery location is included.
In Figure 3.4 on the right side, the possible transitions and distancesfrom state m9 are shown. Let the total distance to that state be a9 . To
calculate the heuristic distance to goal, all remaining items can either be
added to bin 1 or 2, but not to both. The heuristic distances are thusd¯1 (m9 , V) = 1 (as only one item can be added to order 1) and d¯2 (m9 , V) =2 + 4 = 6 (as two items can be added to order 2). d¯2 (m9 , V) is chosen forthe heuristic as per (3.19). Then the estimate of the total cost of the goalstate is ē(m9 , V) = a9 + d¯2 (m9 , V). The estimate for the total distance togoal is very optimistic. If a better estimate were to be used, the runningtimes of the algorithm might decrease.
The A∗ -algorithm is used to route batches of two customer orders, aswell as the ﬁnal batches after batching has been ﬁnished, in steps 4 and
7 of the PCES algorithm, as explained in Figure 3.2, respectively.

45

Order Batching when Items Have to Be Picked in a Strict Sequence

3.2.3 Batching
After all savings and travel distances for batches of size two have been
determined, orders have have to be batched in batches of size up to N .
The batching algorithm presented in this section, SA-REMIX, is a Large
Neighborhood Search algorithm (see Section 1.4.5) consisting of a construction phase followed by improvements embedded in a simulated annealing scheme to avoid local optima. The decision to use a simulated
annealing metaheuristic in the algorithm was made based on it being
previously incorporated into Large Neighborhood Search-type methods in
Ropke and Pisinger (2006a) and Pisinger and Ropke (2010). C&W(i) isused to form an initial solution. This starting solution gives an upperbound for the total travel distance. Next, a REMIX (a random destroyrandom repair) procedure is applied.
REMIX

This procedure randomly selects a ﬁxed number of batches fromthe ones that are currently selected (steps 5 and 6 in Figure 3.2). The
orders contained in the batches are mixed and randomly sampled again
to form a number of new batches. The REMIX procedure is embedded
in a simulated annealing routine. The size of the new batches is limited
to N . The total travel distance savings of these new batches (Snew ) are
estimated and compared with those of the original batches (Sorg ). Let
Q ∼ U [0, 1] be a uniformly distributed random variable and Tn be thecurrent temperature, where n is the current iteration. New batches are
accepted if one of the following two criteria is met:
1. if their savings are larger than those of the original batches, i.e., Snew >Sorg ;
2. or, otherwise, if Q < exp((Snew − Sorg )/Tn ).
The initial temperature T0 is set to equal the highest savings value of thebatches with two customer orders. Updates to Tn are done using Lundy’s(Lundy and Mees, 1986) temperature schedule: Tn+1 = Tn /(1 + Tn /n).

3.3

Validation and Results

In this section, the results of batching with the method presented in Section 3.2.3 are compared both with company results and with optimal and46

Order Batching when Items Have to Be Picked in a Strict Sequence

heuristic batching, while using the A∗ -routing. The company data is used
as input for the experiments. Data consisting of voice-picking logs of almost 35,000 batches (approximately 100,000 customer orders) was analyzed, each consisting of up to three customer orders and an approximate
total of 1.8 million log entries (i.e. pick locations to be visited), spanningthe ﬁrst three months of 2011. Each log entry includes a slot address, picktime, batch number, customer order number and picker number. Also thesequence in which each customer order should be picked is given. Currently the batches are created by grouping the customer orders per truck
delivery route in a FCFS sequence. The pick routes are constructed using
the S-shape heuristic, preserving single-directional aisle travel and strict
item precedence.
The experiments are carried out with 1,536 customer orders sampledfrom the data set (about one day’s work). The customer orders of each day
are grouped in waves. Each wave represents a number of truck delivery
routes with scheduled departure times within a certain time interval. To
make a meaningful comparison between different algorithms, the numberof customer orders per wave, or wave size R, is set to vary between 12 and150. As most of the orders arrive early before the morning shift starts,
it is possible to do order batching for much larger instances, especially iftruck departure times overlap. Due to the relatively long computationaltime of the routing, such instances are impractical to solve with most of
the tested algorithms from the literature. The results show, that after
the wave size of 126 the differences in relative performance between the
different algorithms do not change much, only the computational time
grows.
All experiments were run single-threaded on an Intel(R) Core(TM) i7
920 CPU running at 2.67 GHz and with 12 GB of memory. Due to the
exponential complexity of most algorithms that solve the order batchingproblem for a maximum batch size of three or more orders, the time tosolve a given problem instance is considered in addition to the solution
quality. The space of possible batch combinations when maximum batch
size is three or more increases exponentially with the wavesize R. Thus
the number of simulated annealing iterations run with PCES is increased
with wavesize. A good number of iterations for the SA-routine that scaleswith both wave and batch sizes was found to be J = 500, 000 + (2N − 5)R2 .
This was determined ad-hoc fashion. Other amounts should work as well,
but this amount produced consistently good results and proved scalable.

47Order Batching when Items Have to Be Picked in a Strict Sequence

Case Description

The warehouse is operated by a retail store chain sell-

ing non-food products and supplying 148 stores mostly on a daily basis.
Orders arrive online and can be processed in ﬁve or six waves per day of R
customer orders each. Orders are picked manually from pallet racks into
customer orders using a picker-to-parts order picking system (De Koster
et al., 2007). Pickers have to travel long distances as the warehouse has
57 aisles of 104 meters each and three cross-aisles, see Figure 3.1. The total warehouse area is 75,000 m2 , the average pick tour is 880 meters, anda warehouse picker travels an average of 7300 meters a day. The company
sells and stores over 50,000 SKUs. Each warehouse aisle has 52 pick slots,
which are divided into four subslots and two storage levels. This amounts
to a total of 23,296 different pick locations in the warehouse. These locations are condensed into 2912 locations and corresponding items. The
distance between two locations is measured using Manhattan distances,taking into account that pickers can only travel along the aisles in one
direction, for safety reasons, and to avoid congestion. If the next pick isfrom a slot which is in the same aisle but in opposite travel direction, thetravel distance is calculated respecting these trafﬁc rules.
The warehouse uses a family grouping strategy to allocate items to storage slots, similar to the stores. Such a storage strategy is common in
many retail warehouses. Efﬁciency losses in the warehouse by not usingclass-based storage are traded-off against efﬁciency gains in the storesby faster customer order unloading and shelf replenishment (De Koster
et al., 2007). The customer orders must be ﬁlled in a predeﬁned sequence,
primarily due to the family grouping. Each item has a unique pick location in the warehouse. The contents and pick sequences for each of
the customer orders are generated a priori by a Warehouse Management
System (WMS). Every pick route contains up to three customer orders (abatch), currently all for the same store, in one tour of the warehouse (seeFigure 3.1), but this might change, depending on potential travel time savings. A picker drives an order picking truck (see Figure 1.2) and follows
a tour determined by the on-board customer orders, directed by a truckmounted RF terminal. In addition to the truck terminals, the company
uses a pick-by-voice system, which is not used for guiding the pickers, but
only for pick conﬁrmation (De Koster et al., 2007). Pickers can overtake
each other. The warehouse layout can be seen in Figure 3.1. The delivery
location of each completed customer order is known in advance: it is theone closest to the dock door from which the outbound truck will depart to48

Order Batching when Items Have to Be Picked in a Strict Sequence

the store.
Each order picking tour related to a single batch starts from the depot,
see Figure 3.1. Before returning to the depot, all relevant pick locations
and drop-off points must be visited. There are seven drop-off locations in
the warehouse.
Validation of the Estimation Method

To speed up batch evaluation, allbatches of three customer orders or more are not explicitly routed in PCES
but estimated using (3.16) or (3.18). To validate this approach, a set of
10,000 batches (P) of three customer orders is sampled from the previously sampled data set of 1,536 customer orders. The 10,000 batches are
all routed, and real and estimated savings are calculated for all r ∈ Pusing (3.14) and (3.16), respectively. Thus when calculating real savings, explicit routing is done for batches size three or more as well. Allbatches with negative estimated savings are discarded, as they would
never be chosen by SA-REMIX, resulting in just over 5000 batches or
the set P + ⊂ P. Both sets of savings, i.e., real and estimated (Pv+ and
Pe+ , respectively), are then ordered in descending order to establish realand estimated integer batch ranks vvr and ver , ∀r ∈ P + . For validation of
the estimate of a batch r, a proportional batch rank error calculated byr = |ver − vvr |/|P + | is used.
Figure 3.5 shows a moving average over 20 values of r , ∀r ∈ P + , sorted
by real batch rank (with 1 being the batch with the highest savings). For
the best batch, the average is over batches ranked 1,..,20, for the second
best, it is over 2,...,21 and so on. In SA-REMIX, the interest lies in batcheswith the highest savings, i.e., ﬁnding the best feasible combination of thebest ranked batches. Figure 3.5 shows that the savings estimates for thebest batches are much better than for the rest of the batches. This validates the performance of PCES as seen later in the results.
Comparison to Batching Heuristics from Literature

The solution quality

and the computation time of PCES are both compared to the following
heuristics for maximum batch sizes three and four: C&W(i) and C&W(ii)from De Koster et al. (1999); the VNS-batching algorithm by AlbaredaSambola et al. (2009); ABHC by Henn and Wäscher (2012) which uses either C&W(i) (a new variant of the ABHC used in this chapter) or C&W(ii)(as in Henn and Wäscher (2012)) as the initial solution; and the generally
best performing seed algorithm combination from De Koster et al. (1999)(seed selection by ﬁnding the order with the longest travel time, and order addition by the order that minimizes the sum of the distances from49

Order Batching when Items Have to Be Picked in a Strict Sequence

Figure 3.5. Error % in comparing batch ranks, ordered by the batch with the highest
total rank for 10,000 batches of three customer orders sampled from a set of
1,536 customer orders. Only 4,900 are shown here so that small ranks canbe distinguished.

the closest item in the seed to every item in the order). A 240 minute time
limit is imposed for each wave size.
The results are shown in Figure 3.6. Figure 3.6a shows the percentage improvement compared to C&W(ii) in travel distance for wave sizes
between 12-150, for all the other heuristics for a maximum batch size of
three customer orders. C&W(i) is the weakest heuristic and always inferior to C&W(ii), as has been shown in the literature previously, see, e.g.,
De Koster et al. (1999) and Albareda-Sambola et al. (2009). However,
C&W(i) is computationally light. The seed selection heuristic did not perform well, so it was left out of the results. Both the VNS and ABHC
base their good results on detailed exploration of the combinatorial batch
space, thus requiring a lot of batch evaluations. As the space gets larger
with larger wave sizes, the number of batch evaluations increases exponentially. Initially, up to wave size 30, the VNS heuristic comes out on
top — after this, the time limit for solving the problem kicks in, and thesolution quality degrades rapidly, see Figure 3.6a. This is partly due tothe fact that no initial solution is used - the search starts from each order
being in a separate batch.
A similar effect imposed by the time limit and the computationally expensive batch evaluation can be seen with ABHC. Initially, it is the second
best algorithm, but within the given time limit it can only make some im-50Order Batching when Items Have to Be Picked in a Strict Sequence

provements to the initial solution found with C&W(ii). Henn and Wäscher(2012) suggested that they got comparable results with initial solutions
found with a FCFS-algorithm. However, during the tests it was foundthat using C&W(ii) to obtain a decent initial solution requires less computation time and produces better results. In general, PCES has better solution quality than C&W(ii) with less computation time (see Figure3.6b). Also, for medium to large wavesizes, its solution quality is better
than those of all the other methods, while simulataneously being much
faster.
In Figure 3.7a the solution qualities of PCES and ABHC (initial solution with C&W(i)) are compared to C&W(i) for batches of four customer
orders. The number of possible batch combinations is much larger thanfor batches of size three (see (6.9)). The routing is also computationallymuch more expensive. C&W(ii), VNS and thus also ABHC with C&W(ii)all fail to reach a proper solution result withing the time limit and arethus left out from the results. Using PCES yields a mean improvementof 2.3% compared to the batching of C&W(i), while ABHC is able to improve the C&W(i) initial solution by 0.01% on average, mostly because of
the timelimit. Computational times of C&W(i) and PCES are similar forsmall wave sizes. PCES nears the timelimit as the number of iterations of
SA-REMIX is increased quadratically with wave size, while ABHC always
hits the time limit.
In both instances of three and four customer orders, PCES exhibits similar behaviour. For small wavesizes (up to 24) the error in the batch savings
estimate dominates, and a good result is not reached reliably. However,recall that the estimate (Equations (3.16) and (3.18)) is quite good in estimating batches with high savings (see Figure 3.5), so as the number ofsuch batches increases exponentially with both wave and batch size, theprobability of ﬁnding a good solution increases with it.
Furthermore, it was found that the solution time is proportional to thenumber of batch evaluations of maximum batch size (three or four in thisstudy).
Comparison to optimal batching

It is possible to compare the results ofthe batching procedure presented in this chapter, SA-REMIX, with optimal batching for up to 27 customer orders and batches of size three. The
optimal routing algorithm is based on the A∗ -routing method presented
earlier in this chapter, while the order batching part uses dynamic programming, with all possible customer order combinations as state space.

51

Order Batching when Items Have to Be Picked in a Strict Sequence

For the optimal batching the wave size R is varied between 12 and 27,
in steps of three customer orders. Figure 3.6a shows the results. For
the optimal algorithm (opt in Figure 3.6a, see also Section 2.2.1 for thealgorithm), the solution time increases exponentially with the wave size
— calculating the optimal solution for the whole dataset of 1,536 orders
with 27 customer orders per wave and N = 3 took over 23 hours. When
compared to optimal solutions of waves of sizes between 12 and 27, PCES
based batching gives results with less than a 1.2% error on average compared to the optimal solution for N = 3. Using PCES, when R increases,
the time to solve the problem instance decreases. This is due to a reduction in calls to SA-REMIX. Note also the differences in computation timein Figure 3.6b.
Effect of wave and batch size

For batches of one customer order, wave

size does not matter as each customer order is picked separately. For
the batches of size two, travel distances are calculated with the package
CPLEX 12.4 using the order-batching model presented by Gademann and
Van de Velde (2005).
Figure 3.8 shows the beneﬁts of different wave and batch sizes. Having more customer orders per wave reduces the total travel distance. For
waves larger than 150 customer orders, there are no major changes inthe different algoritms’ relative performance. When comparing optimalbatching for two customer orders to three customer orders, with a wave
size of 27 customer orders, two customer order batches appear to lead to19.4% longer travel distances. If the batch size is increased to four, there
are still considerable savings to be had in travel distance when comparedto those of batches of three customer orders. However, routing for thesefour-order instances is much more time consuming, and in practice thisrules out the computation of an optimal baseline, as well as the application of most of the tested order batching heuristics from literature, apartfrom C&W(i).
Re-routing and rebatching of original batches

Now the effect of PCES

batching is compared with the company’s current method, including the
mandatory travel directions. Since the company’s routing method is not
optimal, an additional comparison is made to see how total travel distanceis affected by just rerouting the original batches using the A*-routing
method presented in Section 3.2.2. The results can be found in Table
3.1. By rerouting the original batches with a wave size of 126 customer
orders, 3.1% or 980 km is saved, while rebatching and rerouting together

52Order Batching when Items Have to Be Picked in a Strict Sequence

(a)

(b)
Figure 3.6. Comparison of PCES to other heuristics and an optimal solution for batch
sizes of 3. The top ﬁgure shows the quality of the solution when compared to
the C&W(ii), while the bottom one shows the computation time for each wave
size for solving an instance of 1,536 orders.

53Order Batching when Items Have to Be Picked in a Strict Sequence

(a)

(b)
Figure 3.7. Comparison of PCES to other heuristics and an optimal solution for batch
sizes of 4. The top ﬁgure shows the quality of the solution when compared to
the C&W(i), while the bottom one shows the computation time for each wave
size for solving an instance of 1,536 orders.

54Order Batching when Items Have to Be Picked in a Strict Sequence

Figure 3.8. Percent of travel distance saved for optimal (solid lines) and PCES based
batching (dashed lines) as a function of the wave size R in customer orders,for batches of size N = 2, N = 3, and N = 4, in comparison to N = 1. The
mean error of PCES batching compared to optimal is 1.2% when N = 3.
Table 3.1. Results from rerouting and combined rerouting and rebatching of the original
batches.
Month
Jan
Feb
Mar
Total

ORG [km]1 RR [km]2 RRS [km]3 RRS %4 RB [km]5 RBS [km]6 RBS %7
9698
9381
317
3.3%
8168
1530
15.8%
10189
9722
467
4.6%
8437
1752
17.2%
12060
11873
187
1.6%
10313
1748
14.5%
31947
30967
980
3.1%
26918
5029
15.7%

1 ORG is the travel distance resulting from the original batching and routing.
2 RR is the rerouted travel distance with the original batching.
3 RRS is the savings using the rerouted travel distance with the original batching.
4 RRS % is the percentage of saved travel distance using the rerouted travel distances with theoriginal batching.
5 RB is the rerouted travel distance with PCES batching.
6 RBS is the savings using the rerouted travel distance with PCES batching.
7 RBS % is the percentage saved in travel distance with rerouting and PCES batching.

give savings of 15.7% or 5029 km.

3.4

Discussion

This chapter presents PCES batching based on optimal precedenceconstrained routing in combination with a simulated annealing algorithm
which uses estimated travel distance savings to group orders in batches.
This chapter’s contribution are in developing a generic method for solving
the difﬁcult problem of batching orders with strict sequence constraints.
The method uses an estimate for savings generated by combining orders, which appears to perform remarkably well. The estimate uses exact
savings from forming batches of size two to generate savings for larger

55Order Batching when Items Have to Be Picked in a Strict Sequencebatches. By using the estimate, the number of routings to be calculated
reduces dramatically, particularly for larger batches.
The PCES algorithm is able to generate signiﬁcant savings of up to
15.7% compared to the original batching, or over 5000 km for a threemonth period. These savings can be achieved without having to makeany changes to the warehouse layout or customer order composition. If
the cost of order picking is 55% of total warehouse expenses (Drury, 1988)
and 50 % (Tompkins et al., 2003) of this is due to travelling between picks,
the savings on the total operating costs of the sample warehouse would be
4.3%.
The core idea behind PCES is to minimize the number of computationally complex batch evaluations. This makes it very competitive with otherstate-of-the-art batching heuristics, especially in solution time, but alsoin quality. Compared to other algorithms, which use a set of deterministic neighbourhood heuristics, the randomness inherent to the REMIX
heuristic allows it to be ﬂexible and perform well. For batches of threecustomer orders, PCES generally ﬁnds better solutions than the C&W(ii)algorithm, while its computational complexity is between C&W(i) andC&W(ii), and beating out other algorithms from literature for medium
to large sized waves. When comparing performance to C&W(i), previous
studies (Albareda-Sambola et al., 2009, Henn and Wäscher, 2012) haveachieved proportionally higher savings than presented in this chapter.
This is probably due to a larger maximum batch size and the fact that
real orders are used in this study, whose composition is already somewhat
optimized. However, as can be seen from the results in this chapter, PCES
compares favorably to both the VNS and ABHC algorithms presented in
Albareda-Sambola et al. (2009) and Henn and Wäscher (2012). For fourorder batches, PCES is peerless in solution quality with a 2.0% improvement over the next best algorithm. However, with small wave sizes of
up to 30 customer orders, one would get better results using either of the
VNS or ABHC heuristics, due to the inconsistent behaviour of PCES inthat range (see Figure 3.6). This inconsistency is most likely due to factthat the probability for the estimate to be accurate is much smaller forsmall waves. As the number of good batches increases exponentially withwave size, so does the probability of the estimate being accurate.
When compared to optimal batching for small wave sizes of up to 27
customer orders and a maximum batch size of three, PCES performs decently, with a 1.2% mean error compared to the optimal solution. Savings

56Order Batching when Items Have to Be Picked in a Strict Sequencein travel distance that are comparable to optimal batching and routing
can be achieved with much less computation time. When compared to justpicking one bin at a time, Figure 3.8 shows that increasing the batch size
can result in substantial travel distance savings. Depending on the wave
size, for batches of two orders, 30% to 40% is saved, while for batchesof size three and four the savings range from 40% to 54% and 45% to61%, respectively. Larger wave sizes seem to have a considerable effect on
the total travel distance for up to about 100 customer orders. Wave sizes
larger than this do not seem to affect solution quality very much.
Another beneﬁt of the PCES batching method is that it allows us to evaluate the effect of different batch sizes in any warehouse with precedenceconstrained routing. In case a different routing scheme is followed, it israther straightforward to replace the routing algorithm presented in section 3.2.2 with this other method.
In practice, there may be further opportunities to reduce travel distanceeven more. New storage policies were not considered (see Dekker et al.
(2004)), which might lead to additional savings. In addition, pickers’
skills appear to vary signiﬁcantly. Assigning the right pickers to the rightbatches (orders and products) will also impact total time. This is dealt
with in Chapters 5 and 6.

57

Order Batching when Items Have to Be Picked in a Strict Sequence58

4. Batching when Routing with a
Generic TSP Algorithm

This chapter explores the performance of the batching part of PCES, the
SA-REMIX algorithm, (see Chapter 3) when the precedence constrained
routing is relaxed. The A* routing is replaced with a general-purpose
TSP-solver, the Lin-Kernighan-Helsgaun algorithm (LKH) (Helsgaun,2000). Theys et al. (2010) show that LKH is applicable to order picking
problems irrespective of the warehouse layout, and that it provides goodsolutions (with a gap of 0.1%). SA-REMIX uses estimated savings forbatches composed of three customer orders or more, as shown in Section
3.2.1. The results presented in Chapter 3 show that the algorithm performs very well when compared to other heuristics and optimal solutions.
Central to this performance are the estimated savings, which are accurate for the most important part of the batches — those with high savings
values (see Section 3.3). This chapter deals with the followin researchquestion: how will SA-REMIX perform if the precedence constraints arerelaxed and the A* algorithm is replaced with another routing algorithm?
The results in this chapter show that SA-REMIX batching can be applied
to any picker-to-parts order picking process in any warehouse not constrained by a particular layout or a speciﬁc routing method.
This chapter is organized as follows. First, the (non-precedence constrained) joint routing and batching problem is presented. Second, theexperimental setup is detailed, followed by the results. The chapter concludes with a summary. LKH can be downloaded from Keld Helsgaun’s
website: http://akira.ruc.dk/~keld/research/LKH/.

4.1

Problem Description

In this section, the model presented in Section 3.1 is modiﬁed by eliminating the strict internal precedences. Furthermore, it is now assumed that59

Batching when Routing with a Generic TSP Algorithm

no set drop-off locations exist, but all orders are dropped off at the depot(see Figure 3.1). This is done because LKH does not handle precedences.
The model presented in Section 3.1 is modiﬁed as follows. Equation (3.7)
is left out, and sets of ﬁrst of and last items Gf and Gl , used in (3.4) and(3.5), are not substracted from from the set of all items G. Gf existed because each tour began by visiting one of the ﬁrst items in each order, andthis condition no longer applies due to the lack of precedence constraints.
Gl is not used as it was composed of the drop-off locations, which are notused in this model. As with the previous model, each pick truck has anumber of bins, to which orders can be assigned to. The objective is tominimize total travel distance.
Model

Similarly as in Section 3.1, in order to model the batching androuting problem, all of the orders have to assigned to the batches. However, this time the sequence of items to be picked to each order is not
relevant. To that end, let P be a set of potential batches. The number ofpotential batches should be sufﬁcient to carry out all orders. Each batch
can contain up to N bins. The set of bins is denoted by B. The batch towhich bin b ∈ B belongs is denoted by r(b) ∈ P.
For each order o ∈ O and bin b ∈ B a binary decision variable Yo,b is introduced to indicate whether or not order o is picked in bin b. Furthermore,
let Xg,k,r be a binary decision variable indicating whether batch r ∈ P
goes directly from item g ∈ G to item k ∈ G. Finally, let Ug ∈ N indicate
the relative position of item g ∈ G in the tour that picks the order of itemg.

60

Batching when Routing with a Generic TSP Algorithm

 

min

(4.1)

Yo,b = 1∀o ∈ O(4.2)

Yo,b ≤ 1∀b ∈ B(4.3)

r∈P

g∈G k∈G\{g}subject toXg,k,r

dg,k

b∈B



o∈O


k∈G\{g}Xg,k,r =Yo(g),b∀g ∈ G ∀r ∈ P(4.4)

Yo(k),b∀k ∈ G ∀r ∈ P(4.5)

b∈B:r(b)=r

Xg,k,r =g∈G\{k}b∈B:r(b)=r

Ug − Uk + |G| ·Xg,k,r ≤ |G| − 1∀g, k ∈ G : g = k

(4.6)

r∈R

Yo,b ∈ {0, 1}∀o ∈ O ∀b ∈ B

Xg,k,r ∈ {0, 1}
Ug ∈ N∀g, k ∈ G : g = k ∀r ∈ P∀g ∈ G

(4.7)
(4.8)
(4.9)

Constraints (4.2) enforce the assignment of each order to exactly one
bin, and (4.3) states that each bin covers at most one order. This results
in a correct assignment of orders to bins.
Each location related to an item g ∈ G is exited so that the correspondingorder is assigned to one of the bins in the batch as in (4.4). Equation(4.5) enforces that each item is entered so that the corresponding orderis assigned to one of the bins in the batch. Thus, all items of an order
are picked by the same batch. Constraints (4.6) are based on the subtour
elimination constraints of (Miller et al., 1960). Constraints (4.7), (4.8) and(4.9) enforce non-negativity and binary qualities of the decision variables.

4.2

Results

Similarly as in Section 3.3, a set of 1536 orders is sampled from the data
set of the Finnish retailer. The customer orders are grouped in waves,with the wave size R varying between 12 and 150 orders. The maximum
batch size N is set to three. From (6.9), the total number of possible batch
combinations per wave ranges from 298 (12 orders per wave) to 562625(150 orders per wave). C&W(i) is used as the baseline algorithm for comparisons. It is compared to C&W(ii), SA-REMIX, the VNS by Albareda-

61

Batching when Routing with a Generic TSP Algorithm

Sambola et al. (2009), and an optimal solution for small instances. C&W(i)is used as the initial solution for SA-REMIX. The number of iterations J
used in the experiments is J = 500, 000 + R2 , as the number of possiblecombinations increases with wavesize. For the algorithms to process all
1536 orders, a maximum timelimit of 10h is imposed.
Solution quality

Figure 4.1a shows the solution quality of the algorithms.

As in Section 3.3, VNS is the algorithm with the best solution quality forwave sizes smaller than 60. From that point on, SA-REMIX provides
the best solutions, being on average 2.6% better than C&W(i). For waves
equal or larger or than 60 orders is the next best algorithm is VNS (up to
114 orders, after which the timelimit kicks in), followed by C&W(ii). When
compared to optimal batching and routing with LKH, the combination of
LKH and SA-REMIX provides solutions that are on average 1.2% worse,
which is in line with the results in Chapter 3.
Computational time

Due to the complexity of the algorithms that areused to solve joint picker routing and order batching problem, computational time was considered. Figure 4.1 shows the results for the tested
algorithms by wave size. The C&W(i) and C&W(ii) algorithms are thefaster ones, in that order, followed by SA-REMIX. This is a bit surprising,as previously in Chapter 3 the computational effort of SA-REMIX wasseen to be somewhere in between the two C&W-algorithms. Undoubtedly
the reason for this is the changing of the routing algorithm, at least indirectly. Moreover, the C++ implementation of the SA-REMIX algorithm
can be further optimized, but this is left as future work. The VNS exhibits
similar behavior as in Chapter 3, with steadily increasing computation
time. Eventually ﬁnding a good solution to the batching problem with
VNS becomes impossible due to the time constraint of 10h.

4.3

Discussion

The results show that using SA-REMIX, which is the batching part of
the PCES algorithm, can also be used to solve batching problems without
precedence constraints. Moreover, the algorithm provides good overallsolution quality and the best performance within the comparison groupwith large wave sizes.

62

Batching when Routing with a Generic TSP Algorithm(a) Solution quality per wave size for different order batching algorithms compared to C&W(i).

(b) Computational time per wave size for different order batching algorithms.
Figure 4.1. Comparison of PCES to other heuristics and an optimal solution for batch
sizes of 3. The top ﬁgure shows the quality of the solution when comparedto the C&W(i), while the bottom one shows the computational time for eachwave size.

63

Batching when Routing with a Generic TSP Algorithm

64

5. Forecasting Batch Execution Timefor Individual Pickers

Most warehouse management systems (WMS) store order picking logs,
which are captured at a very detailed level by advanced picking tools. In
this section, such log data are used to construct models to forecast thebatch execution time for individual pickers. Batch execution time may
depend on many factors, such as the details of the batch to be picked, butalso on behavioral factors, such as intrinsic and extrinsic motivation andability (Larco Martinelli, 2010). However, rather than explicitly including
behavioral factors, they are implicitly included by only considering thepast performance of each picker, as this is what can be found in the WMS
data. WMS data of a picked order in a batch typically contain: picker ID,
roll cage IDs in which the items are picked, drop-off locations of the rollcontainers, time-stamp of each order line (a part of customer order thathas to be picked, also called an item or pick line), slot address per line,
item IDs, number of units picked, etc. These data are used in a multilevel model, where the pickers form the "groups". Multilevel modeling isnaturally suited for distinguishing between-group (i.e., between-picker)
differences. It also scales well for groups of different amounts of input
data. A group’s random effects are allowed to deviate more from the ﬁxed
effects if it has a large amount of input data. As potential independent
variables, the following are selected: number of pick lines in a batch, total batch travel distance, total pick item mass, and volume in a batch, aswell as the mean pick height level at which items are picked during thepicking tour.
First, the dataset is brieﬂy described, followed by an explanation of theused regression method and the data cleaning. Finally, the regression results are shown. Cluster analysis is done on the picker models to illustrate
the differences between the pickers.

65Forecasting Batch Execution Time for Individual Pickers

5.1

Warehouse and Dataset Description

The sample warehouse in this chapter is a large picker-to-parts warehouse in Finland. Three months of extensive pick data from the warehouse’s pick-by-voice system were obtained and processed. The warehouse
has three cross-aisles, 57 aisles, seven drop-off locations, unidirectional
travel in the aisles, and a single depot storing empty roll containers (seeFigure 3.1 for the layout). Each customer order contains a maximum of 50
order lines, and is picked to a single roll container using a motorized truckthat can transport a maximum of three roll containers. Thus batches consist of a maximum of three customer orders.
The data-set allows the extraction of original customer orders andbatches and the identiﬁcation of the picker who did the work. The batch
execution time depends on: (1) picker ID, (2) total number of lines, (3) total
travel distance calculated from the distance between all sequentially visited locations in a tour [m], (4) total mass [kg], (5) mean pick level, where
1 is low and 2 is high, and (6) total volume [m3 ].
The models are then used in an ALNS algorithm (as explained in sections 6.1, 6.3, and 6.4) to assign orders and batches to pickers. The wholedataset consists of 37,841 batches handled by 229 pickers during a threemonth period.

5.2

Data Cleaning

To ﬁnd possible outliers, the data were preprocessed. The results are
summarized in Table 5.1.
Three different ways of picking were found: (i) multiple orders in a batch(pick tour), (ii) individual orders, and (iii) multiple orders with no apparent routing, and all lines with the same timestamp. The second and third
methods can occur in exceptional situations and represent non-standard
ways of working. The second occurs when picking heavy items to pallets
instead of roll containers (Pallet), and the third occurs when the warehouse manager signs off previously picked orders by a single press of abutton (Timestamp). The large majority of the batches (83%), representedby (i), is used in the regression modeling; (ii) and (iii) were not considered
due to their representing a different way of working.
In some cases, trucks seemed to exceed the maximum speed of 10 km/h.
Such batches are ﬁltered out (Speed). Some items could not be mixed

66Forecasting Batch Execution Time for Individual Pickers

PalletTimestamp

Speed

Mass

Long time

Many pickers

ZeroesPicker out

Remaining

Batches
Percentage

Total

Table 5.1. Summary of data cleaning. The cleaning categories are in the order they areapplied — thus Speed is ﬁrst and Picker out is last.

37841
100%

3257
8.6 %

3148
8.3 %

852
2.3%

553
0.1%

146
<0.1%

146
<0.1%

34
<0.1%

4070
10.8%

24669
65.2%

Table 5.2. Mean values and coefﬁcients of variation per batch before and after data
cleaning.
Before clean-upμCV
After clean-upμCVTime [min]LinesTravel [m]Mass [kg]Level

Vol [m3 ]

29.6
1.21

46.3
0.77

598.5
0.96

206.5
0.78

1.1
0.12

1.0
0.53

31.7
0.69

51.8
0.65

620.3
0.63

224.0
0.67

1.1
0.11

1.0
0.40

with other products, e.g., because they were highly fragile. The weightof such line items is artiﬁcially inﬂated to reach the threshold weight of
600 kg per order. Orders containing such lines were discarded (Mass). In
addition, batches that took too long to execute were also discarded. Thiscan be caused by breaks, shift changes in the picking process, or task
interruptions with other causes. Based on experience in warehouses withsimilar picking processes, the threshold was set at a maximum of 2 hours
net picking time per batch (Long time). Batches that were picked by morethan one picker were also excluded (Many pickers).
Neither the output nor the inputs were allowed to contain zero values— all such batches were omitted from consideration (Zeroes). Finally, toconstruct a reasonably large set of cross-validation data, it is required foreach picker included in the model to have at least 75 batches of input data
(80% of which are used for tuning the model) (Picker out).
After the data cleaning and the requirement on the number of data lines,
99 pickers out of 229 qualiﬁed. Table 5.1 shows a summary of the resultsof the data cleaning. After completing the data cleaning, 24,669 batchesremained to be used for the modeling of the pickers.
Table 5.2 lists the means and coefﬁcients of variation (CV) for the inputs
and outputs. Most of the CV values are smaller than one: this suggeststhat many of the tasks (batches) that the pickers have done are quite similar, or at least comparable. In the non-cleaned data, batch execution time
can vary quite a bit, most probably due to breaks, etc., that are ﬁlteredout in the data cleaning. Batch parameters Lines, Travel and Mass seemto vary the most, while Level seems to vary very little: most of the picking
is done from the low level.

67Forecasting Batch Execution Time for Individual Pickers

5.3

Multilevel Modeling

Let w ∈ W denote the groups (the pickers) and r ∈ R are indices to thedata (the batches). Furthermore, let vector β w contain the model parameters for group w and let tr,w be the output for a forecasting model, which
forecasts total batch processing time, for some r and w. Using the notationof Bryk and Raudenbush (1992), a linear multilevel forecasting model canbe formulated as in equations (5.1).

tr,w = β0,w +n

βi,w xi,r + κw,r

(5.1a)∀i ∈ {0, ..., n},(5.1b)

i=1

βi,w = γi + ui,w

where xi,r is the i:th element of input vector xr , κw,r is the within grouperror term, and γi is the slope effect of the dependent variable i ∈ {1, ..., n}.
For group w, u0,w is the intercept error and ui,w are slope errors ∀i ∈{1, ..., n}.
A two-level multilevel model is used to forecast each picker’s pick time
per batch. Following Bliese (2002), three important sources of variation
are found: within group (σ 2 ), between-group variation in intercepts (τ0 ),and between-group variation in slopes (τ1 ). The models are built using R
software and its packages, multilevel and nlme.
The data is divided in two parts. For each picker, a random 80% of thebatches are used for regression, and the rest (20%) for cross-validation
purposes only. By making this division, values for the cross-validation
residuals can be found for all batches used in the validation. These residuals give us an idea of the accuracy of the model along with other measures,such as the Akaike Information Criterion (AIC) and applicable versions of
R2 . As the forecasting models are directly used to calculate estimates ofbatch execution times in Section 6.4, the forecasts should be as accurate
as possible. In regressing, the batches are ordered by timestamp to accommodate possible learning-based autocorrelations in the data.
Plots of batch execution times and independent variables show thatvariance increases with total batch execution time, i.e., there is heteroscedasticity in the data, which appears to be largely remedied by applying a logarithmic transform of the output data (batch execution times).
First a level-1 model is constructed, i.e., a model with between-pickers
variability only in the intercept. Second, the level-1 model is extended to68Forecasting Batch Execution Time for Individual Pickers

level-2 by allowing between-picker variability in the slopes.
Step 1: intercept variation.

In this section, a test is made for signiﬁcanceof between-picker variation in the intercept. If the variation is not signiﬁcant, it will not matter which picker gets a job — everyone will perform ina similar manner. Hence, solving any assignment problem will be pointless. However, if there are between-pickers differences, one should strive
for a better assignment of work.
Independents are added to the model stepwise and tested whether there
is signiﬁcant variation between pickers. The variability of the intercept
term is examined with a level-1 Null Model (Bryk and Raudenbush, 1992):

tr,w = e(β0,w +κw,r ) E(eκw,r )(5.2a)

β0,w = γ0 + u0,w .

(5.2b)

where tr,w is batch execution time of batch r for picker (group) w, γ0 is thecommon intercept, κw,r the within-group error term and u0,w the betweengroup error term, and E(eκw,r ) is the Smearing Estimate (Duan, 1983) ofpicker w (i.e., the expected value of the retransformed residuals) used tocorrect the model bias resulting from the non-linearity of retransforming
the logarithmic dependent tr,w . In combined form, the model is tr,w =e(γ0 +u0,w +κw,r ) E(eκw,r ).
The Null Model has two possible sources of variance, τ0 for how muchthe groups’ intercept varies from the overall intercept (γ0 ), and σ 2 for thewithin-group variance. The Intraclass Correlation Coefﬁcient ICC(1) =τ0 /(τ0 + σ 2 ) (Bryk and Raudenbush, 1992, Kreft and De Leeuw, 1998)
equals 0.103, implying 10.3% of the total variance in the natural logarithm of time is due or related to differences between pickers. This suggests that it is beneﬁcial to assign the right batches to the right pickers.
Analysis of variance (ANOVA) is used to test for the signiﬁcant difference in -2 log-likelihood ratios between a model with and a model without
a random error term in the intercept (Bliese, 2002). A model with a high-2 log-likelihood is better than one with a low one, and the statistical signiﬁcance of the difference is tested. To achieve this, a comparison is madebetween the Null Model and a generalized least squares (GLS) ﬁt of asimilar intercept-only model which does not contain the between-grouperror term. The null hypothesis that no signiﬁcant difference in the -2
log-likelihood scores exists is rejected (p < 0.0001), so individual aspects

69Forecasting Batch Execution Time for Individual Pickers

explain up to 10.3% of the total variance in log-transformed batch execution time.
A Null Model can be built for non-transformed time, similarly as above.
In this case ICC(1) = 0.131, i.e., 13.1% of the total variance in time isdue to differences among pickers. The null hypothesis is rejected as well,so the result is signiﬁcant. Such differences among pickers can arise
from different physical or mental characteristics of the pickers, such as
strength or height, motivation, experience or stacking ability.
Step 2: slope variation and model selection.

The level 1 multilevel model

is now extended into a level-2 model by including group level errors for
slope. For n independent variables the model can be formulated as inequations (5.3) with the similar notation as in (5.1).

tr,w = e(β0,w +nβi,w = γi + ui,w

i=1 βi,w xi,r +κw,r )E(eκw,r )∀i ∈ {0, ..., n}(5.3a)(5.3b)A model is selected using a stepwise procedure. The selection procedure
is summarized in Table 5.3.
Each iteration, a new independent term is added to the current model.
ANOVA is used to test whether a random effect should be added to thenew model. This is done by comparing the -2 log-likelihoods of two versions of the new model: the ﬁrst with only a ﬁxed effect for the new term,
the second with both ﬁxed and random effects. Models are selected based
ﬁrst on the AIC and second on the sum of squared-error residuals (SS) of
the set of batches used for cross-validation.
The procedure begins by adding the ﬁrst term, number of lines (Lines)
to the Null Model (see Table 5.3). To better ﬁt the data to a normal
distribution, two models are compared: one with log-transformed input
data (Model 1L) and one with non-transformed data (Model 1X). Comparing AIC and SS for both models shows that Model 1L is better. A
random effect cannot be excluded (p < 0.0001). Next, the second independent variable (Travel) is added to Model 1L, resulting in models 2LX(non-transformed variable) and 2LL (transformed variable). Both models
are also tested for inclusion of random effects. The transformed input isbetter in terms AIC and SS for Model 2LL than for Model 2LX; therefore Model 2LL is chosen. The process continues by adding each of the70

Travel(2)
NA
NA
NA
***
***
***
***
***
***
***
***

Lines(1)
NA3
***
***
***
***
***
***
***
***
***
***

***

Mass(3)
NA
NA
NA
NA
NA
***
***
***
***
***
***

Level(4)
NA
NA
NA
NA
NA
NA
NA
***
***
***
***

Volume(5)
NA
NA
NA
NA
NA
NA
NA
NA
NA
***
17780

AIC4
46006
26631
16834
13922
13225
12675
12345
11996
11992
11688
2.5847e+09

SS 55.6719e+09
4.6986e+09
2.2053e+09
2.0428e+09
1.8963e+09
1.9143e+09
1.8727e+09
1.8517e+09
1.8510e+09
1.8290e+09
0.79

2 6
Rm
0
0.59
0.68
0.71
0.64
0.65
0.61
0.55
0.66
0.66

-

Rc2 6
0.10
0.68
0.82
0.84
0.86
0.85
0.87
0.89
0.87
0.87

-

p(rand)7<0.0001<0.0001<0.0001
0.0498<0.0001<0.0001<0.0001<0.0001<0.0001<0.0001of both the OLS and GLS models are the same).

9 An OLS/GLS multivariate regression over the whole data, with some of the variables log-transformed as with the best model 5LLLLX (the coefﬁcients

8 The regression did not converge for this model.

7 ANOVA test result on whether a between-pickers error term (random effect) should not be included in the model.

random (between-picker) factors (Nakagawa and Schielzeth, 2013).

6 R2 is used to describe the proportion of variance explained only by ﬁxed factors, while R2 is used to describe the variance explained by both ﬁxed andmc

5 Sum of squared residuals resulting from the cross-validation of the models.

4 The Akaike information criterion of the model, a measure of information lost if the model is used to represent the data. Lower is better.

3 Independent not included in the model.

2 Conﬁdence in the ﬁxed-effect independent, "***" indicates p ≤ 0.001, "-" that the independent is not signiﬁcant.

model. If an independent is added, the notation is: not included (N); in a non-transformed form with a random effect (X); or log-transformed form witha random effect (L). The addition of an autocorrelation model is denoted by (A). The last addition is the rightmost character.

1 The model number. Zero indicates the intercept-only model. The characters following the model number indicate the modiﬁcation done to a previous

Model1
Int(0)
0
***2
1X
***
1L
***
2LX
***
2LL
***
3LLX
***
3LLL
***
4LLLX
***
4LLLL
***
5LLLLX***
5LLLLL8
OLS/GLS9 ***Table 5.3. Multilevel model selection procedure

Forecasting Batch Execution Time for Individual Pickers

71Forecasting Batch Execution Time for Individual Pickers

log-transformed independent variables and the non-transformed ones, always testing for the signiﬁcance of the random effect. The tested models
are shown after Model 1L in Table 5.3. This results in Model 5LLLLX,
which has logarithmic inputs for all other terms apart from Volume.
The addition of interaction terms between the independent variables toModel 5LLLLX did not signiﬁcantly improve any ﬁt or error measure, or
the regression did not converge. The ﬁnal model in the notation of Bryk
and Raudenbush (1992) isa(β w , xr ) =β0,w + β1,w ln(x1,r ) + β2,w ln(x2,r ) + β3,w ln(x3,r )+ β4,w ln(x4,r ) + β5,w x5,r + κw,r
tr,w =ea(β w ,xr )βi,w =γi + ui,w

(5.4a))(5.4b)∀i ∈ {0, ..., 5}(5.4c)

E(eκw,rTo test for multicollinearity, a single multiple linear regression model
was built with the variables of Model 5LLLLX and run over all data.
The maximum VIF value was 4.4 (for ln(Lines)) suggesting that multicollinearity is not a major issue.
Figure 5.1 shows the cross-validation residuals grouped per picker. The
pickers’ residuals appear to be zero mean, and for approximately half ofthe pickers, the models seem to relatively accurate. However, there is stillclear heteroscedasticity in approximately half of the pickers — batch time
residuals are mainly skewed to the right, as is the case with the originaldata. However, the consistent reduction of the sum of squares of crossvalidation residuals SS as the model selection progresses in Table 5.3 isencouraging, despite the apparent heteroscedasticity.
Multivariate OLS and GLS models are ﬁtted using all independent variables, some log-transformed as in Model 5LLLLX, over the whole data set.
Between-picker differences are not taken into account. The model coefﬁcients in the ﬁtted OLS and GLS models coincide. ANOVA is used to test
the null hypothesis of Model 5LLLLX not having a signiﬁcantly better -2
log-likelihood value compared to the GLS model. The null hypothesis isrejected with p < 0.0001.
The resulting elasticities and standard deviations of model 5LLLLX andthe corresponding random effects can be found in Table 5.4. The betas(elasticities) show that the number of lines in a batch is the most important factor in the model. The second row shows the between-picker

72

Forecasting Batch Execution Time for Individual Pickers−20 20

64
●
●
●
●
●
●
●
●
●
●
●
●

Real time [min]

120
80
40
0

120
80
40
0

55

120
80
40
0

66

●●
●●●●
●
●
●
●
●
●
●

●●
●
●
●
●
●●
●
●
●
●
●●●
●
●
●
●
●
●

56

57

47

48

●
●
●
●
●
●
●●●
●
●● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

●●
●●
●
●
●
●●●●
●●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●

●●
●●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●

37

●
●
●●●
●●
●
●
●
●
●
●

38

●●
●● ●
●
●●
●
●
●
●●
●●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●

28

39

●● ●
●
●
● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

●
●
●
●
●
●
●
●
●
●
●
●

29

30

●

●

● ●
●●
● ●●
●
●●●
●
●
●●●
●●
●
●
●
●
●
●
●
●●

●
●●
●
●
●●
●
●
●
●
●
●
●

19

20
●

21

●

●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

●●
●
●
●
●
●
●●
●

10

11

●
●●
● ●
● ●
●●
●
●●
●
●
●
●
●
●
●
●
●
●

●
●●
●
●
●
●●● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

1

●●
●
●●
●
●
●
●
●
●●
●
●
●
●
●
●
●

−20 20

2

12

●

70

71

●●● ●
●
●
●
●●
●
●●
●
●●
●
●
●
●
●

●
●
●
●●●
●●
●
●
●
●
●
●

●●
●
●●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●

●●
●●
●
●
●
●●
●
●
●
●
●
●●
●
●

●
●
●
●
●●
●
●

120
80
40
0

●
●●
●
●
●
●
●
●
●

●
●
●
●
●
●
●
●
●
●
●

●●●●
●
●
●
●
●
●
●
●
●
●
●

●●
●
●
●
●
●

●●●
●
●
●
●
●
●
●
●
●
●
●
●●●
●
●
●
●
●
●
●
●●
●
●
●

●●
●
●●
●
●
●
●
●
●
●
●

●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

●● ●
●
●
●
●●
●●
●
●
●
●
●

● ●●
●
●
●
●
●
●
●
●
●

●
●
●
●
● ●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●

●●
●●
●●●
●
●
●
●
●
●
●●
●●

●●
●
●
●
●
●
●
●●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●

●●
●●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

●
● ●●
●
● ●
●●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

●
●
●●
●●
●
●
●
●
●
●
●
●●
●
●
●
●●
●
●
●
●
●
●
●●

●
●
●
●●
●
●
●
●
●●
●

●
●●
●
●
●
●
●
●
●
●

●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

●
●●
●●
●
●
●
● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

●
●
●●● ●
●●
●
●
●
●
●
●

●
●● ●
●●
●
●
●
●
●
●●●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●

●● ●
●●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

●
●●
●●● ●
●
●●
●
●
●
●
●
●
●
●●
●
●
●●
●
●
●
●
●●
●

●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

●
●
●
●● ●
●
●
●
●●
●
●
●
●
●●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●

●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●

●●
●
●●
●
●●
●●●●
● ●●
●●
●
●
●
●●
●
●
●
●
●
●
●
●
●

●
●
● ●●●
●
●
●● ●
●●
●
●
●
●
●
●
●
●

●●
●
●●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●

● ●●
●●●●
●● ●
●
●●
●
●
●●●
●
●
●
●
●
●●
●
●
●

59

49

60

50

40

31

●

32

●●
●
●
●
●
●
●●

22

● ●
●●
●●
●● ●
●
●
●
●
●
●
●
●
●
●
●●●

4

●●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●

−20 20

●
●

52

42

43

33

23

●
●●
●
●
●
●
●
●
●
●●

●

34

●
●●
●
● ●
●●
●
●●
●●
●
●●
● ●
●
●
●●●
●
●
●
● ●
●
●
●
●
●●
●
●
●
●
●
●
●
●
● ●
●
●
●
●●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

●

13

61

51

41

3

●
●●
●
●
●●
●
●
●
●●

69

●
●● ●
● ●
●
●
●●
●●●
●
●
● ●●
●●●●
●
● ●
●●●
● ●
●●
●
●
●
●
●
●
●
●●●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●

●●
●
●
●●
●
●
●●●
●
●
●
●●
●
●
●
●●
●
●
●
●
●●
●
●
●
●

−20 20

68

58

●●
●
●●
●●
●
●
●
●●
●

46

−20 20

67

●

●

●
●
●
●● ●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

●
●●
●●
●●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●

120
80
40
0

−20 20

65

14

●●
● ●
●
● ●
●
●●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●

5

●
●
●●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●

●

24

●●●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●
●●●
●
●
●
●
●●
●
●
●
●

15

●
●
●●
●
●
●
●
●●
●
●
●
●
●
●
●
●
●
●
●

6

●●

● ●
●● ●
●
●
●
●
●
● ●
●
●●
●
●
●
●
●
●

−20 20

25

●

16

7

●

−20 20

62

63

53

54

44

120
80
40
0

●●

45

35

●● ●

36

26

120
80
40
0

27

17

18

8

120
80
40
0

9
●

●

−20 20Residual [min]Figure 5.1. Cross-validation residuals for ﬁnal full model vs. batch time, grouped per
picker.

standard deviation of the model. The standard deviation of ln(Lines),ln(T ravel) and ln(M ass) are quite similar. The effect of mean picking
level seems to vary the most among the picker models.
The effect of increasing any of the log-transformed inputs is proportional
to the (re-transformed) output in this model. The dependent variable is
log-transformed, as are most of the independent variables, apart from Volume. In the case, when both are log-transformed, the scaled model beta
can directly be used to calculate the effect of an independent variable on
the dependent variable. For example, for the ﬁxed effect model, an 1%
increase in the number of lines in a batch (Lines) results in an 0.66%
increase in total batch execution time. For the non-transformed independent variable, Volume, a unit increase (1m3 ) will reduce the batch execution time by 0.16%. For Model 5LLLLX, the standard deviation of theresidual is 0.32, which means that if the forecasted time for a picker topick a batch is tw,r , there is an approximate 68% chance that it is accurate by a factor of e0.32 = 1.37, i.e, it is in the range [tw,r /1.37

1.37tw,r ]73

Forecasting Batch Execution Time for Individual Pickers

Table 5.4. Model betas for 5LLLLX with scaled data, effects on ln(T ime).
βRandom effect stdev

Intercept
3.2
0.40

ln(Lines)
0.66
0.09

ln(Travel)
0.16
0.04

ln(Mass)
0.17
0.08

ln(Level)
-0.44
0.38

Vol
-0.16
0.12

(Gelman and Hill, 2007, p. 62). Variance of batch execution time thusincreases proportionally to the value of time.
Table 5.4 shows that the batch execution time of the average picker
grows with the number of order lines, the distance traveled and total
batch mass, all of which seem reasonable: as the amount of work, travel
or mass increases, it is natural that time increases as well. A higher picklevel seems to make picking quicker, which suggests that currently the
pickers might be mostly picking from a non-optimal height. The negative
coefﬁcient in volume is somewhat more difﬁcult to explain. This might
result from the effect that some of the batches that contain relatively few
items of high volume are fast to stack. In some sense, the volume coefﬁcient can also reﬂect the stacking skill of the picker.
Picker clusters

Each of the pickers has their own forecasting model.

Generalizing a bit, and to explain the coefﬁcients in terms of individual
skill, one could characterize the pickers based on on the model coefﬁcientsin the following way:
• a low ln(Lines) represents a generally good worker;• a picker with a low ln(T ravel) is most likely a quick driver;• a low ln(M ass) represents strenght;• a low ln(Level) most likely means that the picker is tall, or has an easiertime picking items from the high level;• a low V ol represents a picker with good stacking ability.
To group and illustrate similarities among pickers, hierarchical clustering
is performed on the forecasting models using R software and the function
hclust. Euclidian distance between the picker coefﬁcients is used as themetric, and chosen mainly for its simplicity. The hierarchical clustering
routine is agglomerative, i.e., it ﬁrst assigns a cluster for each of the pickers and then iteratively proceeds to merge the most similar clusters. In

74

Forecasting Batch Execution Time for Individual Pickers

Figure 5.2. Dendogram of picker groups formed with hierachical clustering

Figure 5.2 the results are illustrated with a dendogram. Picker clusters
are numbered by size: 1 and 2 indicate the two largest clusters. If a cluster is divided, a minor number is added to the current clusters id, e.g.,cluster 1 is split into two subclusters indicated by 1.1 and 1.2 in Figure
5.2.

75

Forecasting Batch Execution Time for Individual Pickers

Table 5.5. Picker clusters, cluster sizes and mean coefﬁcients ordered by cluster heightin Figure 5.2.
Cluster id
1
2
2.1
2.2
2.2.1
2.2.2
1.1
1.2
2.2.2.1
2.2.2.2
1.2.1
1.2.2
2.2.1.1
2.2.1.2
Overall

size
26
73
1
72
35
37
10
16
7
30
7
9
12
23
99

ln(Lines)
0.69
0.65
0.70
0.65
0.65
0.65
0.72
0.67
0.75
0.63
0.60
0.73
0.71
0.62
0.66

ln(Travel)
0.17
0.16
0.15
0.16
0.16
0.16
0.16
0.18
0.17
0.16
0.17
0.18
0.14
0.17
0.16

ln(Mass)
0.18
0.17
0.15
0.17
0.17
0.18
0.17
0.19
0.09
0.20
0.22
0.17
0.17
0.17
0.17

ln(Level)
-0.04
-0.59
-1.42
-0.57
-0.71
-0.45
-0.20
0.05
-0.52
-0.43
-0.05
-0.13
-0.66
-0.74
-0.44

Vol
-0.15
-0.16
-0.48
-0.15
-0.17
-0.14
-0.22
-0.10
-0.11
-0.14
-0.05
-0.14
-0.26
-0.13
-0.16

Table 5.5 shows the mean picker coefﬁcients for each numbered cluster
in Figure 5.2 and the number of pickers present in each cluster. The pickers in cluster 2 are in general faster than those in cluster 1. In particular,
this can be seen from the differences in by the coefﬁcients for ln(Lines)
and ln(Level). It could be said, that the pickers in cluster 2 are more
adept in picking batches with larger numbers of order lines and beneﬁt
from more from the items being on a higher level. Continuing to follow
the dendogram from top to bottom, it can be observed that the next cluster split occurs to cluster 2: a single picker deviates from this cluster toform 2.1. He/she seems adept in picking items from a high elevation,and batches with relatively high volume. Cluster 2.2 contains the restof cluster 2, and is further split to 2.2.1 and 2.2.2, where the pickers in
2.2.1 seem to beneﬁt more from picking from a higher elevation whencompared to 2.2.2. Next, cluster 1, which contains (at least some of) theslow pickers, splits to 1.1 and 1.2: the pickers in 1.1 are slower in picking
batches with many order lines but the pickers in 1.2 are slower to travel
and adversely affected by pick level, where they have the only positive
coefﬁcient in ln(Level) among all the clusters. Cluster 1.2 also beneﬁts
less from batches with high volume. Clusters 2.2.2.1 and 2.2.2.2 are separated next. The pickers in 2.2.2.1 are less affected by batches with highmass than any of the other groups, but are not particularly fast in picking
batches with many order lines, while opposite is true for 2.2.2.2. Cluster
1.2.1 contains pickers that are good with batches with many order lines
but are affected heavily by batch mass and beneﬁt less from pick leveland batch volume. Cluster 1.2.2 contains pickers who are slow in picking
most batches with many lines, and don’t beneﬁt much from the pick levelor batch volume. Finally, clusters 2.2.1.1 and 2.2.1.2 contain the pickers

76Forecasting Batch Execution Time for Individual Pickers

that beneﬁt most from the high pick level but are either slow with batcheswith many orders and fast to travel (2.2.1.1) or are fast with many orderbatches (2.2.1.2).
From the analysis above, it can be seen that differences among thepicker clusters exist in all factors. The pickers have different strengths
and weaknesses that become apparent when comparing to the mean
picker (coefﬁcients).

5.4

Discussion

In this chapter, it has been shown that a signiﬁcant part, i.e., 13.1%, of
the total batch execution time is due to differences in pickers. This implies that it can be possible to ﬁnd the right picker for each batch. The
coefﬁcients of variation of the input data are all less than one suggestingthat batches are in general quite similar. A higher variance in the batch
parameters in the data (see Table 5.2) would most likely lead to a worse
ﬁt. This would result in less accurate picker models, and consequently tomore inaccurate forecasts. However, analysis of variance tests should inthis case indicate that something is amiss. Inaccuracies in the forecastingmodels would carry over to the optimization presented in the next chapter. A multilevel model is ﬁt to the cleaned-up data, using log-transformed
data to counter heteroscedasticity. The model is systematically built using a stepwise procedure based on improvements in both cross-validation
criteria and AIC. As AIC decreases, the cross-validation criteria generally
does as well. This conﬁrms the sanity of the model selection procedure.
All possible inputs (or their transformed versions) are shown to be statistically signiﬁcant. Furthermore, when comparing the sum-of-squares crossvalidation residuals between the OLS/GLS and ﬁnal multilevel models, it
is apparent that the multilevel model is superior. This reinforces the hypothesis that differences among pickers are signiﬁcant. Generally, transformed variables provide a better ﬁt than non-transformed ones. Thisis perhaps partly due to increased error probability in task execution as
the complexity of the task increases. In practice, the forecasting models
should be periodically updated to account for possible learning of pickers.

77Forecasting Batch Execution Time for Individual Pickers

78

6. Assigning Orders and Batches to
Pickers of Varying Skill

Data on the picking process and individual order pickers, readily available
in many warehouses, can be used to signiﬁcantly improve the efﬁciency ofthe order-picking process. This is illustrated for the case of a Finnish retailer operating a single large order picking warehouse supplying all of its
stores. By optimizing travel distance of order pickers, some savings in total batch execution time can be gained. Instead of considering just travel
time, if factors such as setup, search, and picking time are also included,
process efﬁciency can be signiﬁcantly improved. The results of Chapter 5
show that signiﬁcant differences exist among pickers. In this light, and as
no previous work has included picker skills in optimizing order picking operations, the forecasting models built in Chapter 5 are used to determine
which orders should be combined in a pick tour (i.e., a batch) and formulate the problem as a joint order batching and generalized assignment(BatchGAP) model. In this model, batches are assigned to pickers withthe objective of minimizing the sum of total batch execution times. Comparisons to the current state-of-the-art in order picking show that over
9% of total batch execution time can be saved by considering picker skill
and work assignment. These results can be used by warehouse managers
to either dispatch the right order to the best from the available pickers,
or in the case of hiring ﬂex pickers, to consult past log data and hire thosepickers matching the current orders.

6.1

Joint Batching and Generalized Assignment Problem

As the objective of the optimization presented in this section, customer
orders and batches are assigned to those order pickers who have the best
skills to execute them, thus minimizing the total batch execution time(including pick, travel, and setup time).

79

Assigning Orders and Batches to Pickers of Varying Skill

Order batching for batch sizes of three or more orders has been shown
to be NP-hard (Gademann and Van de Velde, 2005). To further assign
batches to pickers, the order batching model from Gademann and Van de
Velde (2005) is extended to include a generalized assignment problem,
which is also NP-hard (Fisher et al., 1986). This results in a joint batching
and generalized assignment problem (BatchGAP), which is sufﬁciently
complex to justify the use of a heuristic.
The problem is deﬁned as follows. Let R be the set of all possible batch
combinations from the set of orders O, each batch consisting of a maximum of N orders. The orders that are to be picked during a day arrive
early each morning. Let Rs ⊆ R be the set of batches present in solutions ∈ S, where S is the set of feasible and complete solutions to the BatchGAP. Furthermore, let W be the set of modeled pickers. The mapping ofeach customer order in O to a batch r ∈ R is characterized with a zero-one
vector ar of length |O|. If aor = 1, order o is included in batch r, otherwise
it is not. Each batch r ∈ R has multiple order lines that need to be picked
in no particular order during a picking tour y ∈ Yr , where Yr is the set of
all possible tours for batch r.
After picking all lines of an order, the order must be left at an orderspeciﬁc drop-off location in the warehouse. All pick and drop-off locations
must be visited at least once. Thus a tour is a solution to a TSP, where thedrop-off locations have to be visited after the last pick in the batch (seeFigure 3.1). All picking tours start and end at a depot. Each batch r ∈ R
has a travel distance of dr associated with the tour-construction method
used. The forecast total batch execution time, tw,r , depends on the parameters of batch r and the forecasting model of picker w. When calculating
the forecast tw,r , a TSP corresponding to the batch r to get the travel distance dr needs to be solved. The TSP is solved separately with a heuristic(see Section 6.3.1). The TSP constraints are not included in the model below, as the model is mainly illustrative of the problem, and the additional
constraints would unnecessarily complicate it. A heuristic is used to solve
the TSP as no fast optimal solution to a precedence-constrained TSP (notethat the last location of each order is ﬁxed) exists for a warehouse with a
middle cross-aisle and unidirectional travel. The tour length of a TSP related to a particular batch is used as one of the inputs when forecasting
tw,r using equations (5.4). Each picker w ∈ W has a maximum working
time during a shift, Mw .
The following optimization model can be formulated for the BatchGAP:

80

Assigning Orders and Batches to Pickers of Varying Skill

min

 

tw,r Xw,r(6.1)

Xw,r ≤ 1 ∀r ∈ R(6.2)

tw,r Xw,r ≤ Mw

(6.3)

w∈W r∈R

subject to


w∈W



∀w ∈ Wr∈R

ao,r ≤ N∀r ∈ R(6.4)

ao,r Xw,r = 1 ∀o ∈ O

(6.5)

o∈O

 

w∈W r∈R

Xw,r ∈ {0, 1}(6.6)

The goal is to ﬁnd a solution s∗ such that the set of batches Rs∗ ⊆ R form
a complete partitioning of O and have a feasible assignment to the pickers
W, such that the sum of costs of all pickers is minimized and constraints(6.2, 6.3, 6.4 and 6.5) are upheld.
The objective function (6.1) minimizes the sum of all the pickers’ batch
execution times. The binary decision variable Xw,r is one if batch r is
assigned to picker w, otherwise it is zero. Constraints (6.2) enforce thateach batch is picked at most once. Each pickers’ maximum working time,
Mw is enforced by (6.3). Constraints (6.4) enforce that no more than Norders are contained in any one batch. Each order must be allocated once
to any chosen batch (6.5) and integrality of the decision variable Xw,r isenforced by (6.6).
As a by-product of the optimization, the effect on picker productivity isalso of interest. Let lr be the number of order lines in batch r. Productivityof a picker w ∈ W is deﬁned as
Lw =


r∈R

Xw,r lr /Xw,r tw,r

(6.7)

r∈R

or the average number of order lines picked per unit of time worked (usually minutes). Notice that the total work time differs from maximum shift
length. This study aims to model the maximum the maximum total timethat can be saved — salary models are not considered.

81

Assigning Orders and Batches to Pickers of Varying Skill

6.2

Note on Problem Complexity

The complexity of assigning orders to batches and batches to pickers canbe illustrated by the number of combinations of a related bin-packing
problem, where batches and pickers are both homogeneous. The number of possible batch combinations of up to the maximum batch size of N
can be found from (6.9). Let |P| be the total number of batches, and let L
be the number of batches a picker can complete during a shift. Assuming|P| = |W|L, the total number of possible bin-packing combinations for |W|
pickers is
Cassign =

|W|


⎛
⎝

|P| − iL

⎞
⎠.

(6.8)

L

i=1

A regular shift at the sample warehouse deploys by 30 pickers. Let one
picker take half an hour to process one batch, and work 8h each shift, fora maximum of 16 batches/shift. If all pickers and batches are to ﬁlled tocapacity, there are 480 batches and 1440 orders (the maximum batch sizeis thus 3), and the number of joint bin-packing and order combinations is⎛
Ctotal = ⎝

1440
3

⎞
⎠

29

i=0

⎛
⎝

480 − i16
16

⎞
⎠,

(6.9)

which is a very big number, approximately 10152 . To ground this number
somewhat, the total number of atoms in the universe is estimated to be1082 . Optimally solving the joint problem is impractical even for a small
number of pickers.

6.3

Solving the Joint Batching and Generalized Assignment
Problem

6.3.1 Routing Heuristic
To calculate total batch time, total travel distance of a picking tour isneeded. The travel distance is assumed to be independent of the picker.
The sample warehouse has multiple drop-off locations (more than one ofwhich can be visited during a tour) and a middle cross-aisle (see Figure3.1). The drop-off locations can only be visited after all items of an order
have been collected. Each aisle can only be traveled in a single direction,
while the cross-aisles are bidirectional. Tours start and end at the depot.

82

Assigning Orders and Batches to Pickers of Varying Skill

The number of potential batch evaluations for the problem introducedin Section 6.1 is potentially very large. For each batch evaluation, therouting problem needs to be solved. To calculate total travel distance, acomputationally light heuristic is introduced.
In Algorithm 6, a version of the aisle-by-aisle routing heuristic(Vaughan, 1999) for warehouses with unidirectional travel in the aisles
is presented. In this version of the algorithm, aisles do not need to be
traveled completely if there are no picks past the middle aisle. It also
incorporates drop-offs before the eventual return to the depot. The basic
idea is to take advantage of the unidirectional travel in the warehouseand the middle cross-aisle. Once a path to pick all order lines has been
formed (steps 1-4 and 6-8), the tour is optimally completed by adding thedrop-off locations with Dijkstra’s algorithm.
Input: A batch r, order drop-off locations Lr , depot D, distance matrix
B, set of aisles A
Output: A tour y, travel distance dr
1 Sort all orders in r according to the aisle index.
2 Make a vector v a of corresponding order lines for each relevant aislea ∈ A to batch r.
3 Sort each aisle vector v a according to the travel direction of the aisle.
4 Insert all vectors from the smallest aisle index to the largest index to apath ysmall .
5 Use Dijkstra’s algorithm to make ysmall into a tour by adding Lr and D
to it.
6 Store the total path of ysmall cost in dsmall using B.
7 As in 4 but iterate from big to small aisle index and store the path inybig .
8 Store the total path of ybig cost in dbig using B.
9 Use Dijkstra’s algorithm to make ybig into a tour by adding Rr and D
to it.
10 if dbig < dsmall then11y ← ybigdr ← dbig
12
13 else14y ← ysmall
15
dr ← dsmall
16 endAlgorithm 6: Middle aisle multi-dropoff routing heuristic

6.3.2

Initial Assignment of Batches

In many warehouses, the current practice is to assign the next batch in
the queue to a picker who is ﬁrst available to execute it. No planning orscheduling is involved. This method of assigning jobs to pickers is hence-83Assigning Orders and Batches to Pickers of Varying Skill

forth called the ﬁrst-free assignment method. It assumed that ﬁrst-free isused in the example warehouse.
Before assigning batches, an initial feasible solution must be found.
This is done with a combination of a batching and an assignment algorithm. As the batching algorithm, either C&W(i) (Clarke and Wright,
1964) or VNS (Albareda-Sambola et al., 2009) is used. For assigning
batches to pickers, ﬁrst-free is applied.

6.3.3 Assignment of Customer Orders Based on Average Picker
Productivity
Based on the data, the average picker productivity can be calculated using(6.7). This statistic is tracked in many warehouses already. A warehousemanager can choose to prioritize assigning work to the most productive
pickers. To reﬂect this, and to try to justify the beneﬁts of skill based
assignment, the following heuristic, fastest-ﬁrst, is used to assign work
to the most productive pickers ﬁrst. It will usually result in the slowest
pickers being left out of the assignment. This algorithm assumes that
there are as many pickers available as needed.
Step 1. Batch all orders based on travel distance using VNS.
Step 2. Sort all available picker based on productivity in descendingorder.
Step 3. Sort all batches based on the number of order lines indescending order.
Step 4. Assign the ﬁrst unassigned batch to the picker who has thehighest productivity value and so that the picker’s maximum
worktime is not exceeded. Repeat until all batches are assigned,otherwise exit.

6.3.4 Adaptive Large Neighborhood Search Algorithm
After forming a feasible initial solution, Adaptive Large Neighborhood
Search (ALNS) (Ropke and Pisinger, 2006a) is used to search for good
local optima. Pisinger and Ropke (2010) offer guidelines for designing
ALNS algorithms, which have been mostly followed in this section. The
ALNS algorithm uses a set of neighborhood heuristics H. For the ALNS
presented in this section, H is composed of ﬁve different heuristics. See
Section 2.2.6 for more details of ALNS heuristics.

84Assigning Orders and Batches to Pickers of Varying Skill

Following Ropke and Pisinger (2006a), a simulated annealing-type hill
climbing scheme is implemented to complement the search. This simulated annealing scheme uses a geometric cooling schedule (Cohn and
Fielding, 1999). At each iteration, the temperature is updated by T ← φT ,where T is the temperature variable and φ the cooling coefﬁcient.
A vector collecting all the adjustable parameters of the ALNS algorithm
is deﬁned as v = (σ1 , σ2 , σ3 , λ, δ, φ). The actual parameter set is given in
Section 6.4.1.
Next, the neighborhood search heuristics are detailed. All of them maintain the feasibility of the solution by enforcing tw,Rw ≤ Mw , where tw,Rw
is the total time to pick the set of batches Rw ⊆ R assigned to picker
w and Mw is the maximum working time for picker w. Only those solutions that are feasible according to the model in Section 6.1 are accepted.
These heuristics are used to minimize total batch execution time, and an
allocation can result in some of the pickers being left without work.
1 Random destroy, random repair

Step 1. Randomly choose Q ∼ U (2, 7) batches to be destroyed from anyof the pickers, resulting in the set of Od customer orders, where|Od | ≥ Q.
Step 2. Form a set Wd of pickers currently assigned to Od .
Step 3. Randomly choose a number of orders up to the maximum batch
size of N orders from Od . Form a new batch from the chosen orders
and assign it to a random worker w ∈ Wd if constraints (6.3) can be
upheld. Remove the chosen orders from Od ..
Step 4. Repeat Step 3 until there are no orders left to assign.
Step 5. Calculate the tour costs using Algorithm 6 for the new batches
and form completion time forecasts for each new picker batch pair.
2 Savings-based destroy, random repair

The aim of this heuristic is todestroy those batches that have the least savings value (the time difference of picking each order separately minus picking them together in abatch). Let s be a current solution to the BatchGAP and Stot ∈ R be thesum of time-savings of all batches in s. Furthermore, let Smax ∈ R be themaximum batch savings value in s, and Sr be the savings of batch r ∈ Rs .
Then the probability of destroying a batch r is
pr = (Smax − Sr )/Stot .

(6.10)

85

Assigning Orders and Batches to Pickers of Varying Skill

Step 1. Calculate the batch savings for all batches in s.
Step 2. Use the roulette wheel method from (2.16) with batch
probabilities calculated from Equation (6.10) to choose Q ∼ U (2, 7)
batches to destroy from any of the pickers, resulting in M ≥ Q totalcustomer orders.
Step 3. Form a set Wd of pickers currently assigned to the batches to be
destroyed.
Step 4. Randomly choose a maximum of N orders from the destroyedbatches, form a new batch and assign it to a random worker w ∈ Wd
if constraint (6.3) can be upheld.
Step 5. Repeat Step 4 until there are no orders left to assign.
Step 6. Get the tour costs for the new batches and form completion time
estimates for each new picker batch pair.

3 Move batches

This heuristic tries to move a set of batches T ⊆ R be-

tween pickers while leaving the order composition of the batches intact. A
batch r ∈ T from picker wi ∈ W is moved to another picker wj ∈ W − {wi }.
No routing of batches is necessary, the work time estimates of the involved
pickers just need to be calculated.
4 Empty picker and greedy repair

With this heuristic, all batches are

moved from a single picker to other workers. First, the picker is selected
randomly from all pickers with batches. Then, each of the selected picker’s
batches is iteratively assigned to the picker who based on the forecast will
execute it the fastest, and so that the capacity constraints for each picker
are upheld.
5 Move all batches from a picker to another

This ﬁnal heuristic moves

all batches from one random picker to the picker that executes them thefastest.

Acceptance of new solutions
After applying a neighborhood search heuristic to a current solution s, a
new solution s is found. Let s∗ be the best solution found so far. Furthermore, let f (s) be the cost of solution s. If f (s ) < f (s∗ ), the solution is
accepted setting s∗ = s and s = s The current score of the active heuristic is updated by Ψh = Ψh + σ1 . Otherwise, if f (s ) < f (s), the solution is86Assigning Orders and Batches to Pickers of Varying Skill

accepted as the current best one setting s = s and Ψh = Ψh + σ2 . A worse
solution than the current one can be accepted by the simulated annealing
rule. Let Q ∼ U (0, 1) and T be the current temperature of the algorithm.
(s)), a hill climb is performed and setting s = s and
If Q ≤ exp( f (s )−fTΨ h = Ψh + σ3 .

6.4

Results

In this section, it is shown how original batch times relate to the modeled
ones and compute new batching, routing and assignment solutions which
are compared to the original ones. All algorithms were coded in C++ andrun single-threaded on an Intel Xeon 3.2 GHz with 16GB of memory.

6.4.1

Parameter Calibration

Following Ropke and Pisinger (2006a), the calibration parameter vector
is initialized to v = (σ1 , σ2 , σ3 , λ, δ, φ) = (33, 13, 9, 0.95, 200, 0.999995) withthe exception of δ, which was not used by these authors; its value was set
according to Gharehgozli et al. (2013). In order to calibrate the parameters, two sample instances were run of 3 pickers with 78 orders, and 10
pickers with 360 orders. One parameter was varied at a time on a uniform
interval for each parameter around the initial value. The best performing parameters were chosen and set as v = (65, 13, 9, 0.95, 91, 0.999995).
The scaling factor of the computational complexity of heuristic h, ζh from
(2.18), is not used, as no signiﬁcant improvement in solution quality or
time was found when it was used. The starting temperature T0 is set
such that the probability of an uphill climb is 0.5 if the new solution is 3%
worse than the sum of batch execution times in the starting solution s0 by
T0 = −0.03f (s0 )/ ln 0.5.

6.4.2

Data Preparation

During a daily 8h shift in the sample warehouse, the number of available pickers can vary between 20 and 40. Since the data do not contain
due times for the batches, and only a part of all the pickers appearing in
the data can be used due to data cleaning, comparing the allocations andbatches of the real shifts and the ones reconstructed by the ALNS is notnecessarily meaningful. To still allow a comparison, virtual days are constructed by partitioning the data of 24,669 batches into 12 (about) equally87Assigning Orders and Batches to Pickers of Varying Skill

sized subsets. Each partition forms one virtual day, resulting in approximately 2056 orders (685 batches) per day. A picker qualiﬁes for inclusionin the workforce of the virtual day if both the sum of real execution timesof the batches he or she performed and the sum of the forecast batch execution times exceed the minimum threshold of Mmin . This is done toguarantee that each picker included carries out sufﬁcient work. The total
execution time of all the batches executed by a picker should not exceed
the maximum threshold time Mmax . If adding a batch were to exceed thisthreshold, it would be discarded for this picker.
The used data only contain Hamiltonian Paths starting from the ﬁrst
pick line to the last line, leaving out the roll container drop-offs and eventual return of the picker to the depot. Each order has its own uniquedrop-off location. Figure 3.1 shows the different locations of these dropoffs. Dijkstra’s algorithm is applied to ﬁnd the combined minimum distance of visiting the order drop-off locations and further traveling fromthe last drop-off location to the depot and from the depot to the ﬁrst picklocation. Using this distance, the extra time to complete the pick tour
is now calculated, assuming an average truck speed of 1 m/s, and taken
into consideration when calculating sums of batch execution times in constructing a virtual day.
Now the real batch execution times and their forecasts for these newly
created virtual days are known. These can be used to compare the ALNS
batching, routing, and assignment results with the original solutions.

6.4.3 Experimental Setup
Virtual days are constructed as detailed in Section 6.4.2. For each virtual
day, the constraints in (6.3) are set to 8h. Two different instances are run
by setting Mmin to 7.15h or 7.4h and Mmax to 7.75h in both real time andforecast time. This results in 29 and 20 pickers working on the virtualdays on average, respectively. Out of the possible 99 pickers, 93 qualify
for the simulations with Mmin = 7.15h and 86 for Mmin = 7.4h. The ALNS
is run for 5,000,000 iterations in both cases.
Work is currently issued according to the ﬁrst-free allocation.

Twobatching algorithms, C&W(i) and the Variable Neighborhood Search(VNS) by Albareda-Sambola et al. (2009), are used to generate an initialsolution for the ALNS. C&W(i) provides a quick solution to the batching problem, while the VNS method is shown to provide good resultswith the tradeoff of exploring large parts of the solution space (Albareda-88Assigning Orders and Batches to Pickers of Varying Skill

Sambola et al., 2009, Matusiak et al., 2014). Using these initial solutions, the ALNS is used in two ways: (1) including all neighborhood searchheuristics, and (2) allowing only those heuristics which do not break downbatches (i.e., excluding heuristics 1 and 2).
In the ALNS heuristic presented in this chapter, a computationally lightrouting heuristic (see Section 6.3.1) is used, as it has to be invoked numerous times. In order to understand how far this routing is off a nearheuristic, it is compared with LKH (Helsgaun, 2000, Theys et al., 2010)
for 10,000 random batches constructed from customer orders from thedata-set. LKH has a gap of about 0.1% with optimal routing (Theys et al.,2010), at the expense of much longer computation times. LKH does nothandle precedence constraints, so drop-offs are left out of the tours for thiscomparison. That is, in this experiment the tours of both the light routing heuristic LKH lack precedence constraints. On average, LKH leads tosolutions of about 11.7% shorter tours.

6.4.4

Time Savings and Comparison to Original Solution

In this section, the algorithmic solutions are compared to the original andinitial batching allocation, as well as to each other. The following abbreviations are used for the solutions in Figures 6.1 and 6.2. Unless otherwise
noted, batch execution times are forecast:
BF batching done with VNS and allocation using fastest-ﬁrst;
CWI batching done with C&W(i) and allocation using ﬁrst-free;CWR improvement of CWI using the ALNS algorithm with no further
rebatching (only heuristics 3, 4, and 5 are used);
CWF as CWR, but with all neighborhood search heuristics used in the
ALNS;
VNSI batching done with VNS and allocation using ﬁrst-free;
VNSR improvement of VNSI using the ALNS algorithm with no further
rebatching (only heuristics 3, 4, and 5 are used);
VNSF as VNSR, but with all neighborhood search heuristics used in the
ALNS.
The comparison baselines in ﬁgures 6.1 and 6.2 are the forecast batch
execution times with real batches during a virtual day. Both ﬁgures are89Assigning Orders and Batches to Pickers of Varying Skill

13
12% time saved compared to forecasted

11
10
9
8
7
6
5
4
3
2
1
0

BF

CWI

CWR

CWF
VNSI
Solution type

VNSR

VNSF

Figure 6.1. Savings generated by different algorithms, compared to the original batchingwith forecast execution times during a virtual day. 7.15h minimum time toqualify, 29 pickers. The real batch execution times for the original batchingare 3% smaller than the forecast in this case.

13
12% time saved compared to forecasted

11
10
9
8
7
6
5
4
3
2
1
0

BF

CWI

CWR

CWF
VNSI
Solution type

VNSR

VNSF

Figure 6.2. Savings generated by different algorithms, compared to the original batchingwith forecast execution times during a virtual day. 7.4h minimum time toqualify, 20 pickers. The real batch execution times for the original batchingare 1% smaller than the forecast in this case.

90

Assigning Orders and Batches to Pickers of Varying Skill

quite similar, so the amount of spare time in the virtual day set-ups or
the number of pickers does not affect the performance of the algorithms
much. The forecast in Figure 6.1 gives a pessimistic estimate, since thereal times are about 3% smaller. In Figure 6.2 the real times are 1%
smaller. If jobs are assigned based on the average productivity (categoryBF), savings of around 6% can be achieved in both cases. The initial solution CWI gives improvements of over 2% compared to forecast batch
execution times with real batches. The best solution is provided by the
combination of VNS batching and the ALNS, i.e., solution VNSF. However, CWF is very close in terms of solution quality and much faster to
execute. In the 29-picker case, CWF took on average 80 minutes to solve
the problem of batching and assigning 2048 orders, whereas VNSF took
120 minutes. For the 20-picker case, these times are 35 minutes and 70
minutes, respectively. In the sample warehouse, orders arrive early inthe morning a couple of hours before the morning shift starts, so there is
ample time to run the algorithm. Both ﬁgures show that it is beneﬁcial torebatch using the neighborhood search heuristics 1 and 2 in addition to 3,
4, and 5, which just move batches between pickers.
Most importantly, Figures 6.1 and 6.2 show that signiﬁcant savings in
total batch execution time can be obtained by assigning the right batches
to the right pickers. Almost 10% more time can be saved when comparing VNSF to CWI and to VNSI. Compared to the forecast and real time
for the original allocation, between 10% and 12% can be saved. Finally,
when comparing VNSF with BF, i.e., assigning based on picker skill vs. on
productivity, savings of 6% in total batch execution time can be achieved.

6.4.5

Impact on Picker Productivity

Figures 6.3 and 6.4 show the effect of the VNSF solution compared to a
VNSI solution, which implies state-of-the art batching with a ﬁrst-freepicker assignment. Pickers are divided into three categories based on
their productivity (lines picked per time unit worked) in the VNSI solution: the slowest 20%, the medium 60%, and the fastest 20%, with a totalof 20, 53, and 20 pickers, and 15 pickers (19, 48, and 19 in Figure 6.4),
respectively. However, as noted below, the slowest pickers receive muchless work. Here productivity is calculated by the total number of assignedlines divided by the total time to process all assigned batches, i.e., actual
work time.
All categories improve productivity, while most improvements can be ob-91Assigning Orders and Batches to Pickers of Varying Skill

VNSF improvement in productivity over VNSI [%]

25

20

15

10

5

0

Slow (20)Medium (53)Fast (20)
Slowest 20%, medium 60%, fastest 20% picker

Figure 6.3. Change in average productivity by picker category. 7.15h minimum time toqualify, 29 pickers

served in the slow and medium picker categories. For pickers in the slow,
medium, and fast categories for the case of 29 pickers (Figure 6.3), theinitial average productivity values are 1.14, 1.43, and 1.81 lines/minute,respectively. The VNSF solution improved productivity by 19%, 7%, and
12% on average compared to the VNSI for the same categories. For the
20 picker case, VNSI productivities are 1.12, 1.44, and 1.83 lines/minute,which improve by 17%, 5%, and 15% (Figure 6.4). The improvement in all
categories is a result a reduction in total working time compared to theinitial allocation as picker skills are taken into consideration. Additionally, some pickers, particularly those in the slow category, are assigned
less work, and in many cases do not get to work full shifts.
Since VNSI is used as a reference, this improvement in productivity
must result from better picker assignment. However, the slowest pickers
receive fewer orders to be picked, resulting in 40% to 60% fewer pick lines,
as can be seen from ﬁgures 6.5 and 6.6. The medium and fast pickers tend
to get more work in both the 29- and the 20-picker cases.
To better understand which of the pickers are included in the ALNS
assignment over the entire horizon of eight virtual days, the picker categories used above are divided into subcategories, indicated by two characters, XY. X stands for the main categories S, M, F (Slow, Medium, Fast),

92

Assigning Orders and Batches to Pickers of Varying Skill

VNSF improvement in productivity over VNSI [%]

25

20

15

10

5

0

Slow (19)Medium (48)
Fast (19)
Slowest 20%, medium 60%, fastest 20% picker

Figure 6.4. Change in average productivity by picker category. 7.4h minimum time toqualify, 20 pickers

20% change in no. of pick lines assigned

10
0
−10
−20
−30
−40
−50
−60
−70

Slow (20)Medium (53)Fast (20)
Slowest 20%, medium 60%, fastest 20% picker

Figure 6.5. Change in number of lines assigned by picker category. 7.15h minimum timeto qualify, 29 pickers.

93

Assigning Orders and Batches to Pickers of Varying Skill20% change in no. of pick lines assigned

10
0
−10
−20
−30
−40
−50
−60
−70

Slow (19)Medium (48)
Fast (19)
Slowest 20%, medium 60%, fastest 20% picker

Figure 6.6. Change in number of lines assigned by picker category. 7.4h minimum timeto qualify, 20 pickers.

based on the productivity in the VNSI solution. Y has three categories,
0, 1, 2, indicating whether the picker: is left out of the assignment completely, or receives less than a ﬁfth of his or her total number of order lines
assigned in the VNSI solution, or receives more than a ﬁfth of his or her
total number of order lines assigned in the VNSI solution, respectively.
Figures 6.7 and 6.8 show the division of pickers over these categories
compared to VNSI. The pickers in the slow category are the most affected:
in the 29-picker per day case, two out of 15 pickers are completely left
without work for all eight days. In the 20-picker case, one picker is leftout. In the medium category, ﬁve and four pickers left out of the assignment, whereas all pickers in the fast category are assigned to work. Apparently, it pays off not to use some of the workers at all and to use someof the workers as little as possible.

6.4.6 Effect of Skill on Batch Assignment
Forecasting the time for a picker to pick a batch is done by using the model
described in Section 5.3. By including the knowledge of the batch parameters in the optimization process of the joint batching and generalized
assignment problem, a better result can be achieved than optimizing by

94

Assigning Orders and Batches to Pickers of Varying Skill

50
45

Number of pickers in category

40
35
30
25
20
15
10
5
0

S0

S1

S2

M0

M1

M2

F0

F1

F2

Figure 6.7. Picker categories after ﬁnal assignment. 7.15h minimum time to qualify, 29
pickers on average.

45
40

Number of pickers in category

35
30
25
20
15
10
5
0

S0

S1

S2

M0

M1

M2

F0

F1

F2

Figure 6.8. Picker categories after ﬁnal assignment. 7.4h minimum time to qualify, 20
pickers on average.

95

Assigning Orders and Batches to Pickers of Varying Skill

travel distance and statistical productivity (picking speed). The productivity of a picker is a statistic already tracked in many picking processes.
However, to the author’s knowledge, it is not extensively used in assigningwork: but neither are other statistics. It is common to treat order pickers
as a homogeneous unit.
Setup

This following experiment tries to provide an insight to how muchtotal order picking time can be saved by taking into account the additional information of batch parameters and picker models when compared
to good travel distance-based batching and assigning work based on the
pickers’ productivity. If slow pickers are left out, it will improve the overall performance of the process, as seen previously in sections 6.4.4 and6.4.5.
To start, 12 virtual days with a 7.4h initial threshold are constructed,
similarly as explained in Section 6.4.2. Based on this, a combination of
two algorithms is used to optimize the order batching and to assign work.
For the batching, the state-of-the-art VNS algorithm of Albareda-Sambola
et al. (2009) is used. Fastest-ﬁrst, which optimizes total batch time basedon the most commonly measured statistic of picker productivity, is applied
to assign batches to the most productive pickers. This setup will leave out
the least productive pickers. Finally, if any pickers exist that work forless than 7.4h after the fastest-ﬁrst assignment, they are left out of thisexperiment. This results in an allocation where only the statistically mostproductive pickers are selected, with each of the pickers being assigned aminimum of 7.4h of work. It is relatively tight in terms of free time sothat not many pickers will be left out when improving it with the ALNS.
After constructing the initial allocation with the VNS batching andfastest-ﬁrst assignment, the ALNS algorithm is run to improve upon thisinitial solution. Both the ALNS and fastest-ﬁrst try to minimize total
batch execution time over all batches and thus are comparable. The resulting allocation takes into account the additional information of batch
parameters and picker coefﬁcients that affect the output time.
Results

The average number of qualifying pickers per 2048-order virtual

day is 19. Total time savings after applying the ALNS are over 5.5%. As aresult of applying the ALNS, two pickers of the 80 that participate in thesevirtual days are always left out work: the ﬁrst on two separate days, andthe second on one day.
On average, productivity of the remaining pickers increases by 4.6%.
The average time to complete an order line is decreases from 39.8s to96Assigning Orders and Batches to Pickers of Varying Skill

37.8s, or 5.5%. As picker productivity is deﬁned as the number of orderlines picked per time unit in (6.7), its increase fails to explain the savings
completely. Thus it can be concluded that 0.9% of the time savings result
from leaving out the two pickers. The increase in picker productivity must
therefore come from a better assignment of orders (and batches) between
the participating pickers.
To illustrate how orders and batches are assigned to the most suitablepickers, the batch parameters (i.e., Lines, Travel, etc.) of each batch assigned to the pickers are logged for both the initial allocation constructedwith fastest-ﬁrst and the ﬁnal allocation after applying the ALNS. These
data are then used to calculate the average batch assigned to each picker.
A low forecasting model coefﬁcient corresponds to that particular picker
being skilled in that category (see Section 5.3). The hypothesis was madethat savings that are found by applying the ALNS-algorithm are at leastin part due to difference in picker skill: pickers with low coefﬁcients ina particular category will generally be assigned batches with high values in the same category. For example, if a picker is more skilled inpicking heavy items, i.e., he/she has a low model coefﬁcient in ln(M ass),and he/she should be generally assigned heavier batches. For each of themodel inputs (batch parameters), this hypothesis can be veriﬁed from ﬁgures 6.9-6.13. In these ﬁgures, each of the model coefﬁcients of the participating pickers is plotted against the average batch parameter of thebatches that are assigned to those pickers. Figures 6.9 show how ln(Lines)
relates to the number of lines of in a batch and ﬁgures 6.10 show the relation of ln(T ravel) to the average travel distance of a corresponding batch,and so on. In the top ﬁgures, marked (a), batch assignment is done withfastest-ﬁrst, and in the bottom ones, marked (b), ALNS is applied to the
BatchGAP. The line in each of the ﬁgures shows an OLS ﬁt with the picker
model coefﬁcients as the dependent variables and the batch parameter as
the independent variables. These ﬁgures indicate that the ALNS tends toassign batches with high values in a particular parameter to those pickers with a corresponding low model coefﬁcient, while the combination of
VNS order batching and fastest-ﬁrst does not take into account skill. Mostof the savings gained from applying the ALNS result from assigning theright orders and batches to the right pickers, as they cannot come solely
from leaving out the slowest pickers.

97

Assigning Orders and Batches to Pickers of Varying Skill

●

0.8

●

●

●
●

●
●

0.7

●

●

●

●
●

●

●
●

●

●
●

●●
●

●
●●●
●

●

●

●

●

●
●

●

●

●

●

●
●

● ●
●
●●
●

●

●

●

●
●

●

●

●

●

●

●

●

●

●

●

●

●

●
●
●●

●

●
●

●

●

●

0.6

log(lines) picker model coefficient

●

●

●

●
●

●

●

0.5

●

●

0

20

40

60

80

100

120

140

Number of order lines, fastest−first allocation of batches(a) VNS batching and fastest-ﬁrst assignment

●

●

●

●
●

0.7

●
●
●

●

●

●

●

●

●

●

●
●●
● ●
●

●
●
●

●

●

●

●

●

●●

●● ●
●

●

●
●

●

●

●
●
●
●●

●
● ●

●
●

●

●
●
●

●

●
●
●●

●

●
●
●
●
●
● ●

●

●

0.5

●

●
●

●
●

0.6

log(lines) picker model coefficient

0.8

●
●
●

●

0

20

40

60

80

100

120

140

Number of order lines, ALNS(b) ALNS
Figure 6.9. The average number of order lines in a batch assigned to a picker vs. thecorresponding ln(Lines) forecasting model coefﬁcient of that picker, for two
different ways of solving the BatchGAP.

98

Assigning Orders and Batches to Pickers of Varying Skill

●
●

●

●

●

●

0.20

●
●

●

●

log(travel) picker model coefficient

●

●
●

0.15
●

●

●

● ●

●
●
●
●●

●

●

●

●

●

●
●
●

●

●
●
●

●

●

●

●

●

●

●

●

●

●

●

●

●
●●

●

●

●

●

●

●

●

●

●

●

●
●

●

● ●●
●

0.10

●

0.05

●

●

800

1000

1200

1400

1600

Total traveled distance when picking a batch [m], fastest−first(a) VNS batching and fastest-ﬁrst assignment

●
●
●

0.20

●

●

●
●

●

●

●

●

●● ●

●

●
●
●● ●

●

●

●

●
●

●
●
●

●
●

●
●
●
●●

●
●

●

●

●

● ●●

●
●
●
●

●
●

●
●●
●
● ●●
●
●

●

●

●
●

●
●

0.10

●

●
●

●
●

0.05

●

●
●

0.15

●

log(travel) picker model coefficient

●

●

●

●

●

800

1000

1200

1400

1600

Total traveled distance when picking a batch [m], ALNS(b) ALNS
Figure 6.10. The travel distance to pick a batch assigned to a picker vs. the corresponding ln(T ravel) forecasting model coefﬁcient of that picker, for two different
ways of solving the BatchGAP.

99

0.35

Assigning Orders and Batches to Pickers of Varying Skill

0.30

●

●

●

0.25

●
●

●
●

●
●

●

●

●

0.20
0.15

●

●
●

●
●
●

●

●

●
●

●

●

●

●

●

0.10

●
●
●
● ● ● ●● ● ●
●
●
●
● ●
●
●
●
●
●

● ● ●

●

●

●
●
●

●

●

●
●
●

●

●
●
●●

●

●
●

●

●

●

●
●

●
●

●

●

0.05

log(mass) picker model coefficient

●

0.00

●

●

100

200

300

400

Total batch mass [kg], fastest−first

0.35(a) VNS batching and fastest-ﬁrst assignment

0.30

●

●

●

0.25

●

●
●

●

●

●●
●

●●

●

●

●
● ●

●

●

●

0.15

0.20

●

●

●
●
●

●

●
●

0.10

●

●
●
●

●
●● ● ●

●

●
●

●
●
●

●
●

●
●

●

●

●

●
●

●
●

●

●

●
●

●● ●
●

●
●

●

●●
●

●
●

0.05

log(mass) picker model coefficient

●

●
●

0.00

●

●

100

200

300

400

Total batch mass [kg], ALNS(b) ALNS
Figure 6.11. The average mass of a batch assigned to a picker vs. the correspondingln(M ass) forecasting model coefﬁcient of that picker, for two different waysof solving the BatchGAP.

100

Assigning Orders and Batches to Pickers of Varying Skill

●
●

●

0.0

●

●

●

●

●
●

●

●

●
●

●

●

●

−0.5

●
●

●

●

●
●●

●
●●
●

●

●
●
●

●
●

●

●

●
● ●

●●

●

●

●
●

●

● ●

●●
●

●
●
●

●
●

●

●

●

●

●

●

●
●
●

●

●

●

●

●

●

●

●
●

●

●

●

−1.0

log(lvl) picker model coefficient

●
●

●

1.00

1.05

1.10

1.15

1.20

1.25

1.30

1.35

Mean pick level of a batch, fastest−first(a) VNS batching and fastest-ﬁrst assignment

●

0.0

●●
●

●
●●
●

●

●

●
●

−0.5

●

●
●●

●

●

●

●

●

●●
●
●
●
●

● ●

●
●

●
●●
●●●

● ●

●
●●

●

●
●

●

●

●
●

●
●
●

●
●

●

●

●

● ●
●
●
●

●
●●
●

●
●
●

●
●

●
●

●

●
●

−1.0

log(lvl) picker model coefficient

●

●

1.00

1.05

1.10

1.15

1.20

1.25

1.30

1.35

Mean pick level of a batch, ALNS(b) ALNS
Figure 6.12. The mean pick level of a batch assigned to a picker vs. the correspondingln(Level) forecasting model coefﬁcient of that picker, for two different waysof solving the BatchGAP.

101

Assigning Orders and Batches to Pickers of Varying Skill

0.0

●

●●
●

●

●●

●
●

●

●

●

●

●

●

●

●

● ●

●
●
●
●

−0.3

● ●

●

●
●

●

●
●

●

●

●

●

−0.2

●

●

●

●

●
●

●●
●
●
●●

●

●
● ● ●●
●
●
●

●
●
●

●
●

●
●●
●

●

●

●

●●
●
●

●

●

●

●

●

−0.4

Volume picker model coefficient

−0.1

●
●

●

0.6

0.8

1.0

1.2

1.4

1.6

Total volume of a batch, fastest−first(a) VNS batching and fastest-ﬁrst assignment

0.0

●

●

●
● ●
● ●●

●
●

●

●

●

●

● ●

●

●

●

−0.2

● ●●

●

●
●

●

●

●

●

●

●

●

● ●

●

●
●

●
●

●

●

●

●

●

●

●

●
●

●●

●

●

●

●
●
●

●

●
● ●

●

●
●

●
●

−0.3

●

●

●
●
●

●
●

●

●

●
●

−0.4

Volume picker model coefficient

−0.1

●

●

0.6

0.8

1.0

1.2

1.4

1.6

Total volume of a batch, ALNS(b) ALNS
Figure 6.13. The average volume of a batch assigned to a picker vs. the corresponding
V ol forecasting model coefﬁcient of that picker, for two different ways ofsolving the BatchGAP.

102

Assigning Orders and Batches to Pickers of Varying Skill

6.5

Discussion

This chapter introduces a method to model and solve the combined batching, routing, and picker assignment problem, based on properties of thebatches to be executed. The method consists of two parts: ﬁrst the batch
execution times are estimated with multilevel analysis based on properties of historical batches, such as total number of lines, distance between items, total mass and volume, and pick level. Second, the integrated batching, routing, and picker assignment problem is solved using
an ALNS heuristic.
The method is tested for a large retail warehouse for which three monthsof complete pick data could be obtained. In Chapter 5, multilevel modeling of the pickers shows that 13.1% of the variance in batch execution
time can be explained by differences among pickers. Subsequent application of the ALNS heuristic shows that improvements of 9% in total batch
execution time can be achieved compared to state-of-the art batching with
VNS, but ignoring picker skills. Compared to the current picker assignment, 12% is saved depending on the instance (see ﬁgures 6.1 and 6.2).
The presented method is generally applicable in picker-to-parts warehouses, independent of the layout. However, the choice of the routing
heuristic should be made according to the case. The main requirement
for applying the method is the availability of historical pick data to build
the forecasting models of the pickers. These data are gathered in mostmodern warehouses, by pick-by-voice, pick-by-light, and pick-by-terminal
systems. Using these data, warehouse management systems can easily be
extended by including forecast pick times in the batching and assignmentof batches.
Although the method was tested in only one warehouse, it is expectedthat savings can be achieved in other picker-to-parts picking facilities.
The magnitude of the savings will depend on several factors, such as thenumber of order lines in a batch (more lines yield a higher saving opportunity), the size of the pick force (more workers means more opportunity forsavings), the product characteristics (larger variety in sizes and volume
requires more stacking skills of the workers and potentially larger saving opportunities), and variation in the height levels at which the picking
takes place.
This approach has some limitations. The data obtained are extensive,
but they do not contain the departure times of trucks shipping the orders.

103Assigning Orders and Batches to Pickers of Varying Skill

Therefore, there is no information about time windows for the orders,which are to some extent included in the real batch data. This means thecalculated savings slightly overestimate the real savings as the search
space for selecting the next orders to the batch shrinks. However, this
only affects savings compared to the original batching and allocation.
The results presented here suggest that the impact of the right worker
is over 9% in execution time, or 6% if picker productivity is taken into account when assigning work. Furthermore, it is also shown that the skill ofa picker in a particular category inﬂuences the type of batches he or she isassigned: more skilled pickers tend to get more challenging batches to execute. These results imply that managers should consider employee skills
and capabilities during the decision making process. Currently worker
to work assignment is primarily based on qualiﬁcations. However, when
multiple workers qualify for the job, individual differences might still be
exploited. Currently in most warehouses, when workers are qualiﬁed forthe same work, the next job is assigned to the ﬁrst available worker. It
may pay off to select the right job from the stack to be executed by thebest skilled worker. In addition, when hiring a ﬂex workforce and past
performance data are available, the warehouse manager can use them tohire the best performing pickers given a set of orders.

1047. Conclusions

This thesis has two main contributions to the current state-of-the-art
in warehousing and operations research. First, precedence-constrained
routing is considered in conjunction with the order batching problem.
This type of routing introduces new requirements for the order batching algorithm, as it is computationally heavy. It is possible to estimate
savings of batches composed of precedence-constrained customer orders
quite accurately. In other words, the estimate is able to predict how
beneﬁcial it is to combine sets of cities into a larger Traveling Salesman
Problem. This is the key contribution of Chapter 3, as it leads to a dramatic reduction in calls to routing algorithms, which are integral subroutines of most batching algorithms. The algorithm proposed in Chapter 3,
Precedence-Constrained Estimated Savings-based batching, compares favorably in solution quality and execution time to state-of-the-art batchingalgorithms and to optimal solutions for small instances. As seen in Chapter 4, it is also possible to swap the routing algorithm used for anotherone.
The second main contribution is the study of modeling pickers and assigning the right orders and batches to the right pickers, presented in
Chapters 5 and 6. To the author’s knowledge, it is the ﬁrst work thatmodels directly how job parameters affect workers’ job execution times(the worker skills) and exploits this knowledge directly when assigningwork. The main goal of this study is to show that it is beneﬁcial to consider
picker skills in order picking models. A simulation experiment shows thata total of 9% in total order picking can be saved if the order pickers’ skills
are taken into account in conjunction with batching orders.
This study presents multiple opportunities for further research. First,
the results should be further validated in warehouses of different types,
with varying layouts, order, batch and product properties, and those that105

Conclusions

use incentives for pickers as motivation. To form realistic picker proﬁles,data logs need to be obtained, which can be difﬁcult to obtain from companies. The data can be noisy and contain outliers. Additionally, as is
the case with this study, the data can have multiple picking modes, which
may need to be handled separately.
Second, a straightforward extension to this method can be done with theinclusion of time windows. This would require the inclusion of additional
constraints in the model and modiﬁcations to the neighborhood searchheuristics. Another way of extending the method is to consider a rolling
horizon setting, where orders and batches are assigned online to pickers.
Third, although the models of the pickers explain 86% of the variance,
the forecasts might be made more accurate by including other behavioral
variables. Based on previous studies (Bendoly et al., 2006, De Koster
et al., 2011, De Vries et al., 2013), it is likely that also behavioral variables
may have an additional effect. These can include motivational variables
such as prevention and promotion focus, or personality traits.
Fourth, insights obtained from this study can be applied to other areas in addition to picker-to-parts order picking. Parts-to-picker systems,
where products are picked by machines and delivered to the right picker
for further processing, are a natural extension. In general, any job typethat can be parameterized and for which the execution times can be accurately forecast based on these parameters for the different workers, qualiﬁes for careful assignment. The classical vehicle routing problem is an
example. Parcel carriers have to make many deliveries in different typesof areas, urban as well as rural. Based on the types of deliveries, neighborhood, number of parcels, distances, and weight, etc., it might pay off toselect the best driver to ﬁt the job circumstances.

106

Bibliography

Ahuja, R. K., Orlin, J. B., Sharma, D., et al. (1998). New neighborhood search
structures for the capacitated minimum spanning tree problem. Sloan Schoolof Management, Massachusetts Institute of Technology.
Albareda-Sambola, M., Alonso-Ayuso, A., Molina, E., and De Blas, C. (2009).
Variable neighborhood search for order batching in a warehouse. Asia-Paciﬁc
Journal of Operational Research, 26(05):655–683.
Bartholdi, J. J. and Eisenstein, D. D. (1996). A production line that balances
itself. Operations Research, 44(1):21–34.
Bartholdi, J. J., Eisenstein, D. D., and Foley, R. D. (2001). Performance of bucket
brigades when work is stochastic. Operations Research, 49(5):710–719.
Baumann, H. (2013). Order picking supported by mobile computing. PhD thesis,
University of Bremen.
Beamon, B. M. (1998). Supply chain design and analysis: Models and methods.
International journal of production economics, 55(3):281–294.
Bendoly, E., Donohue, K., and Schultz, K. L. (2006). Behavior in operations management: Assessing recent ﬁndings and revisiting old assumptions. Journalof Operations Management, 24(6):737–752.
Berger, S. M. and Ludwig, T. D. (2007). Reducing warehouse employee errors
using voice-assisted technology that provided immediate feedback. Journal of
Organizational Behavior Management, 27(1):1–31.
Bliese, P. D. (2002). Multilevel random coefﬁcient modeling in organizational
research: Examples using SAS and S-PLUS. Jossey-Bass.
Bowersox, D. J. (2011). Supply chain logistics management. Tata McGraw-Hill
Education.
Bryk, A. S. and Raudenbush, S. W. (1992). Hierarchical linear models: Applications and data analysis methods. Sage Publications, Inc.
Campbell, G. M. and Diaby, M. (2002). Development and evaluation of an assignment heuristic for allocating cross-trained workers. European Journal ofOperational Research, 138(1):9–20.
Cattrysse, D. G. and Van Wassenhove, L. N. (1992). A survey of algorithms for thegeneralized assignment problem. European Journal of Operational Research,
60(3):260–272.

107

Bibliography

Chan, F. and Kumar, V. (2008). A TSSA algorithm based approach to enhance theperformance of warehouse system. In 10th IEEE International Conference onControl, Automation, Robotics and Vision, 2008. ICARCV 2008., pages 1696–
1701.
Clarke, G. and Wright, J. (1964). Scheduling of vehicles from a central depot toa number of delivery points. Operations Research, 12(4):568–581.
Cohn, H. and Fielding, M. (1999). Simulated annealing: Searching for an optimal
temperature schedule. SIAM Journal on Optimization, 9(3):779.
Cordeau, J.-F., Laporte, G., Pasin, F., and Ropke, S. (2010). Scheduling technicians and tasks in a telecommunications company. Journal of Scheduling,
13(4):393–409.
Cornuéjols, G., Fonlupt, J., and Naddef, D. (1985). The traveling salesman problem on a graph and some related integer polyhedra. Mathematical Programming, 33(1):1–27. Available from: http://dx.doi.org/10.1007/BF01582008.
Coyle, J., Bardi, E., and Langley, C. (1996). The Management of Business Logistics. West Publishing Company.
De Koster, R., Le-Duc, T., and Roodbergen, K. (2007). Design and control of warehouse order picking: A literature review. European Journal of Operational
Research, 182(2):481–501.
De Koster, R., Stam, D., and Balk, B. M. (2011). Accidents happen: The inﬂuence of safety-speciﬁc transformational leadership, safety consciousness, andhazard reducing systems on warehouse accidents. Journal of Operations Management, 29(7):753–765.
De Koster, R. and Van der Poort, E. (1998). Routing order pickers in a warehouse: a comparison between optimal and heuristic solutions. IIE Transactions, 30(5):469–480.
De Koster, R., Van der Poort, E., and Wolters, M. (1999). Efﬁcient orderbatching methods in warehouses. International Journal of Production Research,
37(7):1479–1504.
De Vries, J., De Koster, R., and Stam, D. (2013). Safety does not happen by accident: How to manage a safe warehouse? Working paper, Erasmus University,
Rotterdam.
Dekker, R., De Koster, R., Roodbergen, K., and Van Kalleveen, H. (2004). Improving order-picking response time at Ankor’s warehouse. Interfaces, 34(4):303–
313.
Doerr, K. H. and Arreola-Risa, A. (2000). A worker-based approach for modelingvariability in task completion times. IIE Transactions, 32(7):625–636.
Drury, J. (1988). Towards more efﬁcient order picking. IMM Monograph 1, The
Institute of Materials Management, Cranﬁeld, U.K.
Duan, N. (1983). Smearing estimate: a nonparametric retransformation method.
Journal of the American Statistical Association, 78(383):605–610.
Escudero, L. (1988). An inexact algorithm for the sequential ordering problem.
European Journal of Operational Research, 37(2):236–249.

108

Bibliography

Fisher, M. L., Jaikumar, R., and Van Wassenhove, L. N. (1986). A multiplier adjustment method for the generalized assignment problem. Management Science, 32(9):1095–1103.
Gademann, N. and Van de Velde, S. (2005). Order batching to minimize total
travel time in a parallel-aisle warehouse. IIE Transactions, 37(1):63–75.
Gelman, A. and Hill, J. (2007). Data analysis using regression and multilevel/hierarchical models. Cambridge University Press.
Gharehgozli, A., Laporte, G., Yu, Y., and De Koster, R. (2013). Scheduling twin
yard cranes in a container block. Transportation Science, forthcoming.
Goetschalckx, M. and Ashayeri, J. (1989). Classiﬁcation and design of order picking. Logistics Information Management, 2(2):99–106.
Gu, J., Goetschalckx, M., and McGinnis, L. F. (2007). Research on warehouse operation: A comprehensive review. European Journal of Operational Research,
177(1):1–21.
Hart, P., Nilsson, N., and Raphael, B. (1968). A formal basis for the heuristic
determination of minimum cost paths. IEEE Transactions on Systems Science
and Cybernetics, 4(2):100–107.
Helsgaun, K. (2000). An effective implementation of the Lin-Kernighan travelingsalesman heuristic. European Journal of Operational Research, 126(1):106–
130.
Henn, S. and Wäscher, G. (2012). Tabu search heuristics for the order batchingproblem in manual order picking systems. European Journal of Operational
Research.
Hong, S., Johnson, A., and Peters, B. (2012a). Batch picking in narrow-aisle orderpicking systems with consideration for picker blocking. European Journal ofOperational Research.
Hong, S., Johnson, A., and Peters, B. (2012b). Large-scale order batching inparallel-aisle picking systems. IIE Transactions, 44(2):88–106.
Hsieh, L. and Huang, Y. (2011). New batch construction heuristics to optimise
the performance of order picking systems. International Journal of Production
Economics, 131(2):618–630.
Juran, D. C. and Schruben, L. W. (2004). Using worker personality and demographic information to improve system performance prediction. Journal ofOperations Management, 22(4):355–367.
Kohonen, T. (1990).
78(9):1464–1480.

The self-organizing map.

Proceedings of the IEEE,

Kreft, I. and De Leeuw, J. D. (1998). Introducing Multilevel Modeling. Sage,
London, UK.
Kubo, M. and Kasugai, H. (1991). The precedence constrained traveling salesman problem. Journal of the Operations Research Society of Japan, 34(2):152–
172.

109

Bibliography

Kuijt-Evers, L., Bosch, T., Huysmans, M., De Looze, M., and Vink, P. (2007).
Association between objective and subjective measurements of comfort anddiscomfort in hand tools. Applied Ergonomics, 38(5):643–654.
Larco Martinelli, J. (2010). Incorporating Worker-Speciﬁc Factors in Operations
Management Models. PhD thesis, Erasmus Research Institute of Management, Erasmus University Rotterdam. Available from: http://hdl.handle.
net/1765/21527.
Lundy, M. and Mees, A. (1986). Convergence of an annealing algorithm. Mathematical Programming, 34(1):111–124.
MacQueen, J. (1967). Some methods for classiﬁcation and analysis of multivariate observations. In Proceedings of the ﬁfth Berkeley symposium on mathematical statistics and probability, volume 1, page 14. California, USA.
Matusiak, M., de Koster, R., Kroon, L., and Saarinen, J. (2014). A fast simulated annealing method for batching precedence-constrained customer orders in a warehouse. European Journal of Operational Research, 236(3):968– 977. Vehicle Routing and Distribution Logistics. Available from: http://www.sciencedirect.com/science/article/pii/S0377221713004670.
Miller, C., Tucker, A., and Zemlin, R. (1960). Integer programming Formulationof Traveling Salesman Problems. Journal of the A.C.M., 7:326–329.
Nakagawa, S. and Schielzeth, H. (2013). A general and simple method for obtaining R2 from generalized linear mixed-effects models. Methods in Ecology
and Evolution, 4(2):133–142. Available from: http://dx.doi.org/10.1111/j.
2041-210x.2012.00261.x.
Pentico, D. W. (2007). Assignment problems: A golden anniversary survey. European Journal of Operational Research, 176(2):774–793.
Pisinger, D. and Ropke, S. (2010). Large neighborhood search. In Handbook of
metaheuristics, pages 399–419. Springer.
Powell, S. G. and Schultz, K. L. (2004). Throughput in serial lines with statedependent behavior. Management Science, 50(8):1095–1105.
Psaraftis, H. (1980a). A dynamic programming approach for sequencing groupsof identical jobs. Operations Research, 28(6):1347–1359.
Psaraftis, H. (1980b). Dynamic programming solution to the single vehicle
many-to-many immediate request dial-a-ride problem. Transportation Science, 14(2):130–154.
Psaraftis, H. (1983). An exact algorithm for the single vehicle many-to-many
dial-a-ride problem with time windows. Transportation Science, 17(3):351–
357.
Randolph, W. (1993). Distance approximations for routing manual pickers in awarehouse. IIE Transactions, 25(4):76–87.
Ratliff, H. and Rosenthal, A. (1983). Order-picking in a rectangular warehouse:a solvable case of the traveling salesman problem. Operations Research,
31(3):507–521.

110

Bibliography

Roodbergen, K. and De Koster, R. (2001a). Routing methods for warehouses
with multiple cross aisles. International Journal of Production Research,
39(9):1865–1883. Available from: www.tandfonline.com.
Roodbergen, K. and De Koster, R. (2001b). Routing order pickers in a warehousewith a middle aisle. European Journal of Operational Research, 133(1):32–43.
Ropke, S. and Pisinger, D. (2006a). An adaptive large neighborhood searchheuristic for the pickup and delivery problem with time windows. Transportation Science, 40(4):455–472.
Ropke, S. and Pisinger, D. (2006b). A uniﬁed heuristic for a large class of vehicle
routing problems with backhauls. European Journal of Operational Research,
171(3):750–775.
Shaw, P. (1997). A new local search algorithm providing high quality solutions tovehicle routing problems. APES Group, Dept of Computer Science, Universityof Strathclyde, Glasgow, Scotland, UK.
Theys, C., Braysy, O., Dullaert, W., and Raa, B. (2010). Using a TSP heuristic
for routing order pickers in warehouses. European Journal of Operational
Research, 200(3):755–763.
Tompkins, J., White, J., Bozer, Y., and Tanchoco, J. (2003). Facilities Planning,
3rd ed. John Wiley & Sons.
Vaughan, T. (1999). The effect of warehouse cross aisles on order picking efﬁciency. International Journal of Production Research, 37(4):881–897.
Weaver, K. A., Baumann, H., Starner, T., Iben, H., and Lawo, M. (2010). An empirical task analysis of warehouse order picking using head-mounted displays.
In Proceedings of the SIGCHI Conference on Human Factors in Computing
Systems, pages 1695–1704. ACM.
Whittley, I. and Smith, G. (2004). The attribute based hill climber. Journal of
Mathematical Modelling and Algorithms, 3(2):167–178.

111

School of Electrical Engineering

Marek Matusiak

Aalto-DD 62/2014

Optimizing Warehouse Order Batching when Routing Is Precedence Constrained and Pickers Have
Varying Skills

Aalto University
School of Electrical Engineering
School of Electrical Engineeringwww.aalto.fi

Marek Matusiak

BUSINESS +ECONOMYART +DESIGN +ARCHITECTURE
SCIENCE +TECHNOLOGYCROSSOVER
DOCTORALDISSERTATIONS

Aalto University

9HSTFMG*afgheg+

ISBN 978-952-60-5674-6
ISBN 978-952-60-5675-3 (pdf)ISSN-L 1799-4934
ISSN 1799-4934
ISSN 1799-4942 (pdf)

Optimizing Warehouse
Order Batching when
Routing Is Precedence
Constrained and Pickers
Have Varying Skills

DOCTORALDISSERTATIONS