Nested Column Generation Applied to the Crude Oil Tanker Routing and
Scheduling Problem with Split Pickup and Split Delivery
Frank Hennig,1 Bjørn Nygreen,1 Marco E. Lübbecke2
1

Department of Industrial Economics and Technology Management, Norwegian University of Science and Technology,
NO-7491 Trondheim, Norway
2

Chair of Operations Research, RWTH Aachen University, D-52072 Aachen, Germany

Received 25 January 2010; revised 11 March 2012; accepted 12 March 2012
DOI 10.1002/nav.21489
Published online 13 April 2012 in Wiley Online Library (wileyonlinelibrary.com).

Abstract: The split pickup split delivery crude oil tanker routing and scheduling problem is a difﬁcult combinatorial optimization
problem, both theoretically and practically. However, because of the large expenses in crude oil shipping it is attractive to make useof optimization that exploits as many degrees of freedom as possible to save transportation cost. We propose a nested column generation algorithm for this particular split pickup split delivery problem which bears several complexities such as a heterogeneous ﬂeet,
multiple commodities, many-to-many relations for pickup and delivery of each commodity, sequence dependent vehicle capacities,
and cargo quantity dependent pickup and delivery times. Our approach builds on a branch-and-price algorithm in which the column
generation subproblems are solved by branch-and-price themselves. We describe our implementation in the branch-cut-and-priceframework SCIP and give computational results for realistic test instances. The high quality schedules we obtain for these instances
improve on those in previous studies. © 2012 Wiley Periodicals, Inc. Naval Research Logistics 59: 298–310, 2012
Keywords: maritime transportation; crude oil; split pickup; split delivery; routing and scheduling; branch-and-price; nested
column generation

1.

INTRODUCTION

In this study, we propose an exact solution approach to areal-world maritime transportation problem called the “crude
oil tanker routing and scheduling problem.” This is a particularly rich split pickup split delivery problem with time
windows and capacity constraints and has been introducedin [10].
Maritime crude oil transportation is a major internationaltransportation mode. In 2010 the crude oil tanker tonnage
represented 34% of the total world ﬂeet capacity. The total
volume of oil shipped amounted to 2.75 billion tons [17].
A typical vessel engaged in long range, large volume crude
oil transportation is the very large crude oil carrier (VLCC).
These ships have a length of about 330 meters, a breadthof 60 meters, and can transport 300,000 tons of crude oil atonce. Each day by which the duration of a given transportation task may be shortened would save about 50,000 USD attoday’s fuel price of around 650 USD/ton. A port visit avoided
through good routing and scheduling could reduce total costby 100,000 USD [15]. Through more efﬁcient transportation
Correspondence to: F. Hennig (frank.hennig@iot.ntnu.no)© 2012 Wiley Periodicals, Inc.

ﬂeet utilization increases, more cargo can be transported andless harmful emissions per ton mile have to be accepted.
In the type of vehicle routing problem (VRP) considered,vehicles are ships and locations to be serviced are ports. In
contrast to the type of VRP, where vehicles operate from a
depot and do either pickup or delivery, in maritime crude oiltransportation the usually heterogeneous ships start at theirpresent location, pick up and deliver along their route and are
not scheduled for a certain destination. Port visits for loadingor discharging are time constrained and physical and regulatory circumstances set different limits for ship capacity on
different sailing legs. The cargo to be transported is liquidbulk of varying quality or type. Same qualities or types aresupplied and demanded in different locations so that a cargocannot be speciﬁed with is origin and destination alone. Alsoa given transportation task is not necessarily associated witha given shipment volume. Ships may load the quantities they
are suited for and that match an efﬁcient routing. This involves
the possibility of splitting supplied and demanded quantities
in arbitrary ways.
The contribution of our article is twofold. First, we present
an accurate model for a logistics problem of major practical importance and develop a solution approach to obtain

Hennig, Nygreen, and Lübbecke: Split Pickup and Split Delivery

industrially relevant solutions. Second, and probably more
important, we propose a general methodology to exactly
solve split pickup split delivery problems (SPSDP). The splitting of a customer’s demand between several vehicles of agiven ﬂeet is a natural attempt to realize the full optimization potential. This is accounted for in the vehicle routing
literature by studying the split delivery vehicle routing problem (SDVRP), for which a few exact methods were proposed [3, 5, 6, 9, 12, 13, 18]. The literature on split pickup
split delivery problems is even scarcer. We are aware of only
two related papers [14, 16], however, both fail to provide asolution approach for our problem. One [14] does not provide us with a state-of-the-art solution approach. The otherone [16] covers only splits for explicitly paired pickups anddeliveries in a much less complex environment.
Given the limited knowledge on combined split pickup
with split delivery, we build on the two recent papers by
Desaulniers [6] and Ceselli et al. [5] who propose branchand-price algorithms for the split delivery vehicle routing
problem. These authors exploit one crucial property of the
split delivery VRP: There is at most one customer per vehicle route for which only a fraction of the demand is served.
Whereas this is no longer true for pickup and delivery problems which permit splitting in both the pickup and delivery
locations, we are able to generalize the rationale behind theapproach in [5, 6].
Our approach is based on integer programming, and wewill make heavy use of branch-and-price [2, 8], a methodology we assume the reader to be reasonably familiar with. In
particular, we apply a nested column generation approach, aterm used in [19] and others to describe a column generation algorithm in which the subproblem is solved by column
generation itself. In this spirit, we solve a column generation master problem that comprises all ships in the ﬂeet andsolve subproblems, one for each ship, again using column
generation. Throughout this article, we will call the restricted
all-ship master problem level 1 RMP and the restricted singleship master problems level 2 RMP. The term subproblem isused to refer to the level 2 RMP route generation subproblem
solved by dynamic programming.
The article is structured as follows: In Section 2, we
describe the mathematical model thereby introduce somenotation. Section 3 motivates and describes the nested column generation concept and problem decomposition. The
computational algorithm and its implementation in the SCIP
framework is explained in Section 4. Computational results
are shown and discussed in Section 5. Finally, we draw
conclusions and suggest further research in Section 6.

2.

MODEL DESCRIPTION

A heterogeneous ﬂeet V of vessels or ships has to distribute crude oil between a set of ports. Each port supplies and/or

299

demands one or several types of crude oil, so-called grades C.
Each supply/demand has to be picked up or delivered withina certain time window. There may be several distinct time
windows to be serviced per port. Denote the set of all time
windows of all ports by N . For each ship v ∈ V, the set Nv
comprises N , and in addition a time window at the initial andat the ﬁnal location of v, respectively. Formally, the routingof v ∈ V takes place on a directed graph G = (Nv , A) witharc set A.
Each time window i ∈ N calls either for a pickup or adelivery, and is associated with a (supply or demand) quantity Qi of a particular crude grade c ∈ C and a time interval[T i , T i ] in which service has to start. We will say that a shipvisits a time window, meaning that it visits the associatedport within this time window. Supply and demand quantities are speciﬁed in weight units and may exceed the ship’s
capacity. There can be several combinations of pickup time
windows to supply a given delivery time window and viceversa. Pickup and delivery ports are concentrated in certain pickup or delivery regions, respectively. Large distances
between pickup and delivery regions lead to ship routes having what is called a “voyage” structure: On a voyage a ship
may visit one or several pickup regions followed by one orseveral delivery regions. It is unrealistic that a ship carriescargo from a delivery region back to any pickup region. A ship
may undertake several voyages during the planning period,
voyage length permitting. The number of time windows aship services in a region is limited by practical considerations.
Therefore the maximum number of different grades simultaneously onboard is limited as well. Because pickup anddelivery times may amount to several days for large tanker
ships, the quantities picked up and delivered may inﬂuence
the time feasibility of a route. For time feasibility, it is sufﬁcient that a ship arrives at time window i no later than closing
time T i . A maximum waiting time between arrival at time
window i and the time window’s opening time T i may bedesired. Sailing times TijSv for ship v on arcs (i, j ) ∈ A are
assumed to be deterministic. Ships have a volume capacity
based on cargo tank size and a weight capacity due to seaworthiness conditions. The crude grade speciﬁc density Dc ,c ∈ C, allows a conversion between weight quantities andvolumes. The ship capacity may change during the executionof a route due to external factors like, for example, limited
water depth and port regulations. An arc (i, j ) may impose astricter weight capacity limit W ij v , a stricter volume capacity limit V ij v , or both, for ship v. Hence, pickup and delivery
quantities can be inﬂuenced by the order in which different time windows are visited. Set AWv contains all arcs that
impose a weight limit on ship v; set AVv contains all arcs
with a volume limit for v. Arcs may be in both sets. A ship
is allowed to serve any fraction of a time window quantity.
Time for service can be calculated based on the per weight
unit time consumption TiQ . Minimum pickup quantities may
apply.
Naval Research Logistics DOI 10.1002/nav

