The Stochastic Container Relocation Problem
V. Galle ∗

S. Borjian Boroujeni †

V. H. Manshadi ‡

C. Barnhart §

P. Jaillet ¶

arXiv:1703.04769v2 [cs.DS] 12 Oct 2017

Abstract
The Container Relocation Problem (CRP) is concerned with finding a sequence of moves of containers that minimizes the number of relocations needed to retrieve all containers, while respecting a
given order of retrieval. However, the assumption of knowing the full retrieval order of containers
is particularly unrealistic in real operations. This paper studies the stochastic CRP (SCRP), which
relaxes this assumption. A new multi-stage stochastic model, called the batch model, is introduced,
motivated, and compared with an existing model (the online model). The two main contributions are an
optimal algorithm called Pruning-Best-First-Search (PBFS) and a randomized approximate algorithm
called PBFS-Approximate with a bounded average error. Both algorithms, applicable in the batch
and online models, are based on a new family of lower bounds for which we show some theoretical
properties. Moreover, we introduce two new heuristics outperforming the best existing heuristics.
Algorithms, bounds and heuristics are tested in an extensive computational section. Finally, based on
strong computational evidence, we conjecture the optimality of the “Leveling” heuristic in a special
“no information” case, where at any retrieval stage, any of the remaining containers is equally likely
to be retrieved next.

Introduction
With the growth in international container shipping in maritime ports, there has been an increasing
interest in improving operations in container terminals, both on the sea side and land side. The operations
on the sea side involve the assignment of quay cranes to ships, the loading of export containers on vessels,
and the discharging of import containers from vessels onto internal trucks. Import containers are then
transferred to the land side and are stacked in the storage area. Operations on the land side (also called
yard operations) include the routing of internal trucks within the yard, the stacking of containers for
storage, and the delivery of import containers to external trucks for delivery to another location. This
work focuses on the latter problem.
Due to limited space in the storage area, containers are stacked on top of each other. The resulting
stacks create rows of containers as shown in Figure 1. If a container that needs to be retrieved (target
container) is not located at a top most tier and is covered by other containers, the blocking containers
must be relocated to another stack. As a result, during the retrieval process, one or more relocation
moves are performed by the yard cranes. Such relocations (also called reshuffles) are costly for the port
∗

Operations Research Center. Email: vgalle@mit.edu.
Operations Research Center. Email: setareh@alum.mit.edu
‡
Yale School of Management. Email: vahideh.manshadi@yale.edu
§
Operations Research Center and Civil & Environmental Engineering. Email: cbarnhar@mit.edu
¶
Operations Research Center and Electrical Engineering & Computer Science. Email: jaillet@mit.edu
†

1

Figure 1: Illustration of stacks of containers in a storage yard (figure from Tanaka and Takii (2014))
operators and result in delays in the retrieval process. Thus, reducing the number of relocations is one
of the main goals of port operators. The Container Relocation Problem (CRP) (also known as the Block
Relocation Problem) addresses this challenge by minimizing the number of relocations. As this problem
is the main discussion of this paper, we provide a formal definition and an extensive literature review.
The container relocation problem First, it is commonly known that the time to relocate a container
within a row is insignificant compared to the time to relocate a container between two distinct rows.
Therefore, in most cases, port operators tend to avoid relocations between rows. The CRP makes the
assumption that only relocations within rows are allowed, and problems for different rows should be
considered independently. Furthermore, a row usually stores containers of the same type for the sake of
stability and simplicity.
Using these facts, the CRP models one row using a two dimensional array of size (T, S), where S is
the number of stacks, and T is the maximum height, i.e., the maximum number of containers in a stack
limited by the height of the crane. Each element of this array represents a potential slot for a container,
and it contains a number only if a container is currently stored in this slot. Stacks are numbered from left
(1) to right (S) and tiers from bottom (1) to top (T). We refer to this array as a configuration. The common
assumptions of the CRP are the following:
A1 : The initial configuration has T tiers, S stacks, and C containers. In order for the problem to
always be feasible, we suppose that the triplet (T, S, C) satisfies 0 6 C 6 ST − (T − 1).
A2 : A container can only be retrieved/relocated if it is at the top most tier of its stack, i.e., no other
container is blocking it.
A3 : A container can only be relocated if it is blocking the target container. This assumption was
suggested by Caserta et al. (2012), and the problem under this assumption is commonly referred to
as the restricted CRP. Most studies focus on this restricted version, because it is the current practice
in many yards, and it helps decrease the dimensionality of the problem, while not losing much

2

optimality (see Petering and Hussein (2013)). As is common practice, we will not mention the term
“restricted” in the rest of the paper even though we always assume A3 .
A4 : The cost of relocating a container from a stack does not depend on, to which stack the container
is relocated. This allows us to consider the stacks of a configuration as interchangeable. In addition,
it motivates the objective of minimizing the number of relocations, since the cost of each relocation
can be normalized to 1. Note that this assumption is not required for all the results stated below,
hence our approaches could be easily extended to the case when Assumption A4 does not hold.
A5 : The retrieval order of containers is known, so that each container can be labeled from 1 to C,
representing the departure order, i.e., Container 1 is the first one to be retrieved, and C the last one.
The CRP involves finding a sequence of moves to retrieve Containers 1, 2, . . . , C (respecting the order)
with a minimum number of relocations. Figure 2 provides a simple example of the CRP. The CRP with the
above classical assumptions is referred to as static and full information: “Static” because no new containers
arrive during the retrieval process (see Assumption A1 ) and “full information” because we know the
full retrieval order at the beginning of the retrieval process (see Assumption A5 ). This problem was first
formulated by Kim and Hong (2006) in a dynamic programming model.
2
3

6

4

5

1

Reloc 2

−−−−−→

2

6

4

3

5

1

Reloc 4

−−−−−→

4
2

6

3

5

Ret 1,2,3,4

−−−−−−−→
1

6
5

Reloc 6

−−−−−→
6

5

Figure 2: Configuration for the CRP with 3 tiers, 3 stacks and 6 containers. The optimal solution performs
3 relocations: relocate the container labeled 2 from Stack 3 to Stack 1 on the top of the container labeled
3; relocate 4 from 3 to 2 on the top of 6; retrieve 1; retrieve 2; retrieve 3; retrieve 4; relocate 6 from 2 to the
empty Stack 1; retrieve 5; finally, retrieve 6.
Researchers have tackled the static CRP with full information from two point of views. The first
approach aims to find the optimal solution. Primarily, researchers have used Integer Programming
(IP) to address this problem. For example, Caserta et al. (2012) propose an intuitive formulation of the
problem. Petering and Hussein (2013) develop a more tractable formulation, that is, however, unable
to solve real-sized instances efficiently. Zehendner et al. (2015) fix the formulation from Caserta et al.
(2012) and improve it by removing some variables, tightening some constraints, introducing a new upper
bound, and applying a pre-processing step to fix several variables. In all these IP formulations, due to
the combinatorial nature of the problem, the number of variables and constraints dramatically increases
as the size of the bay grows, and the IP cannot be solved for large instances. In order to bypass this
problem, a recent trend has been to look at more efficient ways to explore the branch-and-bound tree,
or even decrease its size using the structural properties of the problem. Ünlüyurt and Aydın (2012) and
Expósito-Izquierdo et al. (2015) suggest two branch-and-bound approaches with several heuristics based
on this idea. Another solution using the A∗ algorithm is explored by Zhu et al. (2012), and built upon by
Tanaka and Takii (2014) and Borjian et al. (2015a). Another solution using branch-and-price is presented
by Zehendner and Feillet (2014b).
As the problem is NP-hard (Caserta et al. (2012)), an alternative approach is to use quick and efficient
heuristics providing sub-optimal solutions. For the sake of conciseness, we only mention some of them
3

that are relevant to this paper. Caserta et al. (2012) introduce a ’MinMax’ policy that is defined and
generalized later in this paper. Wu and Ting (2010) propose a beam search heuristic, and Wu and Ting
(2012) develop the Group Assignment Heuristic (GAH) also generalized in Section 2. Finally, in Kim
and Hong (2006) and Zhu et al. (2012), lower bounds for the CRP are introduced: Kim and Hong (2006)
count the number of blocking containers as a straightforward lower bound, and Zhu et al. (2012) refine
this idea by taking into account additional unavoidable relocations using a family of lower bounds.
Finally, there are many related problems to the CRP. The stacking problem is concerned with how to
store incoming containers in a configuration given an arrival order of containers. The pre-marshalling
problem deals with re-arranging the containers prior to the retrieval process in order to minimize
future relocations, but no container is removed in this process. For both problems, general review and
classification surveys of the existing literature on the CRP can be found in Stahlbock and Voß (2008),
Steenken et al. (2004) and Lehnfeld and Knust (2014). In addition, Assumption A1 can be relaxed, and
this leads to the dynamic CRP where stacking and retrieving are done simultaneously as new containers
are arriving. For this problem, see Borjian et al. (2015b) and Akyüz and Lee (2013).
Finally, the main focus of this paper is an extension of the CRP where the full information assumption
(A5 ) is relaxed. Indeed, Assumption A5 is unrealistic given that arrival times of external trucks at
the terminal are generally unpredictable due to uncertain conditions. Nevertheless, new technology
advancements such as Truck Appointment Systems (TAS’s) and GPS tracking can help predict relative
truck arrival times. Thus, although the exact retrieval order might not be known, some information on
trucks’ arrival times might be available. This leads us to introduce a stochastic version of the CRP.
The Stochastic CRP (SCRP) A common assumption is that, for each container, there is a time window
in which a truck driver will arrive to retrieve it. We refer to a batch of containers as the set of containers
stacked in the same row and with the same arrival time window. This information can be either inferred
using machine learning algorithms, not yet much discussed in the literature, or can be obtained using the
appointment time windows in a TAS, which has gained attention over the last decade. The first TAS was
implemented by Hong Kong International Terminal (HIT) in 1997. It uses 30-minute time slots, where
trucks can register (Murty et al. (2005)). Another TAS was introduced in New Zealand in 2007. Two
other studies, Giuliano and O’Brien (2007) and Morais and Lord (2006), evaluate the benefits of TAS, in
reducing truck idling time by increasing on-time ratio. More recent information can be found in Phillips
(2015) and Bonney (2015).
On the modeling side, Zehendner and Feillet (2014a) formulate an IP to get the optimal number of slots
a TAS should offer for each batch. Very few studies have tackled the SCRP, also referred to as CRP with
Time Windows. This problem was first modeled by Zhao and Goodchild (2010). In the original model
of Zhao and Goodchild (2010), each container is assigned to a batch. Batches of containers are ordered
such that all containers in a batch must be retrieved before any containers in a later batch. Furthermore,
the relative retrieval order of containers within a given batch is assumed to be a random permutation.
From now on, we will refer to the model of Zhao and Goodchild (2010) as the online model. In Section
1, we discuss in more detail how this model assumes information is revealed. For the online model,
Zhao and Goodchild (2010) develop a myopic heuristic (called RDH) and study, in different settings with
two or multiple groups, the value of information using RDH. They conclude that a small improvement
in the information system reduces the number of relocations significantly. Van Asperen et al. (2013)
4

use a simulation tool to evaluate the effect of a TAS on many statistics including the ratio of relocations
to retrievals. Their decision rules are based on several heuristics including the “leveling,” random, or
“traveling distance” heuristics. More recently, Ku and Arthanari (2016) also use the online model. They
formulate the SCRP under the online model as a finite horizon dynamic programming problem, and
suggest a decision tree scheme to solve it optimally. They also introduce a new heuristic called ERI
(Expected Reshuffling Index), which outperforms RDH, and they perform computational experiments
based on available test instances. We will refer frequently to this work, use some of their techniques, as
well as their available test instances to evaluate our algorithms.
In another recent study related to the SCRP, Zehendner et al. (2016) study the Online Container
Relocation Problem, which corresponds to an adversarial model. They prove that the number of relocations
performed by the leveling policy can be upper-bounded by a linear function of the number of blocking
containers and provide a tight competitive ratio for this policy. Moreover, Galle et al. (2016) show that
the ratio of the expected number of relocations to the expected blocking lower bound converges to
one. Finally, Tierney and Voß (2016) study the robust pre-marshalling problem which also considers
uncertainty in the retrieval times of containers.
For a general review of techniques on finite horizon Dynamic Programming, we refer the reader to
Bertsekas (2005) and Sennott (2009). Table 1 summarizes the previous literature review.
Contributions of the paper The contributions of this paper are:
1. A new stochastic model, referred to as the batch model. This new model uses the same probability
distribution as the online model. However, the two models are different in the way new information
on the retrieval order is revealed. The batch model is motivated, described and compared with the
online model.
2. Lower and upper bounds for the SCRP. We derive a new family of lower bounds for which we
show theoretical properties. Furthermore, we develop two new fast and efficient heuristics.
3. A novel optimal algorithm scheme based on decision trees and pruning strategies referred to as
Pruning-Best-First-Search (PBFS), taking advantage of the properties of the aforementioned lower
bounds. The algorithm is explained with pseudocode in Algorithm 2.
4. A second novel algorithm tuned for the case of larger batches referred to as PBFSA (PBFSApproximate). We build upon PBFS and derive a sampling strategy resulting in an approximate
algorithm with an expected error that we bound theoretically. The pseudocode of the second
algorithm is presented in Algorithm 3.
5. We provide extensive computational experiments using an existing set of instances. The first
experiment exhibits the efficiency of PBFS, our lower bounds and two new heuristics for the
batch model based on existing instances, presented by Ku and Arthanari (2016), where batches of
containers are small (2 containers per batch on average). The second experiment illustrates the
advantage of using PBFSA when batches of containers are larger, based on instances obtained by
modifying the existing set. In addition, most of our techniques including lower bounds, heuristics,
and the PBFS algorithm also apply to the online model. The third experiment shows that, in this
5

model, PBFS outperforms the algorithm introduced in Ku and Arthanari (2016) in the sense that it
is faster for instances that Ku and Arthanari (2016) could solve, and it can solve problems of larger
size. Furthermore, our two new heuristics also outperform the best existing heuristic (ERI) for the
online model. Finally, the last experiment is used to test the conjecture about the optimality of the
leveling heuristic in the special case of the online model with a unique batch of containers.

CRP

Static

Dynamic

Caserta et al. (2012)
Petering and Hussein (2013)
Full

Zehendner et al. (2015)

Borjian et al. (2015b)

Information

Expósito-Izquierdo et al. (2015)

Akyüz and Lee (2013)

Zhu et al. (2012)
Tanaka and Takii (2014)
Zhao and Goodchild (2010)
Partial

Ku and Arthanari (2016)

Information

This paper

Table 1: Optimal solutions for the different variant of the CRP
The rest of the paper is structured as follows: Section 1 thoroughly describes the batch model, its
assumptions and objective, the difference with the online model, and the general theory of decision trees
applied to the SCRP. Section 2 gives a good intuition into the problem and defines heuristics and a class of
lower bounds for the SCRP used in subsequent sections. Then Sections 3 and 4 introduce respectively the
PBFS and PBFSA algorithms. Computational experiments for both batch and online models are carried
through Section 5. We conclude the paper by discussing future directions for the SCRP in Section 6.

1
1.1

SCRP and Decision Trees
Motivation

Before stating the general assumptions of the batch model, let us motivate our problem using a typical
example. We consider a port with a TAS offering 30-minute time windows during which truck drivers
who want to retrieve a container can register to arrive at the port. For the sake of the example, we
consider the time window between 9:00 am and 9:30 am. Multiple trucks can be registered in this time
window: in this example, presented in Figure 3, 3 trucks (designated i1 , i4 and i6 ) are registered for this
time window. We assume that all 3 trucks arrive on time (between 9:00 am and 9:30 am) and that their
containers (similarly designated i1 , i4 and i6 ) form a batch to be retrieved. We display the configuration
of interest in Figure 4 (3 tiers, 3 stacks and 6 containers).

6

Appointment time window
9:30am

9:00am

New
information
revealed
Truck
arrival

Truck i4
Truck i6

Truck
arrival

Start batch retrieval
using decisions

Truck
arrival

Service Time

Service
Time

Waitingtime
Time
Waiting

Service time

Waitingtime
Time
Waiting

Waitingtime
Time
Waiting

Truck i1

Service
Time

Service time
Service
Time

Service time