300

Naval Research Logistics, Vol. 59 (2012)

3.

NESTED DECOMPOSITION
3.1.

Motivation

Recently, Desaulniers [6] and Ceselli et al. [5] very successfully applied branch-and-price algorithms to the SDVRP.
They both formulate master problems where each variablecorresponds to a combination of a feasible route and a cargopattern, where the latter represents compatible quantities ofdemand served for each customer on the route. Building on
well-known properties of the SDVRP, they make use of thefact that it is sufﬁcient to consider only cargo patterns where
the demand of at most one customer (the split customer) isserved fractionally, and all other customers are served either
fully or not at all (an analogy holds for the fractional knapsackproblem as well). In the language of Dantzig-Wolfe decompostion this means that an extreme point of the polyhedron
describing feasible cargo patterns has at most one fractional
component. We refer to the quantities associated with these
extreme points as extreme cargo patterns. Given a route anda split customer, only two of these patterns are necessary
to convex combine all possible fractional quantities in themaster problem. This knowledge is instrumental in designing an algorithm for the pricing problem which becomes anelementary shortest path problem (ESPP) where the numberof split customers is limited to one. In fact, this adds a helpful“knapsack ﬂavor” to the pricing problem.
The situation is a bit more complicated for the SPSDVRP.
The slightly more general statement is still true that in an
extreme cargo pattern there is at most one pickup location
and at most one delivery location for which the demand isserved fractionally. Principally, the label setting algorithm for
the ESPP in [6] should not be too difﬁcult to adapt respecting
this limitation to two split customers. In contrast, dominancebetween labels, that is, early elimination of unpromising partial routes is already more complicated. This is because it
relies on being able to compare transported quantities enroute, and two split locations give more freedom and thus
a weaker dominance. In addition, we consider a multicommodity problem, that is, we may have a split pickup andsplit delivery location for each commodity. As in the ordinary pickup and delivery problem, we may ﬁnd ourselves inthe situation that we cannot decide on feasibility of a route
before we see it in total, and this may happen even more frequently with several commodities. This already leads us to
not modifying the existing label setting algorithm because
it is unlikely that we are able to handle the combinatorial
explosion. Even if we had some hope that a combinatorial
algorithm would be the right choice, the weight and volume
constraints on arcs, particular to the COTRASP, ﬁnally burythis. Not only feasibility poses problems, the proﬁts for serving a particular time window can depend on the entire route.
There are situations in which it is neither most proﬁtable
to go to the volume nor to the weight limit of a particular
Naval Research Logistics DOI 10.1002/nav

commodity, see Fig. 1 for an example. Then, reduced costs
depend on the amounts loaded and unloaded further on theroute, and cannot be computed locally. This, in fact, destroys
the “fractional knapsack character” of the pricing problem as
exploited in [5, 6]. Needless to say at this point, that loaded
and unloaded quantities have an inﬂuence on service times
and thus on time window feasibility later on the route. These
complications together render dominance effectively uselessand a label-setting dynamic program is unlikely to be efﬁcient. This is consistent with observations in [4] in a similar
situation.
To overcome these difﬁculties we decided to decouple therouting from the service, that is, we ﬁrst construct an entire
route and then ﬁnd corresponding quantities to be served ateach customer. We retain the key observation by Desaulniers[6] and Ceselli et al. [5] to use extreme cargo patterns, but weneed to obtain them with the help of a mixed integer program(MIP). Whereas in the SDVRP a convex combination of only
two extreme patterns is sufﬁcient for a given route and a split
customer, we may need many more extreme patterns in our
SPSDVRP. Candidates are all the extreme points of a polyhedron describing the feasible quantities, and these are far too
many to explicitly list. Describing feasible cargo patterns via
a polyhedron allows us to formulate even complicated side
constraints; but we need to use column generation to generate
the extreme cargo patterns for a given route. On the positive side, we can use a straightforward implementation of an
ESPP algorithm for routing. In that sense, the decomposition
also leads to simplicity and manageability of our approachin the ﬁrst place, plus it adds a lot more ﬂexibility. As said,on the downside, the overall algorithmic structure gets a little
more complicated because we will nest two column generation algorithms. It should be noted that the approach is exact
in principle, even though we will use it as a heuristic. The
details are given in the following subsection.
3.2.

Concept: Three Levels

The algorithm for solving the COTRASP model is basedon a decomposition into three levels: On the bottom level
we have a route generation subproblem for each single ship.
An intermediate level master problem computes (usually several) cargo patterns for each route. Route and cargo patterns
together form what we call cargo routes, and the cargo routes
from different ships are coordinated on the top level masterproblem. From top to bottom we number the levels 1 through3. Although the feasibility of routes, schedules and cargo patterns is ensured on levels 2 and 3, the correct time window
quantity split is decided on level 1 via a convex combinationof cargo patterns. The restricted master problems at levels 1
and 2 are called RMP 1 and RMP 2, respectively. The following three subsections describe the three levels in more detail.
More background on the underlying practical problem canbe found in [10].

Hennig, Nygreen, and Lübbecke: Split Pickup and Split Delivery

301

Figure 1. Two products, A and B, are transported from 1 to 3. There are weight and volume constraints on arc (2,3). With a density of
2t/m3 , product A is much denser than product B with density 0.5t/m3 . Product B is more proﬁtable regarding weight and one might expect
to load as much as possible of product B. An optimal allocation, however, loads only 10t of B and 70t of A (see ﬁgure on right). [Color ﬁgure
can be viewed in the online issue, which is available at wileyonlinelibrary.com.]

3.3.

Level 1 RMP

Each variable (column) of the RMP 1 belongs to a ship v
and represents the two parts of a cargo route, the route itself
and a route speciﬁc cargo pattern. The set of ship v’s possible
routes is denoted by Rv . R∗v = Rv ∪ {0} comprises a dummyroute 0 which represents that no route is assigned to v. The
set of cargo patterns on route r ∈ Rv is denoted by Pvr . A
cargo pattern k ∈ Pvr consists of feasible quantities QPivrk
for pickup or delivery in time window i. Binary variable λvr
equals 1, if ship v sails route r (possibly the dummy route),
and 0 otherwise. If ship v sails route r, precisely the (extreme)cargo patterns k ∈ Pvr are convex combined to yield theactual pickup and delivery quantities on route r. Continuous
variables µvrk represent this convex combination.
The total cost of ship operation includes route cost and thecosts for picking up and delivering cargo. The constant Cvrk
denotes the cost for ship v on route r including service costs
caused by cargo pattern k. Let Mi be the minimum number ofvisits at a time window. Mi can be derived from the total time
window amount Qi and the ship capacities. Let AN
ivr equal
1 if ship v visits time window i on route r, and 0 otherwise.
The integer variable z reports the number of nonused ships.
The model can be written as follows:
 

min

Cvrk µvrk ,(1)

v∈V r∈Rv k∈Pvr

 

QPivrk µvrk = Qi∀i ∈ N ,(2)∀v ∈ V, r ∈ Rv ,(3)v∈V r∈Rv k∈Pvr

λvr −




µvrk = 0∀v ∈ V,v∈V(9)µvrk + λv0 = 1

 

AN
ivr µvrk ≥ Mi

(4)∀v ∈ V,(4 )∀i ∈ N .

(5 )v∈V r∈Rv k∈Pvr

ANivr λvr ≥ Mi∀i ∈ N ,(5)v∈V r∈Rv



(7)
(8)

r∈Rv k∈Pvr

r∈R∗v



∀v ∈ V, r ∈ Rv , k ∈ Pvr ,∀v ∈ V, r ∈ R∗v ,Objective function (1) minimizes the total cost of sailingand service. Constraint (2) ensures that in each time window
the sum of all pickups or deliveries made by all ships visiting the time window equals the total time window quantity.
The convex combination of (extreme) cargo patterns for eachroute hereby allows the selection of any pickup or delivery
amount feasible on the route. Constraint (3) couples feasible
cargo patterns for a particular route and the route itself: Cargopatterns can be convex combined only when their associatedroute is used. Each ship is allowed to sail at most one route,see constraint (4). Constraint (5) can be used to set a minimum number of visits for each time window. It is used tostrengthen the formulation. Constraints (6) and (9) serve thetechnical purpose to branch on the number z of unused ships(see Section 4.1). Constraints (7) and (8) enforce the correct
domains.
Each time a new route is added to the above formulation, aconstraint (3) needs to be added as well. This is particularly
unfortunate because one can expect the number of routes togrow rapidly. As branching on the binary master variables is
not advisable, we substitute for λvr in constraints (4) and (5)according to constraint (3), relax constraint (8), and drop (3).
Constraints (4) and (5) then take the following form:
 

k∈Pvr

λvr = 1µvrk ≥ 0
λvr ∈ {0, 1}z ∈ Z+ .

λv0 − z = 0(6)

Branching is performed on the original problem variables forarcs and nodes of the underlying network. See Section 4.1
for further discussion.
Naval Research Logistics DOI 10.1002/nav

302

Naval Research Logistics, Vol. 59 (2012)

3.4.



Level 2 RMP

Each RMP 2, one for each ship v, is a column generator
for RMP 1 that provides variables µvrk . Thus, its objectivefunction is the RMP 1 reduced cost of µvrk variables. EachRMP 2 accommodates binary routing variables λvr with thesame meaning as in RMP 1. The set Sv of feasible routes atlevel 2 is a superset of Rv because the feasible quantity split
is decided only at level 1. Set Sv∗ = Sv ∪ {0} again contains adummy route. Continuous variables qiv represent how muchweight is picked up or delivered by ship v in time window i. Aµ-column itself is derived from a selected optimal ship route(binary variable λvr ) and optimal pickup/delivery quantities(continuous variables qiv ).
Parameter δic = 1 signals that time window i supplies ordemands grade c. For all other grades δic = 0. Parameter Ii
indicates the type of time window i, where Ii = 1 means a
pickup and Ii = −1 means a delivery. Continuous variables
lij cv track the weight of grade c onboard ship v on arc (i, j ).
Parameter Aij vr equals 1 when route r for ship v containsarc (i, j ), and 0 otherwise. Continuous variables tiv represent
the start-of-service time of ship v at time window i. Only one
time variable per time window is needed because we assume
that a time window can be visited only once by a particular
ship (see Section 3.5).
To compute the level 1 reduced cost at RMP 2, dual information has to be transferred from the ﬁrst to the second
level. Relevant dual information concerns the route and cargopattern variables λv0 and µvrk :αi ∈ Rβv ∈ Rγi ∈ R+η∈R

for constraint (2),for constraint (4 ),for constraint (5 ), andfor constraint (6).

Dual variables βv and η are constant per ship v. The costEof route r for ship v is denoted by Cvr, and is computed fromarc costs and dual variables γi . The cost per weight unit forQwhichpickup and delivery in time window i is denoted by Civneeds to be corrected by the dual variable αi . An RMP 2 reads
min


r∈Sv

+

E
Cvr
λvr




Q
Civ
− αi qiv − βv − η,

(10)

i∈N

qiv ≤ Qi
lij cv − W ij v
Aij vr λvr ≤ 0∀i ∈ N ,(11)∀(i, j ) ∈ AW
v ,

(12)

Aij vr λvr ≤ 0 ∀(i, j ) ∈ AVv ,
D
c
c∈C
r∈Sv


lj icv + Ii δic qiv −lij cv = 0 ∀i ∈ N , c ∈ C,

(13)

c∈C

 lij cv

j ∈N

r∈Sv

− V ij vj ∈NNaval Research Logistics DOI 10.1002/nav(14)


Aij vr λvr tiv + TiQ qiv

r∈Sv∗


+ TijSv − tj v ≤ 0

T i ≤ tiv ≤ T i
λvr = 1∀(i, j ) ∈ A,

(15)

∀i ∈ Nv ,

(16)

,

(17)

r∈Sv∗

lij cv ≥ 0∀(i, j ) ∈ A, c ∈ C,(18)

tiv ≥ 0i ∈ Nv ,∀i ∈ N ,

(19)
(20)

∀r ∈ Sv∗ .

(21)

qiv ≥ 0λvr ∈ {0, 1}

Objective function (10) minimizes the total reduced cost.
Typically, a single ship services only a subset of time windows
and may service only a fraction of a time window quantity.
Constraint (11) therefore limits only the maximum pickup
and delivery amount in each time window. Constraints (12)and (13) limit the maximum available capacity on used arcs.
If an arc is not in use, the load onboard a ship on the arc
is kept at zero. Constraint (14) is a ﬂow conservation ofloads on arcs. Scheduling constraints (15) and (16) ensure
pickup and delivery amounts dependent on port stay times.
Constraint (15) allows a ship to start service at arc destination only after arrival at arc destination. Service times in anarc’s origin time window are respected in the arc traversal
time. Note that (15) is a nonlinear constraint which can be
re-formulated in a standard way using a “big M” mechanism:tiv + TiQ qiv + TijSv ≤ tj v + M(1 − r∈Sv∗ Aij vr λvr ), where M
is the largest possible duration between any possible arrivaland any possible departure of ship v at ports associated with iand j . Start-of-service has to take place within the time window speciﬁed in constraint (16). Waiting time between arrival
at a time window and time window opening cannot be limited by the above model. However, approximate limits can beapplied during route generation in the subproblem. At mostone route can be selected per ship, as ensured by constraint(17). Variable domain constraints (18) to (21) complete themodel.
3.5.

Subproblem at Level 3

At the lowest level is the column generation subproblem
for the RMP 2 problems. That is, the subproblem furnishes
routes, represented by sequences of time windows. It is arather standard elementary shortest path problem with time
windows (ESPPTW). A problem deﬁnition and discussion ofsolution algorithms based on dynamic programming can befound in [7] and in the survey [11]. The objective function is
the RMP 2 reduced cost for the λvr variables that are generated. It comprises real costs for ship operation and port visits,Hennig, Nygreen, and Lübbecke: Split Pickup and Split Delivery