Figure 3: Timeline of events for the batch model with three trucks
We assume that trucks arrive randomly within the time window, so each truck arrival order is equally
likely to happen. In this example, there are 6 potential arrival orders, each with 1/6 chance of occurring.
At 9:00 am, none of the 3 trucks has arrived and their relative retrieval order is unknown. Consequently, these 3 containers are all labelled 1 in Figure 4a. In Figure 4b, the IDs of all containers and their
locations are depicted.
Between 9:00 am and 9:30 am, trucks arrive in a particular order (e.g., Truck i4 first, then i6 , and i1 last).
In busy terminals, trucks typically queue up as they wait to be served. Their place in line is based on their
arrival order, so the port operator generally retrieves containers based on the arrival order. Processing
in this way, on a first-come-first-served basis, avoids issues with truck unions and maintains fairness
among drivers. Consequently, we take the retrieval order to be exogenously determined and we do not
consider it a potential decision for port operators.
In order to provide a specified level of service to the truck drivers, the terminal operator often sets a
target average waiting time. If the appointment time window is about the same as or shorter than the
target average waiting time, the operator has information about the retrieval order of containers in the
batch before the retrieval of those containers must begin in order to meet the target waiting time. Given
this, we make the simplifying assumption in this work that the retrieval of a batch begins at the end of
the appointment time window associated with the batch, and the retrieval order of all containers in the
batch is known as the retrieval of the batch commences. In our example, the target average waiting
time is 30 minutes. At 9:30 am, the retrieval order of the batch (i4 , i6 , i1 ) is known and the retrieval of the
batch commences soon after. The updated information is depicted in the configuration of Figure 4c.
The assumption that containers to be retrieved are revealed on a batch basis models the reality that
port operators typically know information about all the containers in the same batch before starting
to retrieve them. This is especially true for busy ports that have a TAS. Moreover, we assume that
no information about future batches is available when making decisions for the current batch. Similar
modeling assumptions have been made in previous works (see Zhao and Goodchild (2010) and Ku and
Arthanari (2016)).
The general assumptions we apply to our model are formally stated in Section 1.2 (A∗5 and A∗6 ) and
7

i6

1
1

5

4

5

1

4a Before any truck has
arrived (9:00 am)

i1

i3

i5

i2

i4

4b IDs to match containers
with trucks

2
3

5

4

5

1

4c Before the first container
gets retrieved (9:30 am)

Figure 4: SCRP Example. The left configuration is the input to our problem. The configuration in the
middle denotes each container with an ID il where l = 1, . . . , 6. The configuration on the right denotes
the order of the first batch after it is revealed.
result in the batch model – the main focus of this paper. The goal of the SCRP is to find a sequence of moves
minimizing the expected number of relocations needed to empty the initial configuration.
Labels in Figures 4a and 4c are defined such that two containers only have the same label if they are
in the same batch and their relative order is yet to be revealed. In our example, since Container i5 is the
only container in the second batch and is retrieved after the first 3 containers, it is necessarily the fourth
container to be retrieved (thus labeled 4). Containers i2 and i3 are labeled 5 and when their relative order
is revealed, one will be labeled 5 and the other 6.
The online model The main distinguishing difference between the batch model and the online model
introduced by Zhao and Goodchild (2010) is the information revelation process. The online model
disregards any within batch information available when it plans the moves to retrieve a container. Hence
new information is revealed one container at a time. The online model is especially applicable in
less busy ports where the waiting time is significantly shorter than the appointment time windows. In
this case, because there is a limited number of trucks waiting, limited information about the future is
known. We show through Lemma 1 that ignoring information (if available) results in a potential loss
of operational efficiency as measured by the expected number of relocations. In addition, most of our
batch-based approaches also apply to the online model, and we provide numerical results based on it in
Section 5.

1.2

Assumptions, notations and formulation

In order to define our problem as a multi-stage stochastic optimization problem (the number of
“stages” is the number of batches), we need to define a probabilistic model of the container retrieval order, a
scheme for revealing new information about this order, and an objective function.
Batch model Let us state the assumptions and objective of the stochastic CRP under the batch model.
Assumptions A∗1 , A∗2 , A∗3 and A∗4 are respectively identical to Assumptions A1 , A2 , A3 and A4 .
A∗5 : (probabilistic model) Given an ordering of batches, the probability distribution of the retrieval
order is such that: 1) all the containers in a given batch are retrieved before any containers in
a later batch, and 2) within each batch of container, the order of the containers is drawn from
a uniform random permutation. This paper focuses mainly on this latter assumption, but our

8

model can be extended to the more general case of any probability distribution on permutations
(not necessarily uniform) that respects the order of batches.
A∗6 : (revelation of new information) For each batch w, the full relative order of containers from
the wth batch (i.e., the specific random permutation drawn) is revealed after all containers in
batch 1 through w − 1 have been retrieved.
Given these assumptions, we want to find the minimum expected number of relocations to retrieve
all containers from a given initial configuration. We refer to this problem as the “Stochastic CRP”. Let us
introduce some notations:
- The problem size is given by (T, S, C), respectively the number of tiers, stacks and containers in the
initial configuration.
- The number of batches of containers in the initial configuration is denoted by W. We consider that
the batches are ordered from 1 to W.
- For each batch w ∈ {1, . . . , W}, let Cw be the number of containers in w. By definition

W
X

Cw = C.

w=1

- Each container has two attributes:
• The first attribute, denoted by (cl )l∈{1,...,C} , is the label and is defined as follows: initially,
w−1
X
containers in batch w are labeled by Kw such that Kw = 1 +
Cu , where the sum is empty for
u=1

w = 1. Then, for k ∈ {1, . . . , C}, if a container is revealed to be the kth container to be retrieved,
its label changes to k. Using this labeling, at any point in the retrieval process, two containers
only have the same label if there are in the same batch and their relative order is yet to be
revealed.
• The second attribute is a unique ID denoted by (il )l∈{1,...,C} . This ID is only used to identify
uniquely containers in the initial configuration (see Figure 4b) and for the sake of clarity of
the following probabilistic model. Note that for l ∈ {1, . . . , C}, if Container il is in batch w, then
cl = Kw (until the actual retrieval order of il is revealed).
- For k ∈ {1, . . . , C}, let ζk be a random variable taking values in (il )l∈{1,...,C} , such that {ζk = il } is
the event that Container il is the kth container to be retrieved. According to Assumption A∗5 , the
distribution of (ζk )k∈{1,...,C} is given by


1

 Cw , if w = min{u ∈ {1, . . . , W} | Ku > k} and cl = Kw
P [ζk = il ] = 

 0 , otherwise
Q
In this case, there are a total of W
w=1 (Cw !) orders with equal probabilities. More generally, we
consider the case where the probability of each order within each batch is not necessarily equally
likely. However, we still assume that the batches are ordered, thus P [ζk = il ] can be positive only
if w = min{u|Ku > k} and cl = Kw . In the practical case where probabilities are not considered to be
uniform, a list of potential retrieval orders and their associated probability is given for each batch of
containers (based on historical data), thus P [ζk = il ] can easily be inferred from these probabilities.
9

- An action corresponds to a sequence of relocations to retrieve one container. For k ∈ {1, . . . , C}, we
denote the action for the kth retrieval by ak , and the feasible set of actions is defined according to
Assumptions A∗2 and A∗3 .
- For a given batch w ∈ {1, . . . , W},
1. let yw denote the configuration before the retrieval order of containers in batch w is revealed,
i.e., before ζKw , . . . , ζKw +Cw −1 are revealed. Note that y1 corresponds to the initial configuration.
We denote xKw the configuration after the retrieval order of containers in batch w is revealed,
ζKw ,...,ζKw +Cw −1

and before action aKw is taken. If −−−−−−−−−−−→ represents the revelation of the random variables
ζKw ,...,ζKw +Cw −1

ζKw , . . . , ζKw +Cw −1 , we can write yw −−−−−−−−−−−→ xKw .
2. After the retrieval order of batch w has been revealed, actions to retrieve the revealed containers
have to be made. If Cw > 1, then {Kw , . . . , Kw + Cw − 2} , ∅. In this case, for k ∈ {Kw , . . . , Kw +
Cw − 2}, let xk+1 be the configuration after applying action ak to state xk and before action ak+1 .
ak
ak
Therefore, if −→ represents the application of action ak , we have xk −→ xk+1 .
3. The last container to be retrieved in batch w is the (Kw + Cw − 1)-th container, thus, according
to the previous point, xKw +Cw −1 corresponds to the configuration before aKw +Cw −1 is taken. After
this retrieval, the order of the next batch (batch w + 1) has to be revealed, and according to the
first point, the configuration is yw+1 . The configuration after retrieving batch W will be empty,
aKw +Cw −1
thus we define yW+1 to be the empty configuration. So if −−−−−−→ represents the application
aKw +Cw −1
of action aKw +Cw −1 , then we have xKw +Cw −1 −−−−−−→ yw+1 .
In summary, we have

ζKw ,...,ζKw +Cw −1


yw −−−−−−−−−−−→ xKw




ak
∀ w ∈ {1, . . . , W} , 
xk −→ xk+1 , if Cw > 1 , ∀ k ∈ {Kw , . . . , Kw + Cw − 2}




Kw +Cw −1
 xK +C −1 −a−
−−−−→ yw+1 .
w
w
- Let the function r(.) be such that r(x) is number of relocations to retrieve the target container in
configuration x. It is also equal to the number of containers blocking the target container. This
function is only defined for configurations in which the target container is identified. Specifically, it

is defined for (xk )k=1,...,C (but not for yw w=1,...,W ). For k ∈ {1, . . . , C}, we refer to r(xk ) as the immediate
cost for the kth retrieval.
- Let the function f (.) be such that f (x) is the minimum expected number of relocations required to
retrieve all containers from configuration x. f (x) is commonly referred to as the cost-to-go function

of configuration x. Note that it is well defined for both (xk )k=1,...,C and yw w=1,...,W .
By definition, we have:






f yw =
E
f (xk )


ζk ,...,ζk+Cw −1



 f (x ) = r (x ) + min  f (x )
∀ w ∈ {1, . . . , W}, 
k
k
k+1

ak







f yw+1
 f (xk ) = r (xk ) + min
a
k

10

, where k = Kw ,
, if Cw > 1 and ∀ k ∈ {Kw , . . . , Kw + Cw − 2},
, where k = Kw + Cw − 1,

which can be written as





f yw =
E
f (xKw ) ,



ζKw ,...,ζKw +Cw −1



K +C −1




 wXw

∀ w ∈ {1, . . . , W}, 










 + f yw+1 

(x
)
f
x
=
min
r
,


K

k
w



aKw ,...,aKw +Cw −1 




(1)

k=Kw


and f yW+1 = 0. Therefore, the SCRP is concerned with finding f y1 , where by induction:


 



C




X









 

. . .



 min
(x
)
f y1 =
E
E
E
min
r
.
.
.
.








k









a
,...,a
a
,...,a




ζK ,...,ζK +C −1 K1
ζK ,...,ζK +C −1 KW
K1 +C1 −1 ζK ,...,ζK +C −1
KW +CW −1


1

1

2

1

2

2

W

W

W

(2)

k=1

The online model Using our notations, we briefly present the SCRP under the online model. Instead
of Assumption A∗6 , the online model assumes that for each retrieval k ∈ {1, . . . , C}, only the next target


defined such that k ∈
container is revealed (i.e. ζk ). Therefore, we consider the states yok , xok
ζk

k=1,...,C

ak

is the empty configuration. In this case, if f o denotes the
{1, . . . , C},
−→
−→
where
 
h  i


cost-to-go function, then by definition, we have f o yok = E f o xok (with f o yoC+1 = 0), and ∀ k ∈
ζk
n  

o
 
o
o
o
o o
{1, . . . , C}, f xk = min r xk + f yk+1 . By induction, the SCRP under the online model is hence
yok

xok

yok+1 ,

yoC+1

ak

concerned with finding:

 
 


C



X
 
 








 

o
o






. . .  .
r xk 
E . . . E min 
f y1 = E min 

a1 


 
ζ2  ζC  aC 
ζ1
o

k=1

The next lemma compares the batch and the online models theoretically (the proof can be found in
the Appendix). It states that it is beneficial in terms of the expected number of relocations to use the
batch model compared to the online model, if the first one is applicable. Practically, this suggests that
the operator should always use available information.
Lemma 1. Let y be a given initial configuration, then we have


f y 6 fo y .

1.3

Decision trees

Multi-stage stochastic optimization problems can be solved using decision trees in which chance
nodes and decision nodes typically alternate. A chance node embodies the stochasticity of the model, while
a decision node models the possible actions of the algorithm. In a decision tree for the SCRP, a node
represents a configuration. The root node (denoted by 0) is the initial configuration, and the leaf nodes
are configurations for which we can compute the cost-to-go function.
In our case, we slightly modify the structure of a typical decision tree, in the sense that chance nodes
and decision nodes do not necessarily alternate. A chance node is a configuration for which the target
container is not known yet, and information needs to be revealed (note that this only occurs at the
beginning of each batch). A decision node is a configuration for which the target container is known.

Using our notations, chance nodes correspond to yw w=1,...,W and decision nodes correspond to (xk )k=1,...,C .
Let n be a node corresponding to a configuration in the decision tree. Thus f (n), the cost-to-go
function of n, is defined for all nodes n, and r(n), the immediate cost function of n, is well defined when
11

n is a decision node. We denote by λn the level of n in the decision tree, and define it as the number
of containers remaining in the configuration. We denote the lowest level of the tree by λ∗ = min {λn }.
n∈Tree

It corresponds to the level such that if λn = λ∗ , f (n) can be computed in an efficient way (the empty
configuration being an obvious candidate with a cost-to-go of 0). Moreover,
• If n is a chance node, then there exists a unique w ∈ {1, . . . , W} such that λn = C − Kw + 1. We denote
by Ωn the set of offspring of n, each offspring being a decision node corresponding to a realization
of the random variables ζKw , . . . , ζKw +Cw −1 , i.e., the full retrieval order of containers in batch w. Note
that n has a priori |Ωn | = Cw ! offspring.
• If n is a decision node, then r(n) is well defined and is equal to the number of containers blocking
the target container in configuration n (i.e., the (C − λn + 1)th container to be retrieved). We denote
by ∆n the set of offspring of n, which can either be chance nodes if there exists w ∈ {1, . . . , W} such
that λn + 1 = C − Kw + 1, or decision nodes otherwise. For the sake of simplicity, we compute ∆n
greedily by considering all feasible combinations of relocations of the r(n) containers blocking the
target container in n. Note that |∆n | is of the order of r(n)S−1 , where S is the number of stacks.
Equation (1) provides the relation to compute the cost-to-go by back-tracking. For all n in the decision
tree, we have

1 X



f (ni ) , if n is a chance node,


 |Ωn |
ni ∈Ωn
(3)
f (n) = 





f
(n
)
,
if
n
is
a
decision
node.
r(n)
+
min
i

n ∈∆
i

n

In the case in which the probability of each permutation (in each batch) is not uniform, we mentioned
that in practice, operators provide the probability of potential orders for each batch. Given a chance node
n and one of its offspring ni ∈ Ωn , this input probability is exactly the probability that the actual retrieval

order is the one revealed in ni . For a given node n, we denote these probabilities by pni ni ∈Ωn . In this
case, Equation (3) is replaced by:
 X

pni f (ni ) , if n is a chance node,



 ni ∈Ωn
f (n) = 




 r(n) + min f (ni ) , if n is a decision node,
ni ∈∆n

for all n in the decision tree.
Figures 5 and 6 provide the description of the decision tree corresponding to the example in Figure 4,
using chance/decision nodes and configurations, respectively. Chance nodes are depicted with a circle,
and decision nodes with a square.
To illustrate how to use Equation (3), we derive the calculations using the example in Figure 5. Suppose
that f (25) = f (26) = f (27) = 0.5 are known, then we get f (17) = f (18) = f (19) = 0.5, f (07) = f (09) = 1.5
and f (08) = 2.5, leading to f (01) = 3.5. Similarly, by back-tracking, we can compute f (02) = f (03) =
f (04) = 1.5 and f (05) = f (06) = 2.5 giving us f (00) = 13/6.
As the example shows, considering the full decision tree can become intractable even for the small
examples, hence quickly impossible for larger problems. As previously mentioned, the number of
decision offspring of a chance node scales with Cw !, and the number of offspring of a decision node is of
12

00

07

01

02

03

04

05

06

08

09

10

11

12

13

14

15

17

18

19

20

21

22

23

24

25

26

16

27

Figure 5: Decision tree represented with nodes. The colored arrows represent different values of immediate cost, i.e., the number of containers blocking the target container (dotted yellow: 0, dashed orange:
1, solid red: 2)