and RMP 2 dual variables associated with constraints (12),(13), (15), and (17).
Due to the large distance between pickup and delivery
regions it is not time feasible for a single ship to visit a time
window on several occasions on consecutive voyages. Even
visiting a time window more than once during a single voyage
is considered unrealistic even though it could be beneﬁcial.
The so-called resource extension functions of the subproblem extend the routing cost and time, and information about
grades onboard. Cost is a nonmonotonic resource whereas
time is a strictly increasing resource. Time is extended basedon deterministic sailing times and smallest possible service
times required in real operations. Because pickup and delivery quantities determine the actual service time, the time
window arrival times calculated this way are only approximate. Because routes with extraordinary waiting times are notjudged reasonable in shipping, we enforce maximum waiting time limits for the duration between arrival at a port andstart of service (even though theoretically, this may discard
optimal routes). Again, the time limit can only be approximate. Grade information has to be tracked on partial paths.
An extension of a path to a delivery time window is infeasible
if a demanded grade is not available in already visited pickup
time windows. Similarly, a path cannot be extended from adelivery time window to a pickup time window, that is, to thebeginning of the next voyage, if for a grade picked up on thepresent voyage there is no matching delivery time window
on that voyage.
Label dominance is based on cost, time, visited time windows, and grade information. A label L1 containing all necessary information about a partial path can dominate a label L2
at the same node, if total cost of L1 is no larger than the costof L2 , start of service at the node for L1 is no later than for L2 ,the path of L1 visits a subset of nodes visited on the path ofL2 , and L1 and L2 have identical sets for pickup grades anddelivery grades. Barring grade compatibility, the describeddominance is the same as for the (E)SPPTW.
The main idea of dominance is that all extensions for label
L2 are also feasible for label L1 . In addition L1 may have more
possible extensions. Let us consider a single voyage: for theconsidered pickup and delivery problem with multiple products the visited pickup time windows in a pickup sequence
determine which delivery time windows are allowed to bevisited and which sort of delivery time windows have to bevisited. If, for example, a pickup time window supplying
grade A is visited, then all delivery time windows demanding
grade A are eligible for a visit. In fact, at least one delivery
time window demanding grade A has to be visited. If for
instance L1 visits a subset of pickup time windows visited by
L2 and can only pick up grades A and B in these time windows instead of A, B, and C available in L2 , extensions of
L1 visiting delivery time windows with the demanded grade
C are not allowed. Hence, L1 cannot be extended in the same303

way as L2 . If pickup grades are identical, but delivery grades
are not, L1 has to visit a certain sort of delivery time windowthat L2 already has visited. Also in this case dominance is
not possible.
We omit the rather standard mathematical formulation of
the subproblem and the technical description of the dynamic
program in this article. In our implementation, we use theresource constrained shortest path algorithm provided by the
boost C++ graph libraries (www.boost.org).

4.

COMPUTATIONAL ALGORITHM

The nested column generation algorithm has been implemented in C++ using the noncommercial branch-cut-andprice framework SCIP (scip.zib.de) developed by Achterberg [1]. SCIP can be used as a standalone mixed integer
programming solver as well and provides a widely extensible state-of-the-art tree management. Column generation is
directly supported via SCIP’s plugin-based architecture, andnested column generation is no problem either, even thoughseveral useful functionalities concerning branch-and-priceneed to be (and can be) added for the problem at hand. Even
though there may be a small overhead, for a proof-of-concept
SCIP is a very reasonable choice.
In Section 4.1, we describe our efforts to handle our typeof branching in SCIP. After explaining the initial algorithm
setup in Section 4.2, we describe the level 1 pricing strategy
in Section 4.3. For the latter, a pattern generation heuristic is
applied that is sketched in Section 4.4. Finally, Section 4.5
mentions a few actions taken to accelerate the algorithm.
4.1.

Branching

It is well known that branching (in either master problem) should not be done on the route variables directly
because the down-branch is ineffective and difﬁcult to handle in the subproblems. Moreover, this would lead to an
unbalanced branch-and-bound tree. Instead, we branch on
fractional ﬂows on arcs and nodes, the so-called “originalvariables,” calculated as follows:
xij v =Aij vrµvrk ∀v ∈ V, (i, j ) ∈ A, (22)
r∈Rv

yiv =



r∈Rv

AN
ivr

k∈Pvr



µvrk

∀v ∈ V, i ∈ N .

(23)

k∈Pvr

Parameter AN
ivr equals 1 if ship v on route r visits time window i, and 0 otherwise. We impose integrality on variables
xij v ∈ {0, 1} for arcs (i, j ) and yiv ∈ {0, 1} for nodes i for eachship. We can therefore directly use various variable selection
and branching strategies build-in in SCIP. These additional
constraints, together with their dual variables, are omitted inNaval Research Logistics DOI 10.1002/nav

304

Naval Research Logistics, Vol. 59 (2012)

Section 3.4 for reasons of clarity. Note that branching an xij v
or yiv to zero eliminates the corresponding already generated µvrk variables from the restricted master problem. The
inclusion of the additional dual variables ζij v for constraint(22) and ξiv for constraint (23) in the RMP 2 objective function is straight-forward. Branching in the RMP 2 is done
analogously. There, the arc and node constraints take thesame form as constraints (22) and (23) but λvr replaces termk∈Pvr µvrk .
4.2.

Algorithm Setup

The initial setup consists of an “empty” RMP 1, one“empty” RMP 2 for each ship and a structure for the dynamic
program (DP). Each route set R∗v only contains the dummyroute index 0. All sets Rv are empty. To ensure feasibility of the problem initially and during branch-and-price, we
introduce one artiﬁcial variable for each time window. The
coefﬁcients of these continuous variables equal to Qi fortime window i in constraint (2) and Mi for the minimum
visit constraint (5). These artiﬁcial variables are penalized in
the objective function so that no artiﬁcial variable appears in
an optimal solution. In all RMP 2, sets Sv and Sv∗ , respectively, are initially identical with RMP 1 sets Rv and R∗v .
The dynamic program is identical in structure for each ship.
Because the cost and feasible arc set for the dynamic program have to be updated each time an RMP 2 is solved, we
maintain only one dynamic program the parameters of which
are adjusted each time it is called.
In the following two sections, we explain the pricing
strategy and the actual generation of column for the ﬁrst level.

4.3.

a given ship’s RMP 2 and solves it by branch-and-price, identiﬁes, if desired, a set of cargo patterns, and extends the level
1 RMP with one or several new cargo patterns. Procedure
1 has as input two global (ﬁxed) parameters, t lim and clim .
For each RMP 2 we set an optimization time limit t lim . The
solution of a RMP 2 is stopped as soon as the time limit isreached or an improving column with reduced costs less thannegative threshold clim is found. We reach LP optimality of a
RMP 1 node when no improving column for any ship could
be identiﬁed, and all RMP 2 optimization times are less thant lim .

Level 1 Pricing Strategy

The RMP 1 pricing problem is to identify a minimum
reduced cost column, that is a least cost extreme cargo pattern as described in Section 3.1, with respect to certain dual
variable values. However, an often faster (heuristic) strategy
in branch-and-price algorithms is to just identify a numberof negative reduced cost columns instead of solving the subproblem to optimality. In our approach, this is even more
important because we experienced tail-off effects that canconsume large amounts of time when solving the level 2
column generation algorithm. Therefore, it seems natural
to terminate the level 2 optimization as soon as we ﬁnd
a few negative reduced cost solutions that are routes withcompatible cargo patterns.
The RMP 1 pricing strategy is described in Procedure 1,
which we call SolveNodeLevel1. It illustrates solving the LP
relaxation at a node in the RMP 1 branch-and-bound tree to
optimality. The actual pricing for a particular ship is sketchedin Procedure 2, called ExecutePricing. Procedure 2 initializes
Naval Research Logistics DOI 10.1002/nav