the order of r(n)S−1 . So the size of the tree grows exponentially with the size of the problem. However,
there exist general and specific techniques to reduce substantially the size of this tree, as we discuss now:
First, there are suboptimal approaches. One way to approximate f (n), when n is a chance node, is
to sample from its offspring. When Ωn is large (which can happen in the case of large batches), one
might sample a certain number of realizations, resulting in a set of offspring Ψn ⊂ Ωn . By sampling
“enough,” we show in Section 4, that we can ensure guarantees on expectation for such an algorithm.
Another popular suboptimal approach is to use techniques from Approximate Dynamic Programming.
These techniques can provide good empirical results, but no guarantee on how far from optimality can
be obtained. This direction is not discussed in this paper but can be a future direction of work. Finally,
another approach is to consider heuristics such as the ones described in the next section, which select a
subset of the offspring of decision nodes, and lead to upper bounds on the optimal value f (0).
Second, there exist ways to decrease the size of the tree, while ensuring optimality. One of them
is to reduce the number of nodes using the problem structure of the SCRP. In the online setting, Ku
and Arthanari (2016) propose an “abstraction” technique, which shrinks significantly the size of the tree.
Thanks to Assumption A∗4 , we can consider that the stacks of a configuration are interchangeable, making
many configurations equivalent in terms of number of relocations. For instance, in Figure 5, nodes 20
and 21 are identical in terms of number of relocations.
We describe the “abstraction” step with an example in Figure 7. The five configurations at the top are
all equivalent to the configuration at the bottom. The abstraction transformation first ranks the columns
by increasing height. For stacks with equal height, it breaks ties by ranking them lexicographically
starting from bottom to top. Columns are re-arranged in order to have the first ranked on the left, and

13

1
1

3
2

5

4

5

1

2
1

5

4

5

3

3

4

4

5

3

5

3

5

5

4

2

5

2

5

2

5

5

3

4

5

1

1
2

5

4

5

3

5

4

5

3

5
5

1

3

4

5

2

4

4

5

4

5

4

4

5

2

5

2

5

3

5

2

5

2

3

5

3

5

3

5

4

4

5

3

3

5

3

1

3
4

4

5

2

5

3

5

5

4

5

3

3
2

4

5

5

2

5

2

4

5

5

5

4

5

1

2

4
3

3

5
3

5

5
3

5

5
4

4

5

4
5
5

5
4

5

5
4

5

Figure 6: Decision tree represented with configurations. The colored arrows represent different values
of immediate cost, i.e., the number of containers blocking the target container (dotted yellow: 0, dashed
orange: 1, solid red: 2). Red numbers highlight containers blocking the target container.

14

5

4

4

5

4

5

5

4

5

4

5

3

3

5

3

5

5

3

5

3

4

5

3

5

Figure 7: “Abstraction” procedure

the last on the right. Ku and Arthanari (2016) use a slightly different rule, and add the extra-step of
removing empty columns. Using the abstraction procedure, the proposed algorithm should not generate
twice nodes with identical abstracted versions. Even though Ku and Arthanari (2016) introduce this rule
for the online model, this “Abstraction” step also applies in the batch setting. Throughout the rest of the
paper, we will refer in pseudocode to the function Abstract(n), when applying this method to a given
node n. Finally, we mention that Ku and Arthanari (2016) also suggest caching strategies that could be
added on the top of this simplification step, including caching part of the tree, or using a transportation
table.
Finally, the performance of a decision tree based algorithm depends on the exploration strategy of
the tree. For the online model, Ku and Arthanari (2016) use depth-first-search (DFS), and we propose to
explore the best-first-search (BFS) approach. Note that BFS requires some kind of measure that we define
in Section 3.
In further sections, we explore two other ways to decrease the size of the tree while still ensuring
optimality. The first one is specific to the SCRP, and uses properties of the problem to increase λ∗ . Recall
that λ∗ is the minimum level of the tree at which we can find the optimal expected number of relocations,
without further branching. We show that we can set λ∗ to max{S, CW }, where S is the number of stacks,
and CW is the number of containers in the last batch. Comparatively, Ku and Arthanari (2016) branch
until λ∗ = 0. The second optimal pruning strategy uses lower bounds in a BFS scheme.
After introducing the batch model for the SCRP (as well as the online model) and some preliminary
concepts about decision trees, the next three sections develop approaches to solve the SCRP.

2

Heuristics and lower bounds

Before introducing the two main algorithms, we decribe in this section heuristics and lower bounds
for the SCRP. Indeed, PBFS and PBFSA build upon some of these bounds. In addition, these algorithms
provide good intuition on how to solve the SCRP.
Let n be a given configuration with S stacks and T tiers. We say that a Container c is a blocking container
in n if it is stacked above at least one container which is to be removed before c. Note that all the bounds
15

mentioned below apply both in the batch and online models.

2.1

Heuristics

For the sake of completeness of this paper, we first describe three existing heuristics used in proofs
and/or our computational experiments before describing two novel heuristics.
2.1.1

Existing heuristics

Random For every relocation of a blocking Container c from Stack s, the random heuristic picks any
Stack s0 , s uniformly at random among stacks that are not “full,” i.e., stacks containing strictly less than
T containers.
Leveling (L) For every relocation of a blocking Container c from Stack s, L chooses the Stack s0 , s
currently containing the least number of containers, breaking ties arbitrarily by selecting such leftmost
stack.
Heuristic L is interesting for several reasons. Most importantly, it is an intuitive and commonly used
heuristics in real operations as it uses no more than the height of each stack. It does not require any
information about batches or departure times, which means it is robust with respect to the inaccuracy
of information. In addition, it is optimal for any configurations with S containers or less (see Section
3). Finally, we show strong evidence in the last computational experiment (see Section 5), that this
policy is optimal for the SCRP under the online model with a unique batch (representing the case of
no-information on the retrieval order).
Expected Reshuffling Index (ERI) This index-based heuristic was introduced by Ku and Arthanari
(2016) for the online model. For every relocation of a blocking Container c from Stack s, ERI computes
a score called the expected reshuffling index for each Stack s0 , s that is not full, denoted by ERI(s0 , c).
ERI chooses the Stack s0 , s with the lowest ERI(s0 , c). In the case of a tie, the policy breaks it by
selecting the highest column among the ones minimizing ERI(s0 , c). Further ties are arbitrarily broken
by selecting the leftmost column verifying the two previous conditions. ERI(s0 , c) corresponds to the
expected number of containers in Stack s0 that depart earlier than c. Let Hs0 be the current number of
containers in s0 . If Hs0 = 0, then ERI(s0 , c) = 0. Otherwise, let (c1 , . . . , cHs0 ) be the containers in s0 , then
Hs0
X
1 {ci = c}
.
ERI(s0 , c) =
1 {ci < c} +
2
i=1

2.1.2

First new heuristic: Expected Minmax (EM)

EM considers an idea similar to that of Caserta et al. (2012) for the CRP. Let min(s) be the smallest
label of a container in s (min(s) = C + 1, if Stack s is empty). For every relocation of a blocking Container
c from Stack s, we select the stack to which we relocate c using the following rules:
[Rule 1] If there exists s0 , s such that min(s0 ) > c, let M =

min {min(s0 ) : min(s0 ) > c}. Select a stack such

s0 ∈{1,...,S}\s

that min(s0 ) = M, breaking ties by choosing from the highest ones, finally taking the leftmost stack
if any ties remain.
16

[Rule 2] If for all Stacks s0 , s, min(s0 ) 6 c, let M =

max {min(s0 )}. Select a stack such that min(s0 ) = M.

s0 ∈{1,...,S}\s

If there are several such stacks, select those with the minimum number of containers labeled M.
Further ties are again broken by taking the highest ones, and finally choosing the leftmost one
arbitrarily.
2
4

4

7

7

4

4

7

1

7

3

4 8 3
Rule 1

→
−

2

4

4

7

1

7

3

4 8 2
Rule 1

7

→
−

2

4

4

1

7

7
4

3

2

→
−

4

4

7

7

1

7

4

3

4 4 2
Rule 2

4

→
−
1

2

4

7

7

7

4

3

4 4 2
Rule 2

Figure 8: Decisions of the EM heuristic on an example with 5 tiers, 4 stacks and 9 containers (3 per batch).
Under the batch model, the first batch has been revealed and we present the decisions to retrieve the
first container made by EM. The container with the red label is the current blocking container. Numbers
under the configuration correspond to the stack indices min(s). The green (respectively orange) indices
correspond to the selected stack with the corresponding M when Rule 1 (respectively Rule 2) applies.
Rule 1 says: if there is a stack where min(s) is greater than c (c can almost surely avoid being relocated
again), then choose such a stack where min(s) is minimized, since stacks with larger minimums can be
useful for larger blocking containers.
If there is no stack satisfying min(s) > c (Rule 2), then we have two cases following the same rule. On
one hand, if M = c, then M is the maximum of the minimum labels of each stack, and c can potentially
avoid being relocated again. If there are several stacks that maximize the minimum label, then by
selecting the one with the least number of containers labeled M, EM minimizes the probability of c being
relocated again. On the other hand, if M < c, c will almost surely be relocated again, then EM chooses
the stack where min(s) is maximized in order to delay the next unavoidable relocation of c as much as
possible. We show how EM makes decision on a simple example in Figure 8.
2.1.3

Second new heuristic: Expected Group assignment (EG)

EM is quite intuitive because it tries to minimize the number of blocking containers after each
retrieval. EG aims for the same goal, but uses some more sophisticated rules (althought, as shown in
the experiments in Section 5, EG does not always provide better solutions that EM). EG is inspired by a
heuristic designed by Wu and Ting (2012) for the complete information case, and we generalize this idea
to the SCRP. It is different from ERI and EM because it considers a group of blocking containers together,
while ERI and EM consider them one at a time. EG can be decomposed in two main phases for each
retrieval. The decisions made by EG on the same example as Figure 8 are given in Figure 9.
The first phase assigns the blocking containers for which there exists s0 , s such that min(s0 ) > c. If
this is not the case, the assignments of these containers will be ignored at the first phase. The acceptable
containers are assigned in descending order of labels, i.e., the container with highest label is assigned
first (breaking ties for the highest one first). In order to assign these acceptable containers, the first phase
17

2

2

2

2

4

4

4

4

4

7

2

4

4

7

1

7

3

4

8

3

→
−d

4

4

7

1

7

7

3

4

×

2

→
−

4

4

7

1

7

7

3

4

7

2

→
−
1

4

4

7

7

7

3

4

×

3

→
−
1

4

4

7

7

7

3

Non-assigned

9a. First phase: EG assigns acceptable containers in descending order. The container with the red label is
the next acceptable container that EG tries to assign to a stack. Containers with blue labels are assigned,
gray are unassigned. Below, we show the indices min(s) to apply the first rule of EM (× means that a
container below the considered container has already been assigned to a stack).

4
2
1

4

4

7

7

7

3

4

0

2

4

→
−
1

2

4

4

7

7

7

3

9b. Second Phase, EG assigns all unassigned containers using the index Gmin(s).
Figure 9: Decisions of the EG heuristic on an example with 5 tiers, 4 stacks and 9 containers (3 in each
batch). Under the batch model, the first batch has been revealed and we present the two phase decisions
to retrieve the first container made by EG.
applies the first of the EM rules. Finally, acceptable containers cannot be assigned to a stack if there is a
container below it that was previously assigned to this stack.
The assignment in the second phase for the blocking containers not assigned yet, might lead to
additional relocations. These containers are assigned to other stacks in ascending order of labels. The
second phase first computes a modified min(s0 ) index for each stack denoted by Gmin(s0 ), which is defined
as follows: Let Hs0 be the height of Stack s0 and B(s0 ) be the subset of containers assigned in the first phase
to Stack s0 ,



−1, if |B(s0 )| + Hs0 = T,




0
0

 min(s ), if |B(s )| = 0,
Gmin(s0 ) = 


 B(s0 ), if |B(s0 )| = 1,



 0, otherwise.
If a stack is full after we assign the containers during the first phase, then it cannot be selected. If no
container was assigned, the index remains as the min. If one container was assigned, it is “artificially”
the new minimum of the stack. Finally, if more than one container was assigned, the index becomes
very unattractive by being as low as possible (0). The second phase is similar to the EM heuristic, but
it considers Gmin instead of the min index, breaking ties identically. Note that, after each assignment in
the second phase, we update Gmin accordingly for the remaining containers to be assigned. For more
18

details in the complete information case, we refer the reader to Wu and Ting (2012).

2.2

Lower bounds

After defining heuristics (upper bounds), we are now concerned with defining valid lower bounds
for the SCRP. More specifically, we care about computing lower bounds for decision nodes in the decision
tree defined before. Note that the computation of lower bounds easily extends to chance nodes.
2.2.1

Blocking lower bound

Suppose that the departure order is known, like in the CRP. The following lower bound was introduced by Kim and Hong (2006) and it is based on the following simple observation. If a container is
blocking in n, then it must be relocated at least once. Thus, the optimal number of relocations is lower
bounded by the number of blocking containers.
In the SCRP, the retrieval order is a random variable, so the fact that a container is blocking is also
random. Let us denote the expected number of blocking containers in n by b(n). Therefore, by taking the
expectation on the retrieval order of the previous fact, which holds for every retrieval order, we have the
following observation.
Observation 1. For all configurations n, f (n) is the minimum expected number of relocations to empty n, and
b(n) is the expected number of blocking containers, then
f (n) > b(n).
Lemma 2 shows one way to compute the expected number of blocking containers for one stack, and
b(n) is the sum of the expected number of blocking containers of each stack of n. Mathematically, let bs (n)
S
X
be the expected number of blocking containers in Stack s of n, we have b(n) =
bs (n).
s=1

Lemma 2. Let n be a single stack configuration with T tiers, and H > 0 containers (H 6 T). If H = 0, we have
b(n) = 0.
If H > 1, we denote the label of containers by (ci )i=1,...,H , where c1 is the container at the bottom and cH at the top
(see Figure 10), then we have:
(
)
H 1 ch = min {ci }
X
i=1,...,h
,
b(n) = H −
h
X
h=1
1 {ch = ci }
i=1

where 1 {A} is the indicator function of A.
Proof. Clearly, if H = 0, then b(n) = 0. If H > 1, then by definition, we have
 H

H
X 
 X



b(n) = E 
1 ch is a blocking container  =
P ch is a blocking container .
h=1

h=1

Let us fix h ∈ {1, . . . , H}, and compute the probability that ch is blocking. We consider two cases:
19
















T















cH
...
c2
c1

Figure 10: Example of a single stack configuration
• If ch > min {ci }, then ch is almost surely blocking.
i=1,...,h

• Otherwise ch = min {ci }, and there are
i=1,...,h

h
X

1 {ch = ci } − 1 containers below ch with the same label (or

i=1

batch). Since each departure sequence between containers of the same batch is equally likely, the
Ph

probability that ch is blocking is equal to

1{ch =ci }−1
=1−
i=1 1{ch =ci }

i=1

Ph

1

i=1 1{ch =ci }

Ph

.

Consequently, we get

(
) 
(
)




1
 × 1 ch = min {ci }
P ch is a blocking container = 1 × 1 ch > min {ci } + 1 − Ph