Procedure 1 SolveNodeLevel1(t lim , clim )1. q ← ∅2. for each v ∈ V do
3.
q.enqueue(v)4. end for5. success ← true
6. while success = true do
7.
GetBranching(x, y)
8. SolveRelaxation({Pvr : v ∈ V, r ∈ Rv }, x, y)
9.
p ← q // initialize ship queue with last iteration’s
result queue
10.
t←0
11.
while p = ∅ do12.
p.dequeue(v) // Consider ﬁrst ship in queue
13.
T ← max{t lim , tv }14.
(cmin , t) ← ExecutePricing(v, T , clim )15.
if cmin < 0 then16.
success ← true
17.
q.remove(v)18.
q.enqueue(v) // Queue ship in result queue
19.
break
20.
else
21.
success ← false
22.
if t < t lim then23.
q.remove(v)24.
q.enqueue(v) // Queue ship in result queue
25.
else
26.
p.enqueue(v) // Queue ship at the end ofcurrent queue
27.
tv ← ∞
28.
end if29.
end if30.
end while
31. end while

Procedure 1 is best explained by walking through an example: Given a set {A, B, C, D} of ships, Procedure 1, in lines
1 to 1, ﬁrst deﬁnes a processing order for each RMP 2 to besolved. This is done in forming a queue q = (A, B, C, D). At

Hennig, Nygreen, and Lübbecke: Split Pickup and Split Delivery

the beginning of each pricing iteration the linear relaxation
based on the current sets of columns Pvr for all ships androutes and the current branching situation for variables xij v
and yiv is solved (line 1). Columns representing routes areomitted as mentioned in Section 3.3. Each pricing iteration(line 1) may have several pricing rounds (line 1) depending
on how many level 2 RMPs are solved. A pricing round represents the execution of a level 2 RMP. The goal is to stopa pricing iteration as soon as a level 2 RMP delivers at least
one improving column (line 1) or if no improving columns
can be found for any ship.
The ship queue to be processed in the current pricing iteration is called p and the ship queue for the following iteration
is called q. Queue p is initialized with q coming from theprevious iteration or the initial deﬁnition described above.
The size of p at processing start equals the number of ships|V|, which is also the principle limit of pricing rounds ofthe ship processing. The components tv of vector t represent
ship speciﬁc level 2 optimization time limits. These time limits are initialized with 0. In the ﬁrst pricing round ship v = A
is dequeued from p reducing queue p from (A, B, C, D) to(B, C, D). Queue q = (A, B, C, D) remains unchanged. The
actual RMP 2 optimization limit T is set to t lim . Then pricing is executed for ship A through procedure ExecutePricing.
The pricing procedure returns the minimum reduced cost cminfound by the second level and the actual level 2 optimization
time t.
If cmin is negative, new improving columns have beenfound and the level 1 RMP has been extended by Procedure 2.
Ship A is removed from next iteration’s queue q and enqueued
at the end of the queue: q = (B, C, D, A). The pricing iteration is completed and the whole process starts over again
with LP relaxation solving. In the case where no improving
columns for A exist, that means cmin = 0 and t < t lim , ship
A is also queued at the end of q but a new pricing round has
to be started with the next ship in current queue p, namely B.
Let us assume the latter case and ExecuteP ricing is executed for ship B. At this point p equals (C, D) and q stillequals (B, C, D, A). Now assume that cmin takes on value 0
and the time limit is reached. In this case no conclusion canbe drawn and ship B is enqueued in p = (C, D, B). This way
B may be examined again but next time without the optimization time limit. In the worst case pricing rounds for C and Dlead to no conclusion either. The current queue would thenbe p = (B, C, D) and B would have to be examined again.
Of course, without the optimization limit a clear conclusion
can be drawn after each round.
Procedure 2 facilitates the level 2 branch-and-price. Before
the level 2 problem can be solved, procedure FixLevel2Variables (line 2) ﬁxes all RMP 2 branching variables xij v
for arcs and yiv for nodes already ﬁxed in the RMP 1 to 0 or
1. Procedure UpdateLevel2Objective (line 2) updates the dual
information αi , βv and η in the RMP 2 objective and changes

305

Procedure 2 ExecutePricing(v, T , clim )1. GetBranching(x v , yv )2. GetDuals(α, βv , γ , η, ζ v , ξ v ) // Get dual information
3. FixLevel2Variables(x v , yv )4. UpdateLevel2Objective(α, βv , γ , η, ζ v , ξ v )5. UpdateDPCost(γ , ζ v , ξ v )6. (cmin , t, C) ← SolveLevel2(clim , T )7. if cmin < 0 then
8. for each r ∈ C do
9.
if r ∈/ Rv then10.
Rv ← Rv ∪ {r}11.
Create(Pvr ← PatternGeneration(r))
12.
else
13.
Pvr ← Pvr ∪ GetPatterns(r)14.
end if15.
ExtendLevel1RMP
16. end for
17. end if
return cminE
the route costs Cvr
for all columns already contained in Sv∗
based on the new dual information ζij v on arcs, and ξiv andγi on the nodes. Finally, the same dual information is used byprocedure UpdateDPCost (line 2) to change the arc costs inthe dynamic program. After this initialization step procedure
SolveLevel2 executes the level 2 branch-and-price in line 2.
Procedure SolveLevel2 is interrupted as soon as the ﬁrst negative RMP 1 reduced cost column better than clim is foundor when the time limit is reached. At the outset of Procedure SolveLevel2 a quick test is done to check if any column
already known from previous pricing iterations for the considered ship has a reduced cost smaller than clim . If that is thecase, branch-and-price is not executed. The actual optimization time is recorded in t. This t does not have any impact on
Procedure 2 and is only transferred to Procedure 1 (see line 1
in Procedure 1). All routes from feasible solutions yielding anegative reduced cost column are recorded in set C. Because
SolveLevel2 is a column generation procedure, route sets Sv
and Sv∗ grow during optimization. In general, both sets growduring each pricing round, and of course all routes found inprevious iterations in all predecessor tree nodes are kept andneed not be generated again by the subproblem.
The value clim is negative only if improving RMP 1
columns have been identiﬁed. A level 2 solution returns theroute used and a corresponding feasible cargo pattern, thatcan be used as RMP 1 column. Set C may include one orseveral different routes. For a single route, negative reduced
cost solutions may exist with one or several different cargopatterns. To potentially reduce the number of RMP 1 pricing
iterations we use a heuristic cargo pattern generator, PatternGeneration, to identify, for each route found, several feasible

Naval Research Logistics DOI 10.1002/nav

306

Naval Research Logistics, Vol. 59 (2012)cargo patterns at once. If a route is not yet contained in the
RMP 1 route set Rv we include the route (line 2) and create
a new cargo pattern set Pvr ﬁlled with a number of feasible cargo patterns (line 2). It may happen that the resulting
columns may not all have negative reduced cost but together
span a polyhedron that describes feasible cargo patterns forthe given route. In the best case this region will be sufﬁciently
large to provide all pickup and delivery quantities needed
during level 1 optimization. Should the region be too small,
C will at a future iteration, comprise a route already contained in Rv . To avoid an accidental generation of identicalcolumns we do not execute PatternGeneration in this case andonly include all additional cargo patterns that yielded negative
reduced cost solutions (line 2). Finally, Procedure 2 extends
the level 1 RMP and returns the minimum reduced cost
identiﬁed.
The complete level 2 branch-and-price is solved in SCIP
with default settings. The dynamic program is called in eachtree node as long as negative reduced cost columns can be
identiﬁed. All columns with negative reduced costs provided
by the dynamic program in one pricing iteration are included
in sets Sv and Sv∗ .
4.4.

The cargo pattern generator takes into account all constraints included in level 2. Once a cargo quantity for a pairis decided, time window feasibility along the route must be
ensured. Otherwise a cargo pattern would not be feasible. To
be able to do this test for the second voyage, we assume maximum possible pickup and delivery times on the ﬁrst voyage
and thus can derive a latest start of voyage time for the second
voyage.
The number of cargo patterns generated this way may belarge. A small example may clarify this. A single voyage with
only two pickup and two delivery time windows supplying
and demanding the same grade has four time window origindestination pairs. Hence there are 4! = 24 possible sequenceswith at most 24 different cargo patterns. If a route contains
two such voyages already 24 times 24 cargo patterns are possible. Clearly, many of the single voyage cargo patterns may
be redundant, but especially for longer routes with few grades
the number of unique cargo patterns grows substantially. To
avoid unreasonably large matrices of the linear programs tobe solved in each node of the branching tree, we allow to
specify a maximum number of cargo patterns to be included
into the RMP 1 in one pricing iteration.

Cargo Pattern Generation

Given a route, it is simple to ﬁnd a number of different cargo patterns heuristically. These are not necessarily
extreme but hopefully provide a reasonable approximation.
If needed, extreme cargo patterns are generated if the heuristic fails. We apply a strategy that lends itself from the lesscomplicated subproblem described in [6]. We assume that all
routes contain at most two voyages, which is consistent withthe realistic data available. The procedure works as follows:
1. Find all grade feasible pairs of one pickup and one
delivery time window for each voyage in the known
route.
2. Deﬁne all sequences of pairs for each voyage, wherea sequence is a unique ordering of pairs.
3. For each sequence greedily service as much cargo as
possible for the ﬁrst, then second, then third, etc. pairin the sequence.
4. Combine all non-identical cargo patterns found thatway for the voyages in all possible ways.
This procedure ignores the fact that it may not be optimal
to service as much cargo as possible for a given pair as
explained in Section 3.1. It also ignores that for different timesof beginning the second voyage, different cargo patterns may
be obtained. This enables us to generate cargo patterns foreach voyage independently and combine these patterns afterwards. An empty cargo pattern is also assigned to each route,which represents pickup and delivery quantities equal to 0 ineach time window.
Naval Research Logistics DOI 10.1002/nav

4.5.

Acceleration

SCIP provides a number of heuristics to identify feasible
solutions used at different stages of the algorithm. In addition
to these default heuristics we include two simple heuristics,
one at each level: On level 1 we consider the LP relaxation
solution at each pricing iteration and try to ﬁnd a feasible
solution in the following way: Identify which (fractionally)used route of which ship is used the most. Fix all columns ofthat ship not corresponding to that route to 0 and reoptimize
the LP problem. Continue this until one route for each ship
has been selected and check if the resulting solution is feasible with respect to the total quantity to transport. On level 2
we know that each route in route set S and each route provided by the dynamic program yields a feasible solution in
the RMP 2. However, these solutions may not be identiﬁed
during branch-and-price. Therefore we apply a fast probing
procedure provided by SCIP to solve the LP relaxation witheach route in S at the beginning of Procedure 2 and eachroute provided by the dynamic program one at a time ﬁxed
to 1.
A further acceleration can be achieved by solving the
dynamic program with a heuristic dominance rule. Whenever the heuristic dominance rule does not yield negative
reduced cost columns, the correct dominance rule can beapplied. The heuristic dominance rule we apply only takes
into account costs and timing when comparing two partial
paths. The information about grades picked up and deliveredis ignored.

Hennig, Nygreen, and Lübbecke: Split Pickup and Split Delivery

5.

Table 2.

COMPUTATIONAL RESULTS
5.1.

Instances and Settings

We use the same reality related test instances as in Hennig
et al. (working paper). The test instances are grouped into two
sets. The smaller instances described in Table 1 are of varying
complexity. The time periods of the ﬁrst three instances only
allow routes with single voyages. In contrast, instances 4 to
6 require routes with up to two voyages. Table 2 illustrates
the larger instances. Each instance, starting from instance 8,
has a few additional time windows and ships as compared toits previous instance.
We run SCIP 1.1.0.6 with CPLEX 11.2 as LP solver on
64-bit PCs with Intel Core 2 Duo CPUs, 3 GHz, 6 gigabyte
RAM (only one core is used per instance) and linux operating system. Two instances are run simultaneously on thesame computer as long as there is sufﬁcient memory. We let
the branch-and-price run for 12 hours. Because of the timingof the running time check points, longer branch-and-pricetimes are possible. However, this does not have any signiﬁcant impact on the solution. During branch-and-price we haveto disable SCIP’s preprocessing, cut generation, and propagation functionalities to ensure a correct dual bound. For t lim
we choose 60 seconds and clim = −1. We allow at most 30
cargo patterns per voyage to be selected in the pattern generator. After 12 hours we solve the resulting level 1 masterproblem without pricing with default CPLEX. In this second
run we solve the original problem comprising objective function (1) and constraints (2) to (9) with all routes and cargopatterns found during the ﬁrst run. Preprocessing and cut generation can be applied. Branching is performed on variables
λvr and z.
5.2.

Results and Discussion

The computational results are shown in Tables 3 and 4.
We assess the results with respect to solution quality, algorithm performance, and capability of the underlying model.
In the result tables, we can distinguish between results basedon the ﬁrst branch-and-price run carried out with SCIP and
the second run without pricing performed with CPLEX. The
ﬁrst part of each table, above the broken line, is dedicatedto SCIP. The parts under the broken lines report the CPLEX
information.
Table 1.

Instance set A.

Instances

1

2

3

4

5

6

No. of ships
No. of products
No. of loading time windows
No. of discharging time windows

2
2
2
5

3
2
3
8

5
5
6
14

2
4
5
8

4
5
6
12

6
5
6
16

307
Instance set B.

Instances

7

8

9

10

11

12

No. ships
No. products
No. loading time windows
No. discharging time windows

8
4
9
9

8
5
12
13

9
6
14
15

10
7
15
18

12
8
18
19

15
11
21
26