i=1,...,h
i=1,...,h
1
=
c
}
{c
i
h
i=1
)
(
1 ch = min {ci }
i=1,...,h
= 1 − Ph
.
i=1 1 {ch = ci }
We sum the above expression for h = 1, . . . , H to conclude the proof.



Therefore, one can compute the blocking lower bound as follows: let Hs be the number of containers


in Stack s, and cs1 , . . . , csHs be the containers in Stack s listed from bottom to top, then


(
)

n
o



s
s

s 1 c = min c

H
h
X 
X

i=1,...,h i
 s
 .
b(n) =
H −

h

X

n
o
s=1,...,S 
h=1


s
s

1
c
=
c

Hs >1 
i
h


(4)

i=1

Non-uniform case In the case where probabilities are not uniform across retrieval orders, we still
consider a similar lower bound. For each Container csh , let qcsh be the probability that csh is the first
container to be retrieved among the ones with the with same batch, and positioned below in its stack.
Equation (4) extends to give:

(
)
Hs
X 
X
n o 
Hs −
b(n) =
qcsh 1 csh = min csi  .

i=1,...,h
s=1,...,S
Hs >1

h=1

20

2.2.2

Look-ahead lower bounds

Note that the blocking lower bound b is only taking into account the current configuration. However,
some relocations lead necessarily to an additional relocation. We refer to such relocations as “bad.”
Formally, let s be a stack of a configuration, and min(s) be the smallest label of a container in s. Recall
that, if s is empty, we set min(s) = C + 1. We say that the relocation of Container c from Stack s is “bad” if

c > 0 max 0 min(s0 ) . We propose to construct a lower bound that anticipates “bad” relocations.
s =1,...,S , s ,s

The basic idea is based on a similar one used by Zhu et al. (2012) for the CRP. We consider the 1st
look-ahead lower bound denoted by b1 (n). By definition, we take b1 (n) = b(n) + d1 (n), where b(n) is
the blocking lower bound, and d1 (n) is the expected number of unavoidable “bad” relocations while
performing the first removal. We compute the term d1 (n) by considering all realizations of the first target
container. For each realization, we compute the number of unavoidable “bad” relocations, and average
them. Formally, for a given configuration n, consider Un(the set of potential next
target container in n
)
(which can be a singleton if it is known already), i.e., Un = c c = min (min(s)) . Based on the definition
s=1,...,S

of a bad relocation, we compute the number of unavoidable “bad” moves for each u ∈ Un denoted by
β(n, u), and we take:
1 X
β(n, u),
d1 (n) =
|Un |
u∈Un

or d1 (n) = u∈Un pn,u β(n, u), where pn,u is the probability that u is the next target container in n if the

probabilities considered are not uniform (which can be computed using pni ni ∈Ωn if n is a chance node).
P

4
3
1

3

1

Figure 11: Example for look-ahead lower bounds
For example, in Figure 11, the presented configuration denoted by n is such that b(n) = 2. Now
consider a container u ∈ Un : if u is the container labeled 1 in Stack 1, then there is no blocking container,
so β(n, u) = 0; if u is the other container labeled 1, the relocation of the container labeled 4 from Stack 3 is
necessarily a bad relocation, since min(1) = 1 < 4 and min(2) = 3 < 4, but it is not the case for the blocking
container labeled by 3, hence β(n, u) = 1. Therefore, d1 (1) = 0.5(0 + 1) = 0.5, and b1 (n) = 2 + 0.5 = 2.5,
hence giving a lower bound closer to the optimal solution than b(n). Note that, if n has an empty stack,
then β(n, u) = 0 for all u ∈ Un , and hence d1 (n) = 0.
We can refine this idea, by trying to find unavoidable “bad” relocations for the second removal. In this
case, the configuration depends on the first removal, and the decisions that have been made accordingly.
For the sake of clarity, consider that the first target container has been revealed, and denote it u1 .
After retrieving u1 , only containers blocking u1 have changed from their initial position. It can be very
challenging to detect future unavoidable “bad” moves for these containers. In order to bypass this issue,
we consider that all containers blocking u1 are also removed, resulting in a configuration without u1 and
its blocking containers. Given this new configuration denoted by n(u1 ), we can compute the expected
21

number of unavoidable bad moves d1 (n(u1 )). Since u1 is actually random, we have to consider each
scenario with their associated probability, and compute a new configuration where blocking containers
are retrieved with the target container. We denote the result d2 (n), and it is a lower bound on the expected
number of unavoidable bad relocations for the first two removals starting at n. Finally, our 2nd look-ahead
lower bound is given by b2 (n) = b(n) + d2 (n).
Algorithm 1 Lower bound on the number of unavoidable bad relocations for the k first removals
1:
2:
3:
4:
5:
6:
7:
8:
9:

procedure [dk (n)] = UnavoidableBadReloc (n, k)
if k = 0 or n has an empty stack or n is empty then dk (n) = 0
else let Un = {containers with minimum label in n}
P
if k = 1 then dk (n) = |U1n | u∈Un β(n, u)
else
for u ∈ Un do
Let n(u) be the configuration n without u and all containers blocking u
Compute recursively dk−1 (n(u)) = UnavoidableBadReloc (n(u), k − 1)
P
Compute dk (n) = |U1n | u∈Un β(n, u) + dk−1 (n(u))

This idea can easily be generalized for k > 2 by induction with bk (n) = b(n) + dk (n). Here k is the
number of removals that the lower bound considers to compute the expected number of unavoidable
bad relocations (see pseudocode of Algorithm 1). We mention that we only use the 1st and 2nd look ahead
lower bounds in our computational experiments. However, note that, as k grows, the computational
complexity clearly increases, whereas experiments reveal that the marginal increase of the lower bound,
i.e., bk+1 (n) − bk (n) > 0, decreases.

3

PBFS, a new optimal algorithm for the SCRP

Building upon lower bounds introduced in the previous section, this section introduces, studies and
proves the optimality of one of the main contributions of this paper, the PBFS Algorithm.

3.1

PBFS algorithm

We start by giving the pseudocode of our algorithm, and we derive its optimality in Lemmas 3 and 4.
PBFS takes two inputs, the configuration n for which we aim to compute f (n), and a valid lower bound
l. This algorithm uses a combination of four features to return f (n). The first one is the BFS exploration
of the tree based on a given lower bound l. We first compute f for the “most promising nodes,” because
nodes with small lower bounds are more likely to result in small f . The second technique is stopping to
compute f recursively after level λ∗ = max{S, CW }, by calculating it either using b, or the A∗ algorithm
defined later. The third one is pruning with a lower bound revealing the sub-optimality of some nodes
without actually computing f . Finally, it also uses the abstraction technique described previously.

22

Algorithm 2 PBFS Algorithm
1:
2:
3:
4:
5:
6:
7:
8:
9:

procedure [ f (n)] = PBFS (n, l)
if λn 6 S (n has less than S containers) then f (n) = b(n)
else
if n is a chance node then start with ΨnPBFS = {}
for ni ∈ Ωn do ni ← Abstract(ni )
if there exists m = ni already in ΨnPBFS then pnm ← pnm + |Ω1n |
else if there exists m = ni already in the decision tree then add m to ΨnPBFS and pnm =
else add ni to ΨnPBFS , pnni = |Ω1n | and compute f (ni ) = PBFS (ni , l)
X
pnni f (ni )
Compute f (n) =

1
|Ωn |

ni ∈ΨnPBFS

10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:

3.1.1

else n is a decision node
if λn = CW (the full retrieval order is known) then f (n) = A∗ (n)
else construct ∆n by considering all feasible sets of decisions to deliver the target container
Compute l(ni ) for each ni ∈ ∆n


Sort n(1) , n(2) , . . . , n(|∆n |) in non-decreasing order of l(.)


Compute f (n(1) ) = PBFS n(1) , l
Start with ΓPBFS
= {n(1) } and k = 2
n
n
o
while k 6 |∆n | and l(n(k) ) < min f (n(j) ) do n(k) ← Abstract(n(k) )
j=1,...,k−1

if there exists m = n(k) already in the decision tree then add m to ΓPBFS
n


PBFS
else add n(k) to Γn
and compute f (n(k) ) = PBFS n(k) , l
Update k = k + 1

f (n) = r(n) + min f (ni )
ni ∈ΓPBFS
n

Decreasing the size of decision tree by increasing λ∗ to max{S, CW }

If CW 6 S, then compute f(n) using b(n) Recall that, for every relocation, heuristic L chooses the stack
with the least number of containers, breaking ties arbitrarily by choosing the leftmost one. Note that L
always provides a valid upper bound for the SCRP. So if we denote the resulting expected number of
relocations to empty configuration n using L by fL (n), then we have fL (n) > f (n).
Lemma 3. Let n be a configuration with S stacks, T tiers, and C containers such that C 6 S, then we have
fL (n) = b(n) = f (n)
Proof. Consider a retrieval order of containers from n that has a non-zero probability of occurring. If
there are no blocking containers, then the lemma clearly holds. Otherwise, let c be one of the blocking
containers for this retrieval order, and consider the first removal for which c has to be relocated. For this
removal, there are at most S containers in the configuration, hence there exists at least one empty stack
to relocate c. Since heuristic L chooses always empty columns if one exists, L would move c to one of the
existing empty stacks. Note that in this case, c will never be blocking again, and hence never be relocated

23

again. This observation holds for any blocking containers, thus L relocates each blocking container at
most once.
Since this fact holds for any retrieval orders with non-zero probability, by taking expectation on the
retrieval order, we have fL (n) 6 b(n), thus fL (n) 6 b(n) 6 f (n) 6 fL (n), which concludes the proof.

Lemma 3 states that for configurations with S containers or less, the L heuristic is optimal for the
SCRP. This implies that, for nodes at level S, we have access to the cost-to-go function using b(n), as well
as an optimal solution (provided by heuristic L). Hence PBFS can stop branching at λ∗ = S (line 2 of
Algorithm 2).
If CW > S, then compute f(n) using the A∗ algorithm If n is a decision node at level CW , the full order
of retrieval is known, and computing f (n) reduces to solving a classical CRP, so we can leverage the existence of efficient solutions to the classical CRP such as the A∗ algorithm, and take λ∗ = CW . Throughout
the rest of the paper, A∗ refers to the improved version of this algorithm presented in Borjian et al. (2015a)
and we denote the optimal number of relocations obtained by A∗ (n) (line 11 of Algorithm 2).
Combining with the two previous observations, we can take λ∗ = max{S, CW }.
3.1.2

Decreasing the size of decision tree by pruning using lower bounds

We would also like to reduce the size of the tree before level λ∗ . For a decision node n, PBFS considers
only a subset ΓPBFS
of all the offspring ∆n (line 21 of Algorithm 2). Our goal is to set ΓPBFS
in order to still
n
n
guarantee optimality.

n

...

n(1)




6 ... 6





...

l n(1)

f n(1)

n(|ΓPBFS
|)
n




l n(|ΓPBFS
|)
n


f n(|ΓPBFS
|)
n

n(|ΓPBFS
...
|+1)
n


l n(|ΓPBFS
|+1) 6 . . . 6
n
×

...

n(|∆n |)


l n(|∆n |)
×

Figure 12: Illustration of the pruning rule. First, offspring are ordered by non-decreasing lower bounds.
Then we start computing the objective function starting at n(1) . We stop computing the objective functions
once the pruning rule is reached. In the figure above, green nodes are nodes in ΓPBFS
i.e. f (.) has been
n
PBFS
computed. Orange nodes are nodes in ∆n \ Γn
i.e. f (.) does not need to be computed which is
represented here by ×.

First, PBFS generates all nodes ni ∈ ∆n by considering all feasible sets of decisions to deliver the target
container in n (line 12 of Algorithm 2), and for each of them, compute a lower bound l(ni ), where l is the

24



input lower bound (line 13 of Algorithm 2). Let n(1) , n(2) , . . . , n(|Ωn |) be the list of offspring of n sorted
by non-decreasing lower bound (line 14 of Algorithm 2). The algorithm considers first n(1) , adds it to
and computes f (n(1) ) recursively (line 15-16 of Algorithm 2). Then for k > 2, we consider n(k) ’s
ΓPBFS
n
n
o
and
sequentially, and check if l(n(k) ) < min f (n(j) ) (line 17 of Algorithm 2). If so, add n(k) to ΓPBFS
n
j=1,...,k−1

compute f (n(k) ) recursively. If not, we stop branching on all nodes n(k) , . . . , n(|Ωn |) . An illustration of the
pruning rule is shown in Figure 12 and the next lemma shows the optimality of this rule.
be the subset of nodes considered for this node
Lemma 4. Let n be a decision node in the decision tree, and ΓPBFS
n
in Algorithm 2, and constructed as aforementioned, then we have
min

mi ∈ΓPBFS
n




f (mi ) = min f (ni ) .
ni ∈∆n



Proof. Let n(1) , n(2) , . . . , n(|Ωn |) be the list of offspring of n, sorted by non-decreasing lower bounds. We
consider two cases.
• If ΓPBFS
= ∆n , the statement clearly holds.
n
• Otherwise, there exists k 6 |∆n | such that l(n(k) ) >
∀k0

min

j=1,...,k−1

Note that we have
> k, f (n(k0 ) ) > l(n(k0 ) ) > l(n(k) ) >
n
o

min f (n(j) ) = min f (mi ) .
j=1,...,k−1

n

o
n
o
f (n(j) ) , and ΓPBFS
= n(1) , . . . , n(k−1) .
n
n
o

min f (n(j) ) . Hence min f (ni ) =

j=1,...,k−1

ni ∈∆n

mi ∈ΓPBFS
n


We claim that increasing λ∗ to max{S, CW } together with pruning in a Best-First-Search scheme,
dramatically help in the efficiency of PBFS while keeping the guarantee of optimality. In the case of small
batches, the PBFS algorithm appears to be efficient (see Section 5). However, this algorithm faces the
issue that |Ωn | = Cw ! if n is a chance node. So if Cw is large, typically Cw > 4, the number of nodes to
consider gets too large. We tackle this issue by considering a near-optimal algorithm in the next section.
As a final remark, batches should be as small as possible if information is only at stake. Indeed,
smaller batches correspond to an efficient information system since more information is known about
the retrieval order. But the size of batches is restricted by two intrinsic constraints:
1. Batches should be at least larger than a certain size. Indeed, a terminal offers time slots for trucks
to register, and these slots cannot be too small (in terms of time), as trucks would most certainly
not arrive during their appointed slot due to traffic or other uncertain factors. Therefore, given the
minimum time of a slot, the terminal will allow at least a certain number of trucks to register for
each slot, i.e., the minimum batch size.
2. Batches cannot be too large in order to have the batch model applicable, since in this model, the
appointment time windows are supposed to be the same as or shorter than the target waiting time.
As the target waiting time is limited, there is a limited number of containers that can be retrieved
in a certain batch.
This leads us to consider an alternative to PBFS (see Section 4) in the case of larger batches.

25

4

PBFSA, near-optimal algorithm with guarantees for large batches

Algorithm 3 PBFSA Algorithm
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:

procedure [ f˜(n)] = PBFSA (n, l, )
if λn 6 S then f˜(n) = b(n)
else
if n is a chance node then start with ΨnPBFSA = {}. Let wmin be such that λn = C − Kwmin + 1
n
o
P

∗
Compute δn = min w ∈ {wmin , . . . , W} w
u=wmin Cu > λn − λ to get n =
δn

π fmax (n) − fmin (n) 2
Compute fmax (n) and fmin (n) to get Nn (n ) =
22n
if Nn (n ) 6 Cwmin ! then
for i = 1, . . . , Nn (n ) do
Sample a random permutation, get corresponding ni ∈ Ωn and ni ← Abstract(ni )
if there is m = ni already in ΨnPBFSA then pnm ← pnm + Nn1(n )
else if there is m = ni already in the decision tree then add m to ΨnPBFSA , pnm = Nn1(n )
else add ni to ΨnPBFSA , pnni = Nn1(n ) and compute f˜(ni ) = PBFSA (ni , l,  − n )
else
for ni ∈ Ωn do ni ← Abstract(ni )
if there exists m = ni already in ΨnPBFSA then pnm ← pnm + |Ω1n |
else if there exists m = ni already in decision tree then add m to ΨnPBFSA and pnm =
else add ni to ΨnPBFSA , pnni = |Ω1n | and compute f˜(ni ) = PBFSA (ni , l,  − n )
X
Compute f˜(n) =
pnni f˜(ni )

1
|Ωn |

ni ∈ΨnPBFSA

19:
20:
21:
22:
23:
24:
25:
26:
27:
28:
29:
30:

else n is a decision node
if λn 6 CW then f˜(n) = A∗ (n)
else Construct ∆n by considering all feasible sets of decisions to deliver the target container
Compute l(ni ) for each ni ∈ ∆n


Sort n(1) , n(2) , . . . , n(|∆n |) in non-decreasing order of l(.)


Compute f˜(n(1) ) = PBFSA n(1) , l, 
Start with ΓPBFSA
= {n(1) } and k = 2
n
n
o
while k 6 |∆n | and l(n(k) ) < min f˜(n(j) ) do n(k) ← Abstract(n(k) )
j=1,...,k−1

if there exists m = n(k) already in the decision tree then add m to ΓPBFSA
n


˜(n(k) ) = PBFSA n(k) , l, 
else add n(k) to ΓPBFSA
and
compute
f
n
Update k = k + 1
f˜(n) = r(n) + min

ni ∈ΓPBFSA
n

n

f˜(ni )

o

Building upon PBFS introduced in the previous section, this section describes the randomized algorithm PBFSA and shows some theoretical guarantees on expectation. This new algorithm is identical
to PBFS except when computing the value function of a chance node (lines 4 to 17 of Algorithm 3). In

26

order to decrease the number of decision offspring to consider for each chance node, we sample a certain
number of i.i.d. permutations, and only consider the decision nodes associated with these permutations
as illustrated in Figure 13. The facts that the objective function is bounded and that our problem has
a finite number of sampling stages allow us to independently sample nodes in order to approximate
the objective function. Using concentration inequalities, we can chose the number of samples needed to
control the approximation error.

1
1

5

4

5

1

n

n1
f (n1 )

...

nNn (n )

...



f nNn (n )

ni

...

n|Ωn |

×

×

Figure 13: Illustration of the sampling rule. In this figure, the smallest batch is batch 1 therefore wmin = 1,
and there are 6 containers thus λn = 6. These values allow us to compute the number of samples required
Nn (n ). If Nn (n ) in less than the total number of offspring |Ωn | = Cwmin ! = 3!, then we only compute f (.)
for sampled nodes. ΨnPBFSA represents the subset of sampled nodes collored green and for which f (.)
needs to be computed. Note that ΨnPBFSA = Nn (n ). Orange nodes are nodes in Ωn \ ΨnPBFSA i.e. there
were not sampled and f (.) does not need to be computed which is represented here by ×. Finally, the
approximate value of f (n) is the average of the objective values over all sampled nodes.

Formally, let n be a given chance node, recall that we compute f (n) =

1 X
f (ni ), where each
|Ωn |
ni ∈Ωn

ni ∈ Ωn represents one retrieval order (a random permutation) of batch w (if λn = C−Kw +1). Let Ψn ⊂ Ωn
be defined as the resulting subset of the Nn (n ) (∈ N) offspring drawn i.i.d., where Nn is a function of n ,
1 X
itself a function of n and  > 0 (a target error). The goal is to define Nn (n ), such that f˜(n) =
f (m)
|Ψn |
m∈Ψn

is a “good” approximation of f (n), i.e., f˜(n) − f (n) is bounded by  on average.
PBFSA takes three input arguments, the configuration n for which we want to evaluate f , a valid
lower-bound l and an upper bound  > 0 on the total expected “error” ensured by the algorithm. It
outputs f˜(n), which is a randomized approximation of f (n). Because of the samplings performed in
line 9 in Algorithm 3, the output of PBFSA is random. The average error incurred by the algorithm is
h
i
E f˜(n) − f (n) , where the expectation is taken over the aforementioned samplings. Our main result
h
i
(Lemma 5) states that PBFSA ensures E f˜(n) − f (n) 6 , in other words, PBFSA guarantees an average
error of at most .
Lemma 5. Let n be a configuration with λn > 0 containers, l be a valid lower bound function, and  > 0. If
27

f˜(n) = PBFSA(n, l, ), then
h
i
E f˜(n) − f (n) 6 .

4.1

Hoeffding‘s inequality applied to the SCRP

In order to prove this result, we use Hoeffding‘s inequality to compute the number of samples to
ensure probabilistic guarantees. We first state the well-known inequality, and a direct corollary.
Theorem 1 (Hoeffding‘s inequality). Let X ∈ [xmin , xmax ] be a real-valued bounded random variable with mean
N
1 X
value E [X]. Let N ∈ N and (X1 , . . . , XN ) be N i.i.d. samples of X. If X =
Xi , then we have
N
i=1


∀δ > 0 , P X − E [X] > δ 6 exp

−2Nδ2



(xmax − xmin )2

and


−2Nδ2



∀δ > 0 , P X − E [X] < −δ 6 exp

!

(xmax − xmin )2

,

(5)

!
.

(6)

Corollary 1. Let X ∈ [xmin , xmax ] be a real-valued bounded random variable with mean value E [X]. Let N ∈ N
N
1 X
π (xmax − xmin )2
Xi , then ∀ > 0 such that N >
, we have
and (X1 , . . . , XN ) be N i.i.d. samples of X. If X =
N
22
i=1

+  
E X − E [X]
6 ,
2


−

E X − E [X]
6 ,
2

(7)
(8)

where x+ = max{x, 0} (resp. x− = − min{x, 0}) is the positive (resp. negative) part of x.
o
+ 
 n

Proof. For the first result, define ∆ = X − E [X] = X − E [X] 1 X − E [X] > 0 . Note that ∆ is a nonn
o
negative random variable, and ∀δ > 0, {∆ > δ} = X − E [X] > δ . Let F∆ denote the cumulative distribu



−2Nδ2
tion function of ∆, thus, using Equation (5), 1 − F∆ (δ) = P (∆ > δ) = P X − E [X] > δ 6 exp
,
(xmax −xmin )2
which gives
!
Z ∞
Z ∞
−2Nδ2
(1 − F∆ (δ)) dδ 6
E [∆] =
exp
dδ
(xmax − xmin )2
δ=0
δ=0
√
Z
π (xmax − xmin ) 
(xmax − xmin ) ∞
2
=
exp(−u )du =
6 .
√
√
2
u=0
2N
2 2N

−
The proof of the second result is identical to the first one if we consider ∆0 = X − E [X] =
o
n
o

 n
E [X] − X 1 E [X] − X > 0 and notice that ∀δ > 0, {∆0 > δ} = X − E [X] < −δ , hence 1 − F∆0 (δ) 6


−2Nδ2
exp

2 using Equation (6).
(xmax −xmin )

28

Computing fmin and fmax In order to use Hoeffding‘s inequality, we need to define lower ( fmin ) and

upper ( fmax ) bound functions, such that for each chance node n, fmin (n) 6 min f (ni ) and fmax (n) >
ni ∈Ωn

max f (ni ) .
ni ∈Ωn

Lemma 6. Let n be a chance node, if
fmin (n) = min {b(ni )} ,

(9)

)
  

λn
fmax (n) = min ((λn − S) (T − 1)) + (min {S, λn } − 1) , 2
− 1 max {b(ni )} ,
S
ni ∈Ωn

(10)

ni ∈Ωn

and

(

+

then


fmin (n) 6 min f (ni ) and fmax (n) > max f (ni ) .
ni ∈Ωn

ni ∈Ωn


Proof. Since b(ni ) 6 f (ni ), then we have fmin (n) = min {b(ni )} 6 min f (ni ) .
ni ∈Ωn

ni ∈Ωn

By definition, fmax (n) is the minimum of two valid upper bounds. The first one comes from a basic
observation. If there are λn containers remaining to be retrieved in n, consider two cases:
• If λn > S, take the rth retrieval. If S < r 6 λn , then in order to perform this retrieval, there are at
most T − 1 containers blocking the target container so at most T − 1 relocations are needed. When
there are S or less containers remaining, each container (except the lowest one) is at most relocated
once hence we need at most S − 1 relocations. Combining these two facts, the maximum number
of relocations is bounded by (λn − S) (T − 1) + (S − 1) = ((λn − S) (T − 1))+ + (min {S, λn } − 1).
• If λn 6 S, we know that f (ni ) = b(ni ) 6 λn − 1 = ((λn − S) (T − 1))+ + (min {S, λn } − 1).
This shows the validity of the first upper bound.
For the second upper bound, Zehendner et al. (2016) prove that, in the online case with a unique

 l m
batch, the number of relocations performed by the leveling heuristic (L) is at most 2 λSn − 1 B, where B
is the number of blocking containers. Since L is not using any information about batches (only the height
of stacks), this result holds for both batch and online models with any number of batches. Let ni ∈ Ωn ,
using this result and taking expectation over the retrieval order of containers not unveiled in ni yet, we

 l m
have f (ni ) 6 fL (ni ) 6 2 λSn − 1 b(ni ). By taking the maximum over all ni ∈ Ωn , the latter inequality
results in the second upper bound.

Notice that the previous lemma involves computing min {b(ni )} and max {b(ni )}. The following
ni ∈Ωn

ni ∈Ωn

corollary provides an efficient way of computing these values.
Lemma 7. Let n be a chance node, and wmin ∈ {1, . . . , W} be such that λn = C − Kwmin + 1 (i.e., the minimum
batch in n). For each Stack s of n with Hs > 1 containers, let (csh )h=1,...,Hs be the containers in s, where cs1
is the container at the bottom and csHs at the top (see Figure 10, for the case H = Hs ). Finally, consider
n
o
Cswmin = csh = Kwmin , h = 1, . . . , Hs . Then we have


(
)

n
o



s
s

1
c
=
min
c

i
h
X 
X

i=1,...,h
 s
s

 ,
min {b(ni )} =
(11)
H − Cwmin −


h
ni ∈Ωn
X


n
o
s
s=1,...,S 
h=1,...,H

1 csh = csi 

Hs >1 
csh ,Kwmin
i=1

29

and




X 
X
 s
max {b(ni )} =
H −
ni ∈Ωn

s=1,...,S 
h=1,...,Hs

s

H >1 
csh ,Kwmin

1


)
n o 
csh = min csi 

i=1,...,h
 .

h
X
o 
n
s
s
1 ch = ci 

(

(12)

i=1

Proof. Let n be a chance node, and ni ∈ Ωn be one
X of its decision offspring, such that all containers in batch
wmin have been revealed. Recall that b(ni ) =
bs (ni ), where bs (ni ) is the expected number of blocking
s=1,...,S

containers in Stack s. First, for each Stack s such that Hs = 0, bs (ni ) = 0. Hence b(ni ) =

X

bs (ni ).

s=1,...,S
Hs >1

For each Stack s such that Hs > 1, consider the containers in this stack (csi )i=1,...,Hs . Since all containers
labeled Kwmin , i.e., from batch wmin , are known in ni , we can write
X
X
X
n
o
h
i
h
i
P csh is blocking in ni =
1 csh is blocking in ni +
P csh is blocking in ni .
bs (ni ) =
h=1,...,Hs
csh =Kwmin

h=1,...,Hs

h=1,...,Hs
csh ,Kwmin

h
i
Fix h ∈ {1, . . . , Hs } and csh , Kwmin , then the proof of Lemma 2 uses the fact that P csh is blocking in ni =

n o


1csh = min csi 
X
o
n
i=1,...,h
s
n
o
6 Cswmin . Therefore, we can
.
Finally,
it
is
clear
that
0
6
1
c
is
blocking
in
n
1−
Ph
i
s
s
h
i=1

1 ch =ci

h=1,...,Hs
csh =Kwmin

get the corresponding formulae. As a final remark, note that each of these bounds are tight. Indeed,
consider the offspring of n, in which all containers in batch wmin are in the decreasing (resp. increasing)
order of retrieval from top to bottom, then this offspring has no (resp. Cswmin ) blocking container(s).

Non-uniform case Similar to the blocking lower bound, we can extend Lemma 7 to the case where
probabilities are not uniform across retrieval orders. Recall that qcsh denotes the probability that csh is the
first one to be retrieved among the ones positioned below in its stack and with the same label. Then we
have



(
)
X 
X
n
o

 s
qcsh 1 csh = min csi  ,
min {b(ni )} =
H − Cswmin −


ni ∈Ωn
i=1,...,h

s=1,...,S 
h=1,...,Hs

s
s
H >1

and

ch ,Kwmin





(
)

X 
X
n
o

 s
max {b(ni )} =
qcsh 1 csh = min csi  .
H −


ni ∈Ωn
i=1,...,h

s=1,...,S 
h=1,...,Hs

s
s
H >1

ch ,Kwmin

Now we can prove Lemma 5.
Proof of Lemma 5. The proof is by induction on λn . Throughout the proof, we use the same notations as
the ones introduced in Algorithm 3. We say that f˜ verifies Conditions (A) and (B) at node n, if it verifies
respectively the first and second inequalities below:


−  
+  
˜
˜
6 .
E f (n) − f (n)
6 and E f (n) − f (n)
2
2
30

i
h
Note that if f˜ verifies Conditions (A) and (B) at node n, then E f˜(n) − f (n) 6 , which would prove the
lemma. Given  > 0, and l a valid lower bound, the induction hypothesis is:
If f˜(n) = PBFSA(n, l, ), then f˜ verifies Conditions (A) and (B) at node n.
First, if λn 6 S, then f˜(n) = b(n) = f (n), and therefore, f˜ verifies Conditions (A) and (B) at node n. In
this case, f˜(n) is actually deterministic since no sampling is performed by PBFSA.
From now on, consider n such that λn > S.
If n is a decision node such that f˜(n) = PBFSA(n, l, ) and λn > S. First, if S < λn 6 CW , then
f˜(n) = A∗ (n) = f (n), hence f˜ verifies Conditions (A) and (B) at n.
n
o

If λn > max{S, CW }, consider ñ = argmin f˜(ni ) and n∗ = argmin f (ni ) . Note that f˜(n) − f (n) =
ni ∈∆n

ni ∈ΓPBFSA
n

f˜(ñ) − f (n∗ ) almost surely (a.s.), and by definition, ñ and n∗ are both such that λñ = λn∗ = λn − 1 < λn .
Consider the following measurable event:
n
o
(13)
E = f˜(n) − f (n) = f˜(ñ) − f (n∗ ) > 0 .

−
• Conditioned on E, we have f˜(n) − f (n) = 0 a.s., thus


−
˜
E = 0.
E f (n) − f (n)

(14)

Now let us show that conditioned on E, n∗ ∈ ΓPBFSA
a.s.; we suppose by contradiction that a.s.
n
n
o


n∗ < ΓPBFSA
. If k = ΓPBFSA
, then k < |∆n | a.s., and min f˜(n(j) ) 6 l n(k+1) a.s. By definition f˜ (ñ) =
n
n
j=1,...,k
n
o


∗
PBFSA
, then there exists k∗ ∈ {k + 1, . . . , |∆n |} such that
min f˜(n(j) ) so f˜ (ñ) 6 l n(k+1) a.s. Since n < Γn
j=1,...,k
 




n∗ = n(k∗ ) . Since n(i)
are ordered by non-decreasing l(.), we have l n(k+1) 6 l n(k∗ ) = l (n∗ ).
i∈{1,...,|∆n |}
Therefore f˜ (ñ) 6 l (n∗ ) a.s.; but, conditioned on E, f˜ (ñ) > f (n∗ ) > l (n∗ ) a.s., which leads to a
contradiction. Thus conditioned on E, n∗ ∈ ΓPBFSA
a.s.. Therefore, we have f˜(n∗ ) = PBFSA(n∗ , l, ).
n
By induction, f˜ verifies Condition (A) at node n∗ , thus

E

  
∗
∗ +
˜
f (n ) − f (n )
6 .
2

(15)

n
o
Finally, since ñ = argmin f˜(ni ) and n∗ ∈ ΓPBFSA
, then f˜(ñ) 6 f˜(n∗ ) a.s., so we have f˜(ñ) − f (n∗ ) 6
n
ni ∈ΓPBFSA
n

+

+

+

f˜(n∗ ) − f (n∗ ) a.s. Consequently we have f˜(n) − f (n) = f˜(ñ) − f (n∗ ) 6 f˜(n∗ ) − f (n∗ ) a.s.,
resulting in




+

∗
∗ +
˜
˜
E f (n) − f (n)
E 6 E f (n ) − f (n )
E .
(16)

+
• Conditioned on E, we have f˜(n) − f (n) = 0 a.s., thus


+
˜
E f (n) − f (n)
E = 0.

(17)

Moreover, by definition ñ ∈ ΓPBFSA
a.s., and f˜(ñ) = PBFSA(ñ, l, ), thus the induction hypothesis can
n
be applied to ñ. In particular, we have

−  
˜
E f (ñ) − f (ñ)
6 .
(18)
2
31

Finally, it is clear that f (ñ) > f (n∗ ), then f˜(ñ) − f (n∗ ) > f˜(ñ) − f (ñ) a.s., which is equivalent to

− 
− 
−
f˜(n) − f (n) = f˜(ñ) − f (n∗ ) 6 f˜(ñ) − f (ñ) a.s., resulting in




−
−
˜
˜
E 6 E f (ñ) − f (ñ)
E .
(19)
E f (n) − f (n)
Finally, note the following observation: Let Y > 0 a.s., and F be measurable, then we have
h
i  
E [ Y | F ] P (F ) 6 E [Y] and E Y | F P F 6 E [Y] .
Now we can derive






+  
+
+ 
+
E P (E) 6 E f˜(n∗ ) − f (n∗ )
6 ,
E P (E) 6 E f˜(n∗ ) − f (n∗ )
E f˜(n) − f (n)
= E f˜(n) − f (n)
2
where the first equality comes from Equation (17), the first inequality uses Equation (16), the second one

+
holds thanks to f˜(n∗ ) − f (n∗ ) > 0 a.s., and the last one is Equation (15). Therefore, f˜ verifies Condition
(A) at node n.
Similarly, we have
  

  



−
−  
− 
−
E P E 6 E f˜(ñ) − f (ñ)
E P E 6 E f˜(ñ) − f (ñ)
6 ,
E f˜(n) − f (n)
= E f˜(n) − f (n)
2
where the first equality comes from Equation (14), the first inequality uses Equation (19), the second one

−
holds thanks to f˜(ñ) − f (ñ) > 0 a.s., and the last one is Equation (18). Therefore, f˜ verifies Condition
(B) at node n.
Therefore, we have proven that if n is a decision node with λn > S, f˜ verifies both Conditions (A) and
(B) at node n, which proves the lemma for decision nodes.
X
pnni f (ni ),
If n is a chance node such that f˜(n) = PBFSA(n, l, ), and λn > S. Let us define f (n) =
ni ∈ΨnPBFSA

and show that


+   −  n
−   − n
f˜(n) − f (n)
6
and E f˜(n) − f (n)
6
.
(20)
2
2
Recall that ∀ni ∈ ΨnPBFSA , λni = λn , and ni are decision
nodes such
that f˜(ni ) = PBFSA(n
i , l,  − n ).


−
+   −  n
and E f˜(ni ) − f (ni )
6
Therefore, using the previous result, we know that E f˜(ni ) − f (ni )
6
2
 − n
. We derive the following calculations:
2

+ 



 X
 
 X

+ 




+

 


= E 
E f˜(n) − f (n)
pnni f˜(ni ) − f (ni )   6 E 
pnni f˜(ni ) − f (ni ) 

 


ni ∈ΨnPBFSA
ni ∈ΨnPBFSA


X
X

+
 − n  − n
=
.
=
pnni E f˜(ni ) − f (ni )
6
pnni
2
2
PBFSA
PBFSA


E

ni ∈Ψn

ni ∈Ψn

Similarly, we have

E


− 


 X

− 

 


 X
−


f˜(n) − f (n)
= E 
pnni f˜(ni ) − f (ni )   6 E 
pnni f˜(ni ) − f (ni ) 

 


ni ∈ΨnPBFSA
ni ∈ΨnPBFSA


X
X

−
 − n  − n
=
,
=
pnni E f˜(ni ) − f (ni )
6
pnni
2
2
PBFSA
PBFSA
ni ∈Ψn

ni ∈Ψn

32

which proves Equation (20).
If Nn (n ) > Cwmin !, then f (n) = f (n) so f˜(n) − f (n) = f˜(n) − f (n) a.s., and since
implies that f˜ verifies Conditions (A) and (B) at node n.
Otherwise, we have Nn (n ) 6 Cwmin !. Since

ΨnPBFSA

−n
2

is constructed using Nn (n ) =

6 2 , Equation (20)

π fmax (n) − fmin (n)

samples, thus by using Corollary 1, we have


+  n
−   n
and E f (n) − f (n)
E f (n) − f (n)
6
6
2
2

2

22n

(21)

By combining Equations (20) and (21), we have



+ 
+ 
+   − n n

˜
˜
6 E f (n) − f (n)
E f (n) − f (n)
+ E f (n) − f (n)
6
+
= ,
2
2
2



− 
− 
−   −  n  n

E f˜(n) − f (n)
6 E f˜(n) − f (n)
+
= ,
+ E f (n) − f (n)
6
2
2
2
which shows that f˜ verifies Conditions (A) and (B) at node n and concludes the proof.

5



Computational experiments

Having introduced lower and upper bounds, PBFS, PBFSA, and theoretical guarantees in previous
sections, we present several experimental results in this section to understand the effectiveness of our
algorithms for the SCRP. For clarity, we refer to the set of instances from Ku and Arthanari (2016) as the
existing dataset. We present 4 sets of experiments:
1. Based on instances from the existing dataset, which have relatively small batches, we test the PBFS
algorithm, as well as the two new heuristics and our lower bounds.
2. We slightly modify the existing dataset to obtain the modified dataset, in order to obtain instances
with relatively larger batches. We test the efficiency of PBFSA on this modified dataset.
3. Based on the existing dataset, we show that PBFS improves on the algorithm proposed in Ku and
Arthanari (2016) for the online model. Moreover, the two new heuristics (EM and EG) outperform
the ERI algorithm on expectation for the majority of the instances of the dataset.
4. We change the existing dataset by considering that all containers belong to a unique batch. We
show strong computational evidence to support Conjecture 1, which states that the leveling policy
is optimal for the SCRP under the online model with a unique batch.
All experiments are performed on a MacBook Pro with 2.2 GHz Intel Core i7 processor, 8.00 GB of RAM
and the programming language is MATLAB 2016a. Finally, all results and instances used in this section
are available at https://github.com/vgalle/StochasticCRP.
Implementation of heuristics
1. Computing the number of relocations using b when there are S containers or less: In the retrieval
process, when there are S containers or less remaining in the configuration, the expected number
33

of relocations performed by ERI, EM, EG and L is computed using b. This is motivated by
the following observation: ERI, EM, EG and L are optimal when there are S containers or fewer
remaining in the configuration, and Lemma 3 shows that the optimal expected number of relocations
in this case is equal to b. Therefore, for all heuristics (except Random), instead of running simulations
until there are no containers left, we stop when there are S containers left and compute the expected
number of relocations using b instead.
2. Estimate the expected number of relocations using sampling: In order to estimate the exact objective
value for a given heuristic, one would have to consider all possible retrieval scenarios. Instead, for
each heuristic unless specified otherwise, we report the average over 5000 samples (of retrieval
orders) for each instance, where samples are uniformly drawn at random.
Existing dataset description The full description of the dataset can be found in Ku and Arthanari (2016)
and the original data set is available at http://crp-timewindow.blogspot.com. Note that:
• Configuration sizes vary from T = 3, . . . , 6 tiers, and S = 5, . . . , 10 stacks.
• Two occupancy rates are considered, 50 and 67 percent. The occupancy rate (µ ∈ [0, 1]) is defined

such that the initial number of containers is C = round µ × S × T , where round(x) rounds x to the
closest integer. Therefore, a given triplet (T, S, µ) is equivalent to a given triplet (T, S, C), and note
that if C = round (0.67 × S × T), the condition 0 6 C 6 ST − (T − 1) is satisfied.
• Given a configuration size (T and S) and an occupancy rate (µ) resulting in a given initial number
of containers (C), the dataset includes 30 different initial configurations.
• For all 1’440 instances, the ratio between the number of batches and C is taken to be around half,
i.e., there are on average two containers per batch, which is the smallest size for a batch.
In all our experiments, the time limit is set to an hour, and the 1st look-ahead lower bound b1 is
used as input for both PBFS and PBFSA. All instances are solved by heuristics and lower bounds within
seconds or less.

5.1

Experiment 1: Batch model with small batches

"

Table 2 gives a summary of the results as follows: indicates that all 30 instances are solved optimally
by PBFS. In this case, the average solution time to solve these instances is given in seconds. Otherwise,
the number of instances solved optimally is provided in red and in the form x/30. This table shows the
efficiency of PBFS as it can solve all instances except two, for T = 3 and T = 4. Most importantly, the
average time to solve these instances is under 10 seconds for these problem sizes. Since many ports
today have a maximum tier requirement of 4 and need fast solutions, PBFS could be used in practice
in the case of small batches. However, for T = 5 and 6, PBFS cannot solve all instances optimally in a
timely manner. This suggests that, as the problem grows slightly, some instances become very hard to
solve, which should not be a surprise, knowing the NP-hardness of the problem. In order to avoid such
situations in real operations, heuristics can be used to provide a “good” sub-optimal solution (good in
the sense of being not too far from optimality). Therefore, we want to evaluate the performance of these
heuristics in order to know which one should be used in real operations.
34

S
5

6

7

8

9

10

T
Fill rate

3
50 percent

67 percent

4
50 percent

67 percent

5
50 percent

C
Solved
Time (s)

8

"

10

"

10

"

13

"

13

0.01

0.02

0.03

0.12

0.17

C
Solved
Time (s)

9

"

12

"

12

"

16

"

15

0.01

0.03

0.04

0.86

2.90

C
Solved
Time (s)

11

"

14

"

14

"

19

"

18

0.02

0.04

0.04

0.83

1.37

C
Solved
Time (s)

12

"

16

"

16

"

21

"

20

0.01

0.06

0.16

10.04

6.84

C
Solved
Time (s)

14

"

18

"

18

"

24

0.03

0.10

0.37

8.84

C
Solved
Time (s)

15

"

20

"

20

27
28/30

0.03

0.10

0.54

"

"

"
"

67 percent
17
28/30

6
50 percent
15

"

67 percent
20
15/30

5.17
20
25/30

18

"

24
14/30

15.94
23
24/30

21
23/30

28
5/30

27
20/30

24
22/30

32
5/30

23
29/30

30
10/30

27
19/30

36
2/30

25
29/30

34
12/30

30
22/30

40
2/30

"
"

Table 2: Instances solved by PBFS in the batch model with small batches.
We measure the performance of heuristics and the tightness of lower bounds in Tables 6 and 7.
Concerning lower bounds, b encompasses a significant number of relocations. Adding unavoidable
“bad” relocations in b1 and b2 , improves slightly the lower bound. But experiments seem to confirm
that b2 (n) − b1 (n) 6 b1 (n) − b(n) holds, supporting our intuition that the relative increase of lower bounds
bk (n) − bk−1 (n) decreases with k.
Concerning heuristics, EG and EM clearly outperform ERI as they result in lower expected numbers
of relocations. When we have access to the optimal solutions, both heuristics are on average at most 2%
more than the optimal solution. We expect this behavior to be similar for larger instances, however we
only have access to lower bounds to evaluate their performances. In this case, heuristics are on average
at most 11% more than b2 , hence at most 11% from the optimal solution (even though we believe that this
number is very conservative, as our lower bounds are not “tight”). Therefore, both EG and EM appear to
be good solutions for the SCRP under the batch model with small batches. In this case, we recommend
using EM for its simplicity of implementation and understandability.

5.2
5.2.1

Experiment 2: Batch model with larger batches
Modifying existing instances

For the sake of reproducibility, we use the existing set of instances, but slightly modify it to consider
larger batches. In order to create these instances, for each original instance n, consider n0 with the same
containers in the same configuration. But, if w is the batch of a container c in n, then we take the batch
l m
of c in n0 to be w0 = wγ , where γ > 1, i.e., we merge γ batches together. In these experiments, we take
35

γ = 2, which implies that batches have an average size of 4.

S
5

6

7

8

9

10

T
Fill rate

3
50 percent

67 percent

4
50 percent

67 percent

5
50 percent

C
Solved
Time (s)

8

"

10

"

10

"

13

"

13

0.08

0.29

0.14

4.55

3.20

C
Solved
Time (s)

9

"

12

"

12

"

16

"

15

0.08

0.47

0.25

126.37

14.74

C
Solved
Time (s)

11

"

14

"

14

19
29/30

18

0.13

0.71

0.58

C
Solved
Time (s)

12

"

16

"

16

0.08

1.67

1.26

C
Solved
Time (s)

14

"

18

"

18

0.13

1.49

1.47

C
Solved
Time (s)

15

"

20

"

20

0.17

0.79

3.10

"
"
"
"

"
"
"

67 percent

6
50 percent

17
21/30

"

15

67 percent
20
3/30

72.70
20
18/30

18
27/30

24
1/30

23
9/30

21
14/30

28
0/30

17.74
21
28/30

20
29/30

27
6/30

24
17/30

32
1/30

24
26/30

23
25/30

30
5/30

27
15/30

36
0/30

27
22/30

25
29/30

34
7/30

30
14/30

40
0/30

Table 3: Instances solved by PBFSA in the batch model with larger batches.

5.2.2

Target error 
b(n )

In order to set our target error, we consider the following. Let n0 be a given instance, and set  = 20 .
f (n )
In this case, we know that  6 2 0 , which implies that we are making an error of at most 50%. Note that
this error is very conservative due to 2 major things: first, b(n0 ) is not necessarily representative of f (n0 ),
specially if n0 has many containers. Second, the number of samples given by Hoeffding’s inequality is
also very conservative, probably making our approximation substantially more accurate than what we
can theoretically prove.
5.2.3

Results

"

Results are summarized in Tables 3. Similarly to Table 2, indicates that all 30 instances are solved
approximately by PBFSA within the given expected error. In this case, the average solution time to solve
these instances is given in seconds. Otherwise, the number of instances solved is provided in red and in
the form x/30. This table shows that PBFSA presents several advantages. First, it solves most of instances
with T = 4 and S 6 9 approximately within 2 minutes, while we note that PBFS was not able to solve most
of these. Moreover, as can be seen in Tables 8 and 9, PBFSA still outperforms the best heuristics despite
the fact that we only set the theoretical guarantee to 50% of optimality. Together, these two advantages
show the practicality of PBFSA for problem sizes typically encountered in real ports. Moreover, we note
that increasing the batch size appears to make the problem significantly more complicated to solve as
36

we can solve optimally larger instances in Experiment 1 (see Table 2). Finally, we remark that similar
conclusions of Experiment 1 can be drawn for lower and upper bounds (see Tables 8 and 9).

5.3

Experiment 3: Online model and comparison with Ku and Arthanari (2016)

Table 4 gives a summary similar to the two previous experiments. In addition, we report the results of
Ku and Arthanari (2016) who take a time limit of eight hours for each instance. In this table, X(X) indicates
that all 30 instances are solved optimally by both PBFS and Ku and Arthanari (2016). In this case, the
average solution time in seconds to solve these instances is given for PBFS and for Ku and Arthanari
(2016) in parenthesis. indicates that all 30 instances are solved optimally only by PBFS but not Ku and
Arthanari (2016). In this case, only the average solution time to solve these instances with PBFS is given
in seconds. Otherwise, the number of instances solved by PBFS is provided in red and in the form x/30.

"

S
5

6

7

8

9

10

T
Fill rate

3
50 percent

C
Solved
Time (s)

8
X(X)
0.01 (0.02)

C
Solved
Time

9
X(X)
0.01 (0.01)

C
Solved
Time (s)

11
X(X)
0.02 (0.33)

"

C
Solved
Time (s)

12
X(X)
0.01 (0.33)

"

16

C
Solved
Time (s)

14
X(X)
0.02 (32.24)

C
Solved
Time (s)

15
X(X)
0.03 (58.85)

67 percent
10

4
50 percent

67 percent

5
50 percent

10
X(X)
0.01 (2.51)

13

0.09

13
X(X)
0.16 (2483.30)

12
X(X)
0.06 (139.08)

16

"

15

0.81

1.85

14
X(X)
0.04 (207.62)

19

"

18

0.67

1.38

"

16

21

"

"

0.05

0.10

8.29

5.85

"

18

"

18

24

0.09

0.38

7.26

"

20

"

20

27
28/30

0.08

0.52

"

0.02
12

"

0.04
14

0.04

"

"

"

67 percent

6
50 percent

17
28/30

"

20
25/30

15

67 percent
20
18/30

3.74
18

"

24
15/30

14.92
23
24/30

21
23/30

28
5/30

27
20/30

24
22/30

32
5/30

23
29/30

30
10/30

27
19/30

36
2/30

25
29/30

34
12/30

30
16/30

40
2/30

"

20

Table 4: Instances solved by PBFS and Ku and Arthanari (2016) in the online model with small batch.
Results show strong evidence that our solution is improving significantly the best existing results for
the SCRP under the online model, given that we solve many larger instances optimally. Furthermore, it
also outperforms the most recent algorithm in solution time for the problem sizes it can solve. It appears
that, for problems for which we can solve all (or almost all) instances, most instances are “easy” to solve
as the algorithm finds a solution within seconds. However, as in the batch model, there exists some
instances for which the optimal solution still requires an exponential number of nodes, which makes our
algorithm not tractable.
In Tables 10 and 11, we also report in parenthesis the averages for ERI and Random found by Ku and
Arthanari (2016). The results for random are consistent. However, we find significant better results for our
implementation of ERI. This is unexpected since the only difference between the two implementations
37

is the use of lower bound b, when the configuration has less than S containers remaining. Nevertheless,
ERI should also be optimal in this case, as it reduces to heuristic L. So this should not affect the expected
number of relocations, and we cannot explain this difference. Finally, we point out that the results are
quite similar to those of Experiment 1. Indeed, the existing data set has relatively small batches (on
average 2 containers), which inherently makes the two models, batch and online, very close to each
other.

5.4

Experiment 4: Online model with a unique batch
T
Fill rate

3
50 percent

67 percent

4
50 percent

67 percent

5
50 percent

67 percent

6
50 percent

67 percent

5

C
PBFS
L.

8
2.08
2.08

10
3.33
3.33

10
3.54
3.54

13
6.53
6.52

13
6.56
6.57

17
12.05
12.06

15
9.13
9.14

20
17.28
17.29

6

C
PBFS
L.

9
2.10
2.10

12
4.04
4.04

12
4.23
4.23

16
8.02
8.02

15
7.01
7.01

20
13.48
13.48

18
10.73
10.73

24
20.13
20.13

7

C
PBFS
L.

11
2.69
2.69

14
4.60
4.61

14
4.82
4.82

19
9.55
9.55

18
8.58
8.58

23
14.75
14.74

21
12.22
12.22

28
23.14
23.14

8

C
PBFS
L.

12
2.61
2.62

16
5.19
5.19

16
5.51
5.51

21
9.96
9.95

20
9.12
9.12

27
17.75
17.75

24
13.83
13.83

32
26.04

9

C
PBFS
L.

14
3.31
3.31

18
5.72
5.72

18
6.10
6.10

24
11.58
11.58

23
10.89
10.89

30
19.15
19.14

27
15.40

36
28.84

10

C
PBFS
L.

15
3.36
3.36

20
6.38
6.38

20
6.68
6.67

27
12.98
12.99

25
11.13
11.13

34
22.07
22.06

30
16.87

40
31.77

S

Table 5: Instances solved with PBFS and heuristic L in the online model with a unique batch.
In this experiment, we consider the existing data set, but assign all containers into a unique batch
(W = 1). We consider the SCRP under online model, where containers are revealed one at a time.
Note that, in this case, each container is equally likely to be retrieved, and it is equivalent to know
no-information about containers relative retrieval order. For each instance, we solve it twice: first using
PBFS, and then using heuristic L, for which we sample 10000 scenarios (this is different from the 5000
samples considered in previous experiments). We report the results in Table 5. In this table, for each
problem size, we report the expected optimal number of relocations averaged over 30 instances. ‘-’
means that all 30 instances could not be solved optimally with PBFS within the given time limit of an
hour. Note that, the exepected number of relocations using heuristic L reported in this experiment might
be less than the one of PBFS; this is only due to the fact that we are sampling. Intuitively, L should be
the optimal solution in this setting, and this experiment shows strong evidence that the next conjecture
holds.

38

Conjecture 1. Consider a configuration n with a unique batch. Let f o (n) be the minimum expected number of
relocations to empty n under the online model, and f o,L (n) be the expected number of relocations performed by the
Leveling heuristic under the online model, then
f o (n) = f o,L (n).

(22)

This conjecture could also be made in the dynamic case, when containers arrive to be stacked. These
results would have a strong implication in the port operations: leveling configurations is the optimal policy
when minimizing relocations, and no-information is given in advance.

6

Discussion

Managing relocation moves is one of the main challenges in the storage yard of container terminals,
because it has a direct effect on the costs and efficiency of yard operations. The Container Relocation
Problem, notorious for its computational intractability, addresses this issue. In this paper, we extend the
CRP to the more practical case in which the retrieval order of containers is not known far in advance. First,
we introduce a new stochastic model, called the batch model, show the applicability of this model and
compare it theoretically with the existing model of Zhao and Goodchild (2010). Then, we derive lower
bounds and fast and efficient heuristics for the SCRP. Subsequently, we develop two novel algorithms
(PBFS and PBFSA) to solve the stochastic CRP in different settings. Efficiencies of all algorithms are
supported through computational experiments, for which all results are made available online. Finally,
using our solution methods and based on extensive experiments, we conjecture the optimality of the
simple Leveling heuristic in the online stochastic setting. More generally, the methods developed in this
paper apply to multi-stage stochastic optimization problems, where the number of stages is finite, the set
of feasible actions at each stage is finite, the objective function is bounded and bounds on the objective
function can be easily computed.
Future work could include the proof of Conjecture 1. Another important future work is the optimal
design of time windows for a TAS. On one hand, small time windows imply more information on the
retrieval sequence, hence higher operational efficiency of port operators. On the other hand, large time
windows insure higher flexibility for truck drivers and a high rate of on-time arrivals. In order to balance
this trade-off, one would need to quantify two important metrics with respect to the expected number
of relocations: the “value of information” and the assignment of containers to “wrong” batches. Finally,
in the grand scheme of port operations, the study of stacking and retrieving simultaneously, as well as
the extension in the row dimension of blocks is important for future studies of operations to take into
account.

Acknowledgements
The authors would like to thank anonymous reviewers and the editor for their suggestions that led
to the significant improvements of this paper.

39

References
Akyüz, M. H. and Lee, C.-Y. (2013). A mathematical formulation and efficient heuristics for the dynamic
container relocation problem. Naval Research Logistics, 61(2):101–118.
Bertsekas, D. P. (2005). In Dynamic Programming and Optimal Control, Third Edition, volume 1 of Lecture
Notes in Computer Science. Athena Scientific.
Bonney, J. (2015).
Us ports move toward truck appointment model (04.24.2016).
Information Handling Services IHS.
Retrieved from http://www.joc.com/port-news/us-ports/
port-new-york-and-new-jersey/us-ports-move-toward-truck-appointment-model_20150427.
html.
Borjian, S., Galle, V., Manshadi, V. H., Barnhart, C., and Jaillet, P. (2015a). Container relocation problem:
Approximation, asymptotic, and incomplete information. CoRR, abs/1505.04229. http://arxiv.org/
abs/1505.04229.
Borjian, S., Manshadi, V., Barnhart, C., and Jaillet, P. (2015b). Managing relocation and delay in container
terminals with flexible service policies. CoRR, abs/1503.01535. https://arxiv.org/abs/1503.01535.
Caserta, M., Schwarze, S., and Voß, S. (2012). A mathematical formulation and complexity considerations
for the blocks relocation problem. European Journal of Operational Research, 219(1):96–104.
Expósito-Izquierdo, C., Melián-Batista, B., and Moreno-Vega, J. M. (2015). An exact approach for the
blocks relocation problem. Expert Systems with Applications, 42(17-18):6408 – 6422.
Galle, V., Borjian, S., Manshadi, V., Barnhart, C., and Jaillet, P. (2016). An average-case asymptotic analysis
of the container relocation problem. Operations Research Letters, 44(6):723 – 728.
Giuliano, G. and O’Brien, T. (2007). Reducing port-related truck emissions: The terminal gate appointment system at the ports of los angeles and long beach. Transportation Research Part D, 12:460 – 473.
Kim, K. H. and Hong, G.-P. (2006). A heuristic rule for relocating blocks. Computers & Operations Research,
33(4):940–954.
Ku, D. and Arthanari, T. S. (2016). Container relocation problem with time windows for container
departure. European Journal of Operational Research, 252(3):1031 – 1039.
Lehnfeld, J. and Knust, S. (2014). Loading, unloading and premarshalling of stacks in storage areas:
Survey and classification. European Journal of Operational Research, 239(2):297–312.
Morais, P. and Lord, E. (2006). Terminal appointment system study. Technical Report, Transportation
Development Center of Transport Canada. Retrieved from https://www.tc.gc.ca/media/documents/
policy/14570e.pdf.
Murty, K. G., wah Wan, Y., Liu, J., Tseng, M. M., Lai, E. L. K.-K., and Chiu, H. W. (2005). Hongkong
international terminals gains elastic capacity using a data-intensive decision-support system. Interfaces,
35(1):61 – 75.
40

Petering, M. E. and Hussein, M. I. (2013). A new mixed integer program and extended look-ahead
heuristic algorithm for the block relocation problem. European Journal of Operational Research, 231(1):120–
130.
Phillips, E. E. (2015).
Southern california ports to try trucking appointment system
(08.27.2015).
The Wall Street Journal.
Retrieved from http://www.wsj.com/articles/
southern-california-ports-to-try-trucking-appointment-system-1440711102.
Sennott, L. I. (2009). Stochastic dynamic programming and the control of queueing systems, volume 504. John
Wiley & Sons.
Stahlbock, R. and Voß, S. (2008). Operations research at container terminals: a literature update. OR
Spectrum, 30(1):1–52.
Steenken, D., Voß, S., and Stahlbock, R. (2004). Container terminal operation and operations research a
classification and literature review. OR Spectrum, 26:3 – 49.
Tanaka, S. and Takii, K. (2014). A faster branch-and-bound algorithm for the block relocation problem.
2014 IEEE International Conference on Automation Science and Engineering (CASE), 14:7–12.
Tierney, K. and Voß, S. (2016). Solving the Robust Container Pre-Marshalling Problem, pages 131–145.
Springer International Publishing.
Ünlüyurt, T. and Aydın, C. (2012). Improved rehandling strategies for the container retrieval process.
Journal of Advanced Transportation, 46(4):378–393.
van Asperen, E., Borgman, B., and Dekker, R. (2013). Evaluating impact of truck announcements on
container stacking efficiency. Flexible Services and Manufacturing Journal, 25(4):543–556.
Wu, K.-C. and Ting, C.-J. (2010). A beam search algorithm for minimizing reshuffle operations at container
yards. International Conference on Logistics and Maritime Systems, September, pages 703 – 710.
Wu, K.-C. and Ting, C.-J. (2012). Heuristic approaches for minimizing reshuffle operations at container
yard. Proceedings of the Asia Pacific industrial engineering & management systems conference, pages 1407 –
51.
Zehendner, E., Casserta, M., Feillet, D., Schwarze, S., and Voß, S. (2015). An improved mathematical
formulation for the blocks relocation problem. European Journal of Operational Research, 245:415 – 422.
Zehendner, E. and Feillet, D. (2014a). Benefits of a truck appointment system on the service quality of
inland transport modes at a multimodal container terminal. European Journal of Operational Research,
235(2):461 – 469.
Zehendner, E. and Feillet, D. (2014b). A branch and price approach for the container relocation problem.
International Journal of Production Research, 52(24):7159 – 7176.
Zehendner, E., Feillet, D., and Jaillet, P. (2016). An algorithm with performance guarantee for the online
container relocation problem. European Journal of Operational Research, 245:415 – 422.
41

Zhao, W. and Goodchild, A. V. (2010). The impact of truck arrival information on container terminal
rehandling. Transportation Research Part E: Logistics and Transportation Review, 46(3):327–343.
Zhu, W., Qin, H., Lim, A., and Zhang, H. (2012). Iterative deepening A* algorithms for the container
relocation problem. IEEE Transactions on Automation Science and Engineering, 9(4):710–722.

A

A.1

Theoretical and computational comparison of the batch and the online
models
Theoretical comparison: proof of Lemma 1

Lemma 1. Let y be a given initial configuration, then we have


f y 6 fo y .
Proof. We prove this lemma by induction on the number of batches W. The lemma clearly holds if y is
empty (i.e. W = 0). Now consider W > 1 and C1 > 1. For the sake of clarity of the proof, we define the
following notation:










∀ d ∈ {1, . . . , C1 } , 









ζ1 ,...,ζd

y −−−−−→ xd1
ak

xdk −→ xdk+1 , if d > 1 , ∀ k ∈ {1, . . . , d − 1}
ak

xdk −→ ydk−d+2 , ∀ k ∈ {d, . . . , C}

(23)

ζk

ydk−d+1 −→ xdk , ∀ k ∈ {d + 1, . . . , C}.

These notations corresponds to the following process: the first d containers to be retrieved are all revealed
at once. Then decisions to retrieve these d containers are made. Afterwards, each of the C − d remaining
containers is revealed one at a time (as in the online model). Under this revelation process, the minimum
expected number of relocation is given by
fd


y =




d

X
 
 




d
o
d 
 , ∀ d ∈ {1, . . . , C1 }.
r
x
+
f
y
E  min 
2 

k


ζ1 ,...,ζd a1 ,...,ad 
k=1

Moreover, using the recursion formula from the online model, we have

n 
 

 o
o
d
d
o
f y2 = E min r xd+1 + f yd3 .
ζd

ad



In particular, by definition of the online model, we have f o y = f 1 y .
Using these relations, let us prove that


f d y 6 f d−1 y , ∀ d ∈ {2, . . . , C1 }.

42

(24)

Let d ∈ {2, . . . , C1 }, we have
fd




d


X








d
o
d 

y = E  min 
r
x
+
f
y

2


k

ζ1 ,...,ζd a1 ,...,ad 
k=1



 
d


X
 
 
 






d
o
d 
= E E  min 
min
r
x
+
f
y



2 
k
 ad 



ζ1 ,...,ζd−1 ζd a1 ,...,ad−1 
k=1
 



d


X













d−1
o
d−1 
min 

6 E  min 
E
r
x
+
f
y




3

k
 ζd  a d 



ζ1 ,...,ζd−1 a1 ,...,ad−1 
k=1



d−1 


X
n 


o






d−1
o
d−1
d−1
= E  min 
r
x
+
f
y
r
x
+
E
min

3

d
k
ad


ζ1 ,...,ζd−1 a1 ,...,ad−1 
ζd
k=1



d−1 

X








d−1
o
d−1 
 = f d−1 y ,
= E  min 
r
x
+
f
y

2


k

a
,...,a


ζ1 ,...,ζd−1 1
d−1


(25)

(26)

(27)

k=1

where the equality (27) holds since xd−1
for k ∈ {1, . . . , d − 1} does not depend on ad and ζd . Finally, the
k
inequality holds because we have E [min {Z1 , . . . , Zm }] 6 min {E [Z1 , . . . , Zm ]} for any Z1 , . . . , Zm random
variables. Note that we changed xdk in xd−1
and yd2 in yd−1
. This change is made in order to stay consistent
3
k
with the definition of Equation (23). Indeed, the order between the expectations and the minimums in
Equation (25) implies that the process of the first d retrievals corresponds to
ζ1 ,...,ζd

a1

ad−1

a2

ad

y −−−−−→ xd1 −→ xd2 −→ . . . −−−→ xdd −→ yd2 ,
while the order between the expectations and the minimums in Equation (26) corresponds to the following
process for the first d retrievals:
ζ1 ,...,ζd−1

a1

ad−1

a2

ζd

ad

y −−−−−−→ xd−1
−→ xd−1
−→ . . . −−−→ yd−1
−→ xd−1
−→ yd−1
2
2
3 .
1
d
Recall Equation (1) and apply it with w = 1 (note that K1 = 1 thus K1 + C1 − 1 = C1 ) to get

C

1


X







(x
)
.
r
+
f
y
f y = E  min 

2
k



ζ1 ,...,ζC1 a1 ,...,aC1 
k=1



By induction, for all configuration y2 with W − 1 batches we have f y2 6 f o y2 , thus
C

C



1
1



X
X
 C
 C 













o
 6 E  min 
1
1 

 = f C1 y ,
(x
)
f y = E  min 
r
+
f
y
r
x
+
f
y


2


k




2
k



a
,...,a
a
,...,a




ζ1 ,...,ζC1 1
ζ1 ,...,ζC1 1
C1
C1
k=1

k=1

where we replaced xk by xCk 1 and y2 by yC2 1 because, on the right-hand-side of the inequality, the revelation


process after the first C1 containers is the online model. Finally, since f o y = f 1 y , by applying Equation
(24) for each value of d ∈ {C1 , . . . , 2}, we complete the proof as






f y 6 f C1 y 6 f C1 −1 y 6 . . . 6 f 2 y 6 f 1 y = f o y .

As a final remark, Lemma 1 is tight in the general setting. Indeed, there exists an initial configuration
y for which f (y) = f o (y). For instance consider the configuration in Figure 4a, then we have f (y) =
f o (y) = 13/6.
43

A.2

Computational comparison

There also exist configurations for which f (y) < f o (y). The difference between these two values
represents the value of taking into account available information (if possible).
In order to show a positive difference, we could have compared Experiments 1 and 3. However, since
the average batch size is 2, these experiments do not show a positive difference between both models.
Another possibility would have been to use the instances of Experiment 2. However, as we previously
mentioned, such instances are hard to solve optimally and not approximately.
Instead, we consider another set of simpler instances randomly generated: 100 instances with T = 4
tiers, S = 4 stacks and C = 12 containers. Each instance has W = 3 batches and each batch has Cw = 4
containers (for w = 1, 2, 3). For each of these 100 instances, we solve it under the batch and the online
models. The code and detailed results are available at https://github.com/vgalle/StochasticCRP.
f o (.) − f (.)
We are especially interested about
× 100, which the % difference between the batch and the
f (.)
online models.
On average over the 100 instances, the optimal expected number of relocations under the batch model
is 6.526 and under the online model is 6.61, hence giving a difference of 0.084. We observe here that
this difference represents more than 1.287% of the optimal solution under the batch model, which is
quite significant considering the fact that heuristic EM experimentally lies within 2% above the optimal
solution. In addition, we noticed that for 25 of these instances, this difference was more than 2% and the
maximum was about 4% (see Figure 14a).
0

0.5

1

1.5

2

2.5

3

3.5

0

4

14a T = 4, S = 4, C = 12, W = 3 and Cw = 4
(for w = 1, 2, 3).

2

4

6

8

10

14b T = 4, S = 4, C = 12, W = 2 and Cw = 6
(for w = 1, 2).

Figure 14: Distributions of % difference between the batch and the online models from 100 randomly
generated instances.
We also consider 100 instances for which T = 4, S = 4, C = 12, but now W = 2 and each batch has
Cw = 6 (for w = 1, 2). Figure 14b show that this relative difference appears to increase when the batch
size increases. Indeed, the average difference is about 4.251% (batch: 6.751, online: , 7.038, difference:
0.287) with 25 instances having a difference of more than 5.3%.

44

B

Computational Experiments Tables

S

T

C

Lower bounds
b
b1
b2

PBFS

Heuristics
EG
EM

ERI

L

Rand.

5

3
4
5
6

8
10
13
15

1.64
2.88
4.61
6.28

1.66
2.96
4.88
6.64

1.66
2.99
5.01
7.06

1.70
3.11
5.32
7.59

1.70
3.11
5.40
7.85

1.70
3.13
5.38
7.81

1.70
3.14
5.57
8.09

1.82
3.51
6.17
9.41

2.34
4.62
8.00
12.35

6

3
4
5
6

9
12
15
18

1.68
3.54
5.37
7.19

1.69
3.61
5.57
7.52

1.69
3.63
5.68
7.68

1.74
3.68
5.91
8.23

1.76
3.69
5.97
8.38

1.74
3.68
5.94
8.29

1.74
3.68
6.01
8.63

1.84
4.11
7.21
10.05

2.43
5.59
9.55
13.53

7

3
4
5
6

11
14
18
21

2.82
3.97
6.49
8.82

2.86
4.06
6.65
9.21

2.88
4.10
6.74
9.51

2.88
4.16
6.97
-

2.88
4.17
7.05
10.40

2.88
4.17
7.00
10.35

2.89
4.20
7.07
10.76

2.96
4.65
8.46
12.47

4.12
6.47
11.27
17.69

8

3
4
5
6

12
16
20
24

2.29
4.68
7.20
9.52

2.30
4.73
7.42
9.85

2.30
4.75
7.54
10.09

2.31
4.82
7.85
-

2.31
4.83
7.96
11.10

2.31
4.83
7.93
10.99

2.31
4.83
8.06
11.34

2.43
5.41
9.32
13.29

3.23
7.48
13.44
19.28

9

3
4
5
6

14
18
23
27

2.98
5.63
8.58
10.38

2.98
5.71
8.69
10.78

2.98
5.71
8.77
10.98

3.00
5.73
-

3.00
5.73
9.05
11.59

3.01
5.73
9.02
11.58

3.00
5.73
9.12
11.76

3.19
6.52
11.16
14.62

4.54
9.29
15.57
20.93

10

3
4
5
6

15
20
25
30

3.18
6.20
9.10
11.91

3.18
6.23
9.37
12.28

3.18
6.23
9.39
12.44

3.19
6.28
-

3.19
6.30
9.60
13.01

3.20
6.28
9.60
12.92

3.20
6.28
9.73
13.15

3.27
6.98
11.38
15.93

4.75
10.41
16.64
23.40

Table 6: Results of experiment 1: Performance of PBFS, heuristics and tightness of lower bounds for a
fill rate of 50 percent in the Batch Model, in the case of small batches. Bold numbers highlight the best
heuristic for a given problem size.

45

S

T

C

Lower bounds
b
b1
b2

PBFS

Heuristics
EG
EM

ERI

L

Rand.

5

3
4
5
6

10
13
17
20

2.83
4.69
7.58
9.69

2.97
4.97
8.52
10.79

3.01
5.09
8.92
11.46

3.08
5.58
-

3.08
5.69
10.44
14.53

3.08
5.64
10.48
14.65

3.11
5.75
11.04
15.77

3.33
6.50
12.23
18.47

4.16
8.22
15.27
22.93

6

3
4
5
6

12
16
20
24

3.60
6.20
8.28
11.67

3.70
6.59
8.85
12.22

3.72
6.78
9.16
12.54

3.89
7.28
-

3.89
7.41
10.39
15.19

3.90
7.45
10.38
15.17

3.90
7.61
10.80
16.14

4.32
8.55
12.85
19.33

5.54
11.29
16.43
25.03

7

3
4
5
6

14
19
23
28

3.85
6.25
9.72
13.52

3.89
6.60
10.24
14.41

3.91
6.86
10.55
14.93

3.97
7.29
-

3.98
7.37
11.75
17.72

3.98
7.36
11.71
17.63

4.02
7.54
12.30
18.81

4.40
8.89
14.92
23.10

6.05
11.60
19.65
30.78

8

3
4
5
6

12
21
27
32

4.47
7.62
11.61
15.60

4.57
7.85
12.08
16.39

4.61
7.98
12.52
16.78

4.66
8.29
-

4.66
8.33
13.56
19.28

4.66
8.35
13.47
19.51

4.68
8.43
14.10
20.85

5.14
9.76
17.15
25.89

6.94
13.26
23.11
34.66

9

3
4
5
6

18
24
30
36

4.81
8.98
13.16
16.77

4.96
9.18
13.90
17.36

4.99
9.30
14.29
17.83

5.10
9.58
-

5.10
9.63
15.65
20.38

5.12
9.61
15.79
20.40

5.14
9.76
16.75
21.86

5.66
11.66
20.03
28.12

7.81
16.00
27.41
38.12

10

3
4
5
6

20
27
34
40

5.21
9.18
14.46
19.55

5.21
9.54
14.88
20.24

5.21
9.71
15.16
20.66

5.27
-

5.28
10.27
16.13
23.33

5.28
10.29
16.19
23.20

5.28
10.37
16.69
24.46

5.79
12.16
21.06
32.11

7.86
16.91
29.03
44.07

Table 7: Results of experiment 1: Performance of PBFS, heuristics and tightness of lower bounds for a
fill rate of 67 percent in the Batch Model, in the case of small batches. Bold numbers highlight the best
heuristic for a given problem size.

46

S

T

C

Lower bounds
b
b1
b2

PBFSA

Heuristics
EG
EM

ERI

L

Rand.

5

3
4
5
6

8
10
13
15

1.68
2.96
4.58
6.31

1.68
3.02
4.82
6.72

1.68
3.05
4.93
6.96

1.76
3.33
5.50
7.81

1.77
3.36
5.69
8.28

1.76
3.37
5.65
8.03

1.76
3.38
5.74
8.35

1.87
3.67
6.33
9.49

2.38
4.85
8.22
12.40

6

3
4
5
6

9
12
15
18

1.66
3.61
5.38
7.01

1.67
3.71
5.57
7.26

1.67
3.73
5.64
7.44

1.73
3.93
6.11
-

1.75
3.99
6.23
8.49

1.74
3.99
6.23
8.36

1.74
4.00
6.31
8.61

1.82
4.34
7.16
9.92

2.43
5.82
9.65
13.45

7

3
4
5
6

11
14
18
21

2.76
4.02
6.29
8.69

2.79
4.12
6.39
9.12

2.79
4.15
6.43
9.35

2.85
4.24
6.77
-

2.84
4.31
7.00
10.60

2.84
4.29
6.92
10.52

2.83
4.31
7.01
10.91

2.95
4.73
8.20
12.61

4.06
6.52
11.08
17.79

8

3
4
5
6

12
16
20
24

2.30
4.61
7.31
9.65

2.31
4.62
7.46
9.95

2.31
4.63
7.52
10.12

2.31
4.71
-

2.31
4.74
8.01
11.37

2.32
4.74
8.01
11.37

2.32
4.75
8.09
11.67

2.37
5.25
9.33
13.44

3.19
7.40
13.25
19.51

9

3
4
5
6

14
18
23
27

2.93
5.56
8.49
10.38

2.93
5.58
8.64
10.69

2.93
5.59
8.73
10.90

2.95
5.69
-

2.96
5.74
9.16
11.77

2.96
5.70
9.12
11.75

2.96
5.70
9.16
11.95

3.15
6.33
10.99
14.65

4.48
9.07
15.39
20.95

10

3
4
5
6

15
20
25
30

3.15
6180
9.13
12.09

3.16
6.20
9.31
12.35

3.16
6.21
9.36
12.51

3.15
6.28
-

3.17
6.35
9.66
13.38

3.17
6.34
9.63
13.23

3.17
6.34
9.68
13.42

3.25
6.92
11.44
16.33

4.72
10.27
16.73
23.80

Table 8: Results of experiment 2: Performance of PBFSA, heuristics and tightness of lower bounds for a
fill rate of 50 percent in the Batch Model with larger batches. Bold numbers highlight the best heuristic
for a given problem size.

47

S

T

C

Lower bounds
b
b1
b2

PBFSA

Heuristics
EG
EM

ERI

L

Rand.

5

3
4
5
6

10
13
17
20

2.78
4.71
7.53
9.69

2.87
4.92
8.17
10.56

2.90
5.00
8.44
11.12

3.07
5.70
-

3.08
5.81
10.38
15.00

3.08
5.80
10.32
14.91

3.09
5.89
10.68
16.01

3.36
6.60
11.96
18.22

4.18
8.23
14.9815
22.79

6

3
4
5
6

12
16
20
24

3.56
6.12
8.33
11.77

3.63
6.48
8.72
12.41

3.64
6.61
8.90
12.80

3.90
7.17
-

3.90
7.36
10.36
15.94

3.91
7.41
10.30
15.91

3.92
7.48
10.57
16.86

4.28
8.44
12.54
19.83

5.55
11.02
16.15
25.45

7

3
4
5
6

14
19
23
28

3.95
6.27
9.81
13.64

4.01
6.56
10.27
14.47

4.02
6.77
10.54
14.91

4.14
-

4.17
7.36
12.08
18.36

4.15
7.35
12.10
18.26

4.15
7.52
12.49
19.32

4.56
8.66
14.78
23.10

6.11
11.46
19.67
31.03

8

3
4
5
6

12
21
27
32

4.65
7.58
11.46
15.45

4.74
7.86
11.98
16.28

4.76
7.99
12.28
16.72

4.85
-

4.86
8.38
13.73
19.94

4.85
8.42
13.60
19.83

4.88
8.50
14.11
21.21

5.34
9.82
17.00
25.73

7.14
13.20
22.84
34.70

9

3
4
5
6

18
24
30
36

4.85
8.82
13.15
16.85

4.98
9.00
13.84
17.39

5.02
9.11
14.18
17.80

5.14
-

5.17
9.66
15.91
20.99

5.20
9.60
15.96
20.83

5.20
9.72
16.82
21.97

5.67
11.52
20.16
28.05

7.77
15.70
27.35
38.04

10

3
4
5
6

20
27
34
40

5.19
9.40
14.44
19.49

5.21
9.66
14.83
20.24

5.22
9.82
15.07
20.66

5.31
-

5.31
10.47
16.29
23.83

5.31
10.46
16.29
23.65

5.31
10.54
16.62
24.96

5.79
12.25
21.12
32.04

7.92
16.96
28.86
44.12

Table 9: Results of experiment 2: Performance of PBFSA, heuristics and tightness of lower bounds for a
fill rate of 67 percent in the Batch Model with larger batches. Bold numbers highlight the best heuristic
for a given problem size.

48

S

T

C

Lower bounds
b
b1
b2

PBFS

Heuristics
EG
EM

ERI

L

Rand.

5

3
4
5
6

8
10
13
15

1.64
2.88
4.61
6.28

1.66
2.96
4.88
6.64

1.66
2.99
5.01
7.06

1.70
3.11
5.32
7.59

1.70
3.11
5.38
7.85

1.70
3.13
5.38
7.80

1.71 (1.71)
3.14 (3.20)
5.57 (5.58)
8.08 (8.29)

1.82
3.51
6.16
9.41

2.34 (2.34)
4.62 (4.62)
8.00 (8.00)
12.36 (12.35)

6

3
4
5
6

9
12
15
18

1.68
3.54
5.37
7.19

1.69
3.61
5.57
7.52

1.69
3.63
5.68
7.68

1.74
3.68
5.91
8.23

1.76
3.69
5.96
8.38

1.74
3.68
5.94
8.29

1.74 (1.75)
3.68 (3.75)
6.00 (6.18)
8.62 (8.77)

1.84
4.11
7.21
10.05

2.43 (2.43)
5.59 (5.59)
9.54 (9.54)
13.53 (13.53)

7

3
4
5
6

11
14
18
21

2.82
3.97
6.49
8.82

2.86
4.06
6.65
9.21

2.88
4.1
6.74
9.51

2.88
4.16
6.97
-

2.88
4.17
7.04
10.40

2.88
4.17
7.00
10.35

2.89 (2.88)
4.21 (4.20*)
7.07 (7.18)
10.76 (10.98)

2.96
4.66
8.45
12.46

4.11 (4.11)
6.47 (6.03*)
11.27 (11.27)
17.69 (17.69)

8

3
4
5
6

12
16
20
24

2.29
4.68
7.20
9.52

2.3
4.73
7.42
9.85

2.3
4.75
7.54
10.09

2.31
4.82
7.85
-

2.31
4.83
7.97
11.10

2.31
4.83
7.94
10.98

2.31 (2.32)
4.83 (4.88)
8.07 (8.27)
11.34 (11.61)

2.43
5.41
9.32
13.29

3.23 (3.23)
7.49 (7.49)
13.44 (13.45)
19.29 (19.29)

9

3
4
5
6

14
18
23
27

2.98
5.63
8.58
10.38

2.98
5.71
8.69
10.78

2.98
5.71
8.77
10.98

3.00
5.73
-

3.00
5.73
9.05
11.59

3.00
5.73
9.02
11.58

3.00 (3.00)
5.73 (5.80)
9.12 (9.36)
11.76 (12.09)

3.19
6.52
11.16
14.62

4.54 (4.54)
9.29 (9.29)
15.56 (15.57)
20.94 (20.93)

10

3
4
5
6

15
20
25
30

3.18
6.20
9.10
11.91

3.18
6.23
9.37
12.28

3.18
6.23
9.39
12.44

3.19
6.28
-

3.19
6.30
9.61
13.01

3.20
6.27
9.60
12.92

3.20 (3.20)
6.28 (6.33)
9.73 (9.80)
13.15 (13.51)

3.27
6.98
11.38
15.92

4.75 (4.75)
10.41 (10.41)
16.64 (16.63)
23.41 (23.41)

Table 10: Results of experiment 3: Performance of heuristics and tightness of lower bounds for a fill rate
of 50 percent in the online model with small batches. Bold numbers highlight the best heuristic for a
given problem size. Numbers in parenthesis are taken from Ku and Arthanari (2016).

49

S

T

C

Lower bounds
b
b1
b2

PBFS

Heuristics
EG
EM

ERI

L

Rand.

5

3
4
5
6

10
13
17
20

2.83
4.69
7.58
9.69

2.97
4.97
8.52
10.79

3.01
5.09
8.92
11.46

3.08
5.58
-

3.08
5.68
10.45
14.53

3.08
5.64
10.48
14.65

3.12 (3.10)
5.75 (5.80)
11.04 (11.15)
15.77 (16.14)

3.33
6.50
12.24
18.46

4.16 (4.16)
8.22 (8.22)
15.28 (15.28)
22.93 (22.93)

6

3
4
5
6

12
16
20
24

3.6
6.2
8.28
11.67

3.7
6.59
8.85
12.22

3.72
6.78
9.16
12.54

3.89
7.28
-

3.89
7.41
10.38
15.17

3.90
7.45
10.38
15.17

3.90 (3.92)
7.61 (7.68)
10.80 (10.97)
16.14 (16.65)

4.32
8.54
12.85
19.33

5.53 (5.53)
11.29 (11.28)
16.42 (16.42)
25.04 (25.03)

7

3
4
5
6

14
19
23
28

3.85
6.25
9.72
13.52

3.89
6.6
10.24
14.41

3.91
6.86
10.55
14.93

3.97
7.29
-

3.98
7.37
11.76
17.70

3.98
7.36
11.71
17.64

4.02 (4.01)
7.54 (7.68)
12.30 (12.64)
18.82 (19.49)

4.40
8.89
14.92
23.10

6.05 (6.05)
11.60 (11.61)
19.66 (19.65)
30.77 (30.79)

8

3
4
5
6

12
21
27
32

4.47
7.62
11.61
15.6

4.57
7.85
12.08
16.39

4.61
7.98
12.52
16.78

4.66
8.29
-

4.65
8.32
13.56
19.27

4.66
8.35
13.47
19.51

4.68 (4.7)
8.43 (8.5)
14.10 (14.44)
20.85 (21.72)

5.14
9.75
17.14
25.89

6.94 (6.94)
13.26 (13.25)
23.11 (23.12)
34.64 (34.63)

9

3
4
5
6

18
24
30
36

4.81
8.98
13.16
16.77

4.96
9.18
13.9
17.36

4.99
9.3
14.29
17.83

5.10
9.58
-

5.10
9.63
15.65
20.38

5.12
9.61
15.79
20.40

5.14 (5.19)
9.76 (9.92)
16.75 (16.97)
21.87 (22.73)

5.66
11.66
20.03
28.13

7.80 (7.80)
16.01 (16.00)
27.38 (27.39)
38.11 (38.14)

10

3
4
5
6

20
27
34
40

5.21
9.18
14.46
19.55

5.21
9.54
14.88
20.24

5.21
9.71
15.16
20.66

5.27
-

5.28
10.27
16.13
23.33

5.28
10.29
16.19
23.20

5.28 (5.30)
10.37 (10.50)
16.69 (17.23)
24.46 (25.58)

5.79
12.15
21.07
32.11

7.85 (7.86)
16.92 (16.91)
29.03 (29.03)
44.08 (44.07)

Table 11: Results of experiment 3: Performance of heuristics and tightness of lower bounds for a fill rate
of 67 percent in the online model with small batches. Bold numbers highlight the best heuristic for a
given problem size. Numbers in parenthesis are taken from Ku and Arthanari (2016).

50