Before we start the discussion, the following explanations
may help to read the result tables: For the ﬁrst run the best
lower bound (at optimization time limit) and the objective
value at root node LP relaxation optimum are reported. The
root node gap is the relative gap between the best solution
and the root node LP relaxation. This gap gives a measurefor the tightness of the LP relaxation. In fact, best solution
refers to the best solution identiﬁed by CPLEX. Therefore we
do not report optimality gaps between the best lower bounds
and the best solutions, because best solutions have in mostcases not been identiﬁed in the ﬁrst run. Pricer calls takeplace whenever new dual information is available during thesolving of a node in the branching trees. The number of calls
is reported for the ﬁrst level RMP and for all second level
RMPs in total. Route regenerations refer to the number oftimes level 2 provides a route that is already contained in
the ﬁrst level RMP. The time needed to ﬁnd all columns contained in the best solution is reported under time to best cargopattern set. Under the broken lines in the tables the numberof explored nodes and the total running time for the second
run are reported. Except for instances 1, 2, and 4 the time tobest solution is the sum of the time used for branch-and-priceand the time until CPLEX ﬁnds the best solution.
Let us ﬁrst take a look at the reported (best) solutions: In
all instances except instance 6 and 12 the entire available ﬂeet
has been used. Optimality is only proven for instances 1, 2,
and 4. All other instances may not be optimal. Optimality
gaps are not reported for two reasons. First, for instances 8
to 12 RMP 1 root node LP relaxations have not been solved
to optimality in the given time limit and thus no valid lower
bounds are provided. Second, only instances 1, 2 and 4 provide interesting feasible solutions during the ﬁrst run. For allother instances valuable solutions are ﬁrst identiﬁed in thesecond run. However, where applicable we provide root node
gaps that are essentially identical to the optimality gaps. The
average gap for the multiple voyage instances 4 to 7 is only
4.9%, which is quite small for the freedom provided by themodel. It is conspicuous that the relaxation for instances 1 to
3 allowing only one voyage per route is much tighter than for
the instances with longer time periods. However, a general
conclusion should not be drawn. The quality of the results
greatly depends on the routes and cargo patterns generated
for the ﬁrst level RMP. Hence, the performance of the algorithm especially on the second and subproblem level plays
an important role.
Naval Research Logistics DOI 10.1002/nav

308

Naval Research Logistics, Vol. 59 (2012)
Table 3.

Instances
No. of ships used in best solution
Best solution
Best lower boundRoot node LP relaxation% root node gap
No. of RMP 1 pricer calls
No. of RMP 2 pricer calls
No. of RMP 1 cargo patterns
No. of RMP 1 routes
No. of route regenerations
Time to best cargo pattern set (s)Total branch-and-price time (s)% of B&P time used by level 2% of RMP 2 time used by DP
No. of explored nodes on level 1
No. of explored nodes on level 2
No. of explored nodes in ﬁnal branching
Additional branching time (s)Time to best solution (s)Instance set A.

1

2

3

4

5

6

2
2681
2681
2670
0.4
48
614
540
44
0
1
2
100
4.7
22
506
N/A
N/A
2

3
4112
4112
4088
0.6
439
19,591
12,698
449
0
143
496
96
3.0
161
17,325
N /A
N /A
156

5
6936
6831
6831
1.5
802
13,23929,1221124
60
30,73343,260
100
3.2
2
8499
232,583
10,036
46,915

2
5121
5121
4804
6.6
913
33,63354,399
654
89
2499
4262
96
2.7
392
31,493
N /A
N /A
3550

4
8102
7702
7702
5.2
586
26,97544,745640
58
685
43,200100
1.8
3534,86931,9252201
45,390

4
7587
7416
7408
2.4
254
28,711
21,469345
15
38,17743,380100
3.5
1122,335435
1943,398Performance can be judged by the number of routes andpatterns generated and the time this takes. It turns out that the
number of RMP 1 pricer calls clearly is limited by the timeneeded on level 2. Basically, the entire optimization time isconsumed by the second level. Pricing iterations there are
signiﬁcantly faster. Compared to level 1, many pricer calls
are executed, that means many calls to the dynamic program
are made. However, at most 5.8% of the level 2 optimization time is consumed by the dynamic program. The sumof the RMP 1 routes and of the route regenerations for all
instances exceeds the number of RMP 1 pricer calls. That

Table 4.
Instances
No. of ships used in best solution
Best solution
Best lower boundRoot node LP relaxation% root node gap
No. of RMP 1 pricer calls
No. of RMP 2 pricer calls
No. of RMP 1 cargo patterns
No. of RMP 1 routes
No. of route regenerations
Time to best cargo pattern set (s)Total branch-and-price time (s)% of B&P time used by level 2% of RMP 2 time used by DP
No. of explored nodes on level 1
No. of explored nodes on level 2
No. of explored nodes in ﬁnal branching
Additional branching time (s)Time to best solution (s)Naval Research Logistics DOI 10.1002/nav

means, that with each RMP 1 pricing iteration at least one
route with most often many cargo patterns is included into
the RMP 1. A route is included several times if it is beneﬁcial
to extend the existing polyhedron of cargo patterns for thatroute. The numbers show clearly that there are many uniquecargo patterns for the identiﬁed routes. Cargo pattern generation, which is part of the ﬁrst level, consumes negligible
time. With at most 30 cargo patterns per single voyage or
900 patterns per two-voyage-route we allow the RMP 1 togrow fast with each new route. Still route regeneration takes
place especially for the larger instances and shows that cargo

Instance set B.

7

8

9

10

11

12

8
8471
8056
8041
5.2
975
71,497
12,387888
148
7812
43,200100
1.3
147
59,540817,400
5408
3835

8
10,405
N /A
—
N /A
259
20,7096893
289
6938,12443,260100
2.5
115,5289834
73
43,3209
11,471
N /A
—
N /A
345
13,7889248
340
99
39,91443,320100
2.2
110,32378,290797
44,072

1012,927
N /A
—
N /A
259
4135
7647
299
55
15,082
43,200
100
5.8
1
1807
30,449305
43,485

12
14,523
N /A
—
N /A
294
1686
9155
309
79
24,39847,100
100
3.6
1
761
2,317,267
17,23750,53914
16,998
N/A
—
N/A
360
2302
7921
380
87
35,14843,260
99.6
1.2
1
962
1,028,948
10,96853,385

Hennig, Nygreen, and Lübbecke: Split Pickup and Split Delivery

patterns are not easily obtained. The fractions of total timeused by level 1 underline that producing many patterns poses
no problem with LP relaxation solving.
The tree sizes of the second level reﬂect the complexity ofnested column generation. We try to reduce time consuming
tail-off problems by stopping optimization soon after the ﬁrst
improving solutions have been found. Still we face anothertime problem. If we consider the number of times level 2 pricing has been executed and the number of processed nodes onlevel 2 in total, we see that on average there are only few pricer
calls in each tree node. Subproblem solving consumes comparable little time and so time consumption is mainly causedby solving level 2 LP relaxations. We observe a few cases ineach of the largest instances where solving a single RMP 2
takes more than 1000 seconds. If we could disregard thesecases as problematic exceptions we still would have average
single node LP solving times for the largest instances of up
to 42 seconds.
The second run underlines the complexity of the largest
instances. Although all columns are known, and presolving
and cut generation are applied, CPLEX needs ∼3 to 4 hours tosolve instances 11 and 12. In fact, instances 8 and 12 runningalone on a single computer could not be solved to optimality
due to insufﬁcient memory. The optimizations terminate with
1.6% and 3.2% optimality gap.
There are two key challenges to the proposed algorithm:
time consumption on level 2 and the ability to ﬁnd goodsolutions during the ﬁrst run. A technical way to improve
running times would be to solve all RMP 2 for one pricing
iteration in parallel. This would eliminate the need for parameter tlim and a large speed up may be obtained because often
one of the RMP 2 solves relatively quickly, while others mayneed exceptional running times. A more elaborate patterngeneration may also reduce running time through reducing
the number of route regenerations. At least it is not advisable to replace the proposed pattern generation strategy withsolving an LP for a given route for different sets of dual variable values or omit the pattern generation completely. Bothapproaches are signiﬁcantly slower than the proposed one. It
may also be beneﬁcial to reduce the RMP 1 and RMP 2 matrices through omitting arc and node constraints and implicitly
branching on arcs and nodes through branching on sets ofcolumns. This however, would have prevented us from using
advanced branching techniques and heuristics already implemented in SCIP. Good solutions have to be found heuristically. The employed RMP 1 heuristic mentioned in Section
4.5 helps to ﬁnd solutions but do not work well enough to
replace the second run. A few instances show a considerably
shorter time to best cargo pattern set than to algorithm termination as is the case in many branch-and-bound contexts.
This should be exploited.
The underlying model provides great freedom. We have
limited the routing possibilities only with respect to realistic309

constraints and not to achieve a faster algorithm. If a larger
number of potentially visited time windows on a voyage
is desired, the algorithm would remain the same. Arbitrary
splitting of time window quantities is permitted. In (Hennig
et al., working paper) the same test instances were run withtwo column pregeneration algorithms, one allowing arbitrary
splitting of cargo quantities and the other one based on discretized splits. With the ﬁrst algorithm no feasible solution
was found for many of the larger test instances. For all but oneof the other instances we now improve results slightly and in
one case up to 8%. We improved the results compared to thesecond algorithm in (Hennig et al., working paper) between
1 and 7%. However, we did not ﬁnd better solutions for thetwo largest instances. In both instances we now use one additional ship which leads to a cost increase of 5% for instance
12. Interestingly, costs for instance 11 do not increase due
to one voyage less in the new solution. Note that the costs
reported in Tables 3 and 4 are not directly comparable to
the results reported in Hennig et al. (working paper). The
above discussion is based on a postoptimization step similar
to the one described in Hennig et al. (working paper). It has tobe pointed out that the proposed branch-and-price approachconsiders a signiﬁcantly larger feasible region than the othermentioned models. The same feasible region would not havebeen possible to consider in a column pregeneration context.

6.

CONCLUSIONS AND FURTHER RESEARCH

In this article, we present a rich routing and scheduling problem eminent in today’s industry. Unlike classicalassumptions in the vehicle routing literature, we consider aheterogenous vehicle/ship ﬂeet, multiple commodities, andtime window service quantity dependent service times. The
problem includes many-to-many relations for pickup anddelivery of single commodities and time window sequence
dependent vehicle capacities. We approach this problem witha nested column generation algorithm decomposing the problem into many single ship problems which are themselves
decomposed into a route generation and a cargo decision
problem. This way we avoid a dynamic program, that is possibly extremely difﬁcult to solve. Our results demonstrate thecapabilities of the approach.
For the SDVRP, the largest instances optimally solvable
today comprise about 100 customers [6]. Allowing a split ofloads in both the pickup and delivery time windows implies acomplication of this already tremendously difﬁcult optimization problem. The problem gets even more complex with arc
dependent capacity constraints and quantity dependent time
window service times. As a computational consequence this
implies a further combinatorial explosion of solution possibilities, so it is not surprising that the resulting type of SPSDP
is even harder to solve. It is precisely this explosion (in ﬁnding
Naval Research Logistics DOI 10.1002/nav

310

Naval Research Logistics, Vol. 59 (2012)

optimal combinations of split loads) which consumes most of
our reported computation times. Thus, being able to provide
high quality solutions to instances with almost 50 pickup anddelivery locations (as we do) is noteworthy.
Our research immediately points to further research directions. Even though we have a satisfactory proof-of-concept,
there is still room for improvement. We did not, for instance,
add valid inequalities to any restricted master problem inorder to strengthen the dual bounds, nor did we implement
tailored primal heuristics to be used in our nested decomposition approach. More elaborate branching rules for the
SPSDP could help solving the RMP 2 level quicker. It would
be interesting to study the SPSDP from a generic point ofview. A further very interesting subject of study should be toﬁnd efﬁcient algorithms to solve the cargo pattern generation
problem quickly and in an exact fashion. Efﬁcient branchingrules for the SPSDP are also to be developed.
ACKNOWLEDGMENTS
Most of this work was done while the ﬁrst and third authors
were with the Combinatorial Optimization & Graph Algorithms Group at Technische Universität Berlin, Germany.
The authors thank Dr. Michael Drexl for his comments andsupport in connection with the use of the boost resource
constraint shortest path algorithm.
REFERENCES[1] T. Achterberg, Scip: Solving constraint integer programs,
Math Program Comput 1 (2009), 1–41.
[2] C. Barnhart, E.L. Johnson, G.L. Nemhauser, M.W.P. Savelsbergh, and P.H. Vance, Branch-and-price: Column generation for solving huge integer programs, Oper Res 46 (1998),
316–329.
[3] J.M. Belenguer, M.C. Martinez, and E. Mota, A lower boundfor the split delivery vehicle routing problem, Oper Res 18(2000), 801–810.
[4] G. Bronmø, B. Nygreen, and J. Lysgaard, Column generation
approaches to ship scheduling with ﬂexible cargo sizes, Eur J
Oper Res 200 (2010), 139–150.

Naval Research Logistics DOI 10.1002/nav[5] A. Ceselli, G. Righini, and M. Salani, Column generation forthe split delivery vehicle routing problem, Technical report,
University of Milan - DTI - Note del Polo n. 118, 2009.
[6] G. Desaulniers, Branch-and-price-and-cut for the split delivery vehicle routing problem with time windows, Oper Res 58(2010), 179–192.
[7] J. Desrosiers, Y. Dumas, M.M. Solomon, and F. Soumis,“Time constrained routing and scheduling,” in: M. Ball, T.L.
Magnanti, C.L. Monma, and G.L. Nemhauser (Editors), Network routing, Volume 8, Handbooks in operations research
and management science, Chapter 2, Elsevier Science B. V.,
North-Holland, Amsterdam, 1995, pp. 35–139.
[8] J. Desrosiers and M.E. Lübbecke, “A primer in column generation,” in: G. Desaulniers, J. Desrosiers, and M.M. Solomon(Editors), Column generation, Springer-Verlag, Berlin, 2005,
pp. 1–32.
[9] M. Dror, G. Laporte, and P. Trudeau, Vehicle routing with split
deliveries, Discr Appl Math 50 (1994), 239–254.
[10] F. Hennig, B. Nygreen, M. Christiansen, K. Fagerholt, K.C.
Furman, J. Song, G.R. Kocis, and P.H. Warrick, Maritime
crude oil transportation—A split pickup and split delivery
problem, Eur J Oper Res 218 (2012), 764–774.
[11] S. Irnich and G. Desaulniers, “Shortest path problems withressource constraints,” in: G. Desaulniers, J. Desrosiers,
and M.M. Solomon (Editors), Column generation, SpringerVerlag, Berlin, 2005, pp. 33–65.
[12] M. Jin, K. Liu, and R.O. Bowden, A two-stage algorithm withvalid inequalities for the split delivery vehicle routing problem,
Int J Prod Econ 105 (2007), 228–242.
[13] C. Lee, M.A. Epelman, C.C. White III, and Y.A. Bozer, A
shortest path approach to the multiple-vehicle routing problem
with split pick-ups, Trans Res Part B 40 (2006), 265–284.
[14] M.D. McKay and H.O. Hartley, Computerized scheduling ofseagoing tankers, Nav Res Logist 21 (1974), 255–264.
[15] R. Nersesian, The economy of shipping venezuelan crude tochina, Explor Prod: Oil Gas Rev 2 (2005), 78–80.
[16] M. Nowak, Ö. Ergun, and C.C. White, Pickup and delivery
with split loads, Trans Sci 42 (2008), 32–43.
[17] UNCTAD. Review of maritime transport 2011, United
Nations, New York and Geneva, 2011.
[18] I. Vacca and M. Salani, The vehicle routing problem with discrete split delivery and time windows, In Proceedings of the
9th STRC Swiss Transport Research Conference, September
9, 2009 in Monte Verita Switzerland, 2009.
[19] F. Vanderbeck, A nested decomposition approach to a threestage, two-dimensional cutting-stock problem, Manage Sci 47(2001), 864–879.