Log In Sign Up

Industrial and Tramp Ship Routing Problems: Closing the Gap for Real-Scale Instances

by   Gabriel Homsi, et al.

In a recent study, Hemmati et al. (2014) proposed a class of ship routing problems and a benchmark suite based on real shipping segments, considering pickups and deliveries, cargoes selections, ship-dependent starting locations, travel times and costs, time windows, incompatibility constraints, among other features. Together, these characteristics pose considerable challenges for exact and heuristic methods, and some cases with as few as 18 cargoes remain unsolved. We propose an exact branch-and-price (B&P) algorithm and a hybrid metaheuristic. Our exact method generates elementary routes, but exploits decremental state-space relaxation to speed up column generation, heuristic strong branching, as well as advanced preprocessing and route enumeration techniques. Our metaheuristic is a sophisticated extension of the unified hybrid genetic search. It exploits a set-partitioning phase and uses problem-tailored variation operators to efficiently handle all the problem characteristics. As shown in our experimental analyses, the B&P optimally solves 239/240 instances within one hour, with up to 50 ships, 130 cargoes and therefore 260 pickups and deliveries. The hybrid metaheuristic outperforms all previous heuristics and produces near-optimal solutions within minutes. These results are noteworthy, since the largest instances of this set are comparable in size with the largest problems routinely solved by shipping companies.


page 1

page 2

page 3

page 4


The Electric Vehicle Routing Problem with Nonlinear Charging Functions

This paper outlines an exact and a heuristic algorithm for the electric ...

Arc Routing with Time-Dependent Travel Times and Paths

Vehicle routing algorithms usually reformulate the road network into a c...

Vehicle Routing with Time-Dependent Travel Times: Theory, Practice, and Benchmarks

We develop theoretical foundations and practical algorithms for vehicle ...

Learning to Solve Soft-Constrained Vehicle Routing Problems with Lagrangian Relaxation

Vehicle Routing Problems (VRPs) in real-world applications often come wi...

Solving the Joint Order Batching and Picker Routing Problem for Large Instances

In this work, we investigate the problem of order batching and picker ro...

A hybrid primal heuristic for Robust Multiperiod Network Design

We investigate the Robust Multiperiod Network Design Problem, a generali...

1 Introduction

International trade depends heavily on ship transportation, as it is the only cost-effective means for the transportation of large volumes over long distances. It is common to distinguish between three main modes of operation in maritime transportation: liner, industrial, and tramp shipping. Liner shipping, which includes container shipping, is similar to a bus service: fixed schedules and itineraries must be followed. In industrial shipping, the operator owns the cargoes and controls the fleet, trying to minimize the cargo transportation cost. Finally, a tramp shipping operator follows the availability of cargoes in the market, often transporting a mix of mandatory and optional cargoes with the goal of maximizing profit.

In this work, we focus on a class of ITSRSP. This class of problems typically arises in the shipping of bulk products such as crude oil; chemicals and oil products (wet bulk); and iron ore, grain, coal, bauxite, alumina, and phosphate rock (dry bulk). In , these product types constituted more than of the weight transported in international seaborne trade. Yet, in the wake of the financial crisis in 2008, the freight rates in the dry bulk shipping segment dropped dramatically: the Baltic Dry Index dropped more than , has remained low since then, and experienced record lows in 2016. Additionally, for the fifth year in a row, world fleet growth has been decelerating. Despite this decline, the supply of shipping capacity increased faster than demand, leading to a continuing situation of shipping overcapacity and downward pressure on freight rates (UNCTAD 2017). In this environment, a shipping company can be profitable only if its fleet is routed effectively.

In the ITSRSP, a shipping company has a mix of mandatory and optional cargoes for transportation. Each cargo in the planning period must be picked up at its loading port within a specified time window, transported, and then delivered at its destination port, also within a given time window. The shipping company controls a heterogeneous fleet of ships; each ship has a given initial position and time for when it becomes available for new transportation tasks. Compatibility constraints may restrict which cargoes a ship can transport (for example, due to draft limits in the ports). The shipping company may charter ships from the spot market to transport some of the cargoes. The planning objective in the ITSRSP is to construct routes and schedules, deciding which spot cargoes to transport and which cargoes will be transported by a spot charter, so that all mandatory cargoes are transported while maximizing profit or minimizing costs. The ITSRSP extends the PDPTW with a heterogeneous fleet, compatibility constraints, different ship starting points and starting times, and service flexibility with penalties. The interplay of these complex attributes requires the joint optimization of multiple decision sets.

In a recent work from Hemmati et al. (2014), a set of benchmark instances based on real shipping segments with seven to 130 cargoes (pickup-and-delivery pairs) has been made available to the academic community. The authors also presented a compact mathematical formulation, and solved it with a branch-and-cut algorithm to obtain initial results. However, some instances with as few as 18 cargoes remain unsolved. Clearly, given the current scale of industrial applications, a significant methodological gap must be bridged to respond to practical needs. To compensate for the lack of exact solutions, Hemmati et al. (2014) designed an ALNS heuristic, and subsequently investigated the impact of randomization as well as that of various search operators (Hemmati and Hvattum 2016). However, due to the lack of available lower bounds or optimal solutions, the true performance of these methods is unknown for large problems.

This paper contributes to fill this methodological gap, from an exact and heuristic standpoint.

  • [leftmargin=*]

  • Firstly, we introduce an efficient BnP algorithm for the ITSRSP. It relies on the generation of elementary routes, but exploits DSSR and extensive preprocessing to speed up labeling and pricing, as well as strong branching. Efficient correction strategies allow to maintain the DTI, which can become invalid due to the dual costs but is fundamental for dominance. The BnP is then extended using route enumeration (possible thanks to a sophisticated sequence of completion bounds), inspection pricing, and separation of subset-row cuts. This is first time these methodological building blocks are adapted, improved and combined into a outstanding algorithm for this problem class.

  • Secondly, to quickly generate high-quality solutions, we introduce a HGS. Our approach follows the same principles as the UHGS of Vidal et al. (2014). Yet, UHGS was never applied to heterogeneous fixed fleet and pickup-and-delivery problems as these problems require structurally-different local-search neighborhoods and variation operators (e.g., crossover) to be efficiently handled. Built on a completely new code base, our algorithm bridges these gaps. It uses problem-tailored crossover, LS operators and vehicle-dependent neighborhood restriction strategies to efficiently optimize all aspects of the ITSRSP and take into account its numerous constraints. It is also complemented by a set partitioning intensification procedure so as to stimulate a faster convergence towards good route combinations.

  • We report extensive experimental analyses on the industrial instances from Hemmati et al. (2014) to measure the performance of the new methods. The BnP algorithm is able to solve out of the available instances to optimality within a time limit of hour. The last instance is solved in hours and minutes. Moreover, our HGS finds very accurate solutions within a few minutes, with an average gap of relative to the optima. The quality of these solutions is far beyond that of the previously existing ALNS. These results are remarkable because the instances of Hemmati et al. (2014) were built to withstand the test of time. The ability to solve all of them within five years reflects the considerable progress made by exact methods for rich routing applications. The largest solved instance has 50 ships and 130 cargoes, and therefore 260 pickups and deliveries overall. This size is comparable to the largest problems solved routinely by shipping companies: e.g., Wilson, which is among the largest bulk shipping companies, operates 117 bulk ships between 1.500 and 8.500 deadweight tons and divides its operations into three main segments with 40 ships each (Wilson 2018). We therefore reach a turning point where state-of-the-art exact methods become sufficient for daily maritime practice.

2 Problem Statement and Related Literature

The ITSRSP is defined on a complete digraph , where is the union of a set of pickup nodes , delivery nodes , and starting locations . A tramp or industrial shipping operator owns a fleet of ships , and cargoes are available for transportation. Each cargo is characterized by a load and must be transported from a pickup to a corresponding delivery location . Therefore, for , and . Every node is associated with a time window of allowable visit times . Each ship becomes available at a time , at a location . It has a capacity and can traverse any arc for a cost (including fuel and canal costs) and duration . For every ship and node , there is an associated port service cost and duration . There may be incompatibilities between ships and cargoes (e.g.  due to draft limits in the ports). For each and , the boolean defines whether cargo can be serviced by ship . Finally, a penalty is paid if cargo is not transported by the fleet. This penalty corresponds to the revenue loss (or charter cost) due to not transporting an optional cargo.

The objective of ITSRSP is to form routes that minimize the sum of the total travel cost and the possible penalties in the case where charter ships are used or some cargoes are not transported. The routes begin at their respective starting points but have no specified endpoint, since ships operate around the clock. Every route must be feasible: ships cannot exceed their capacity, cargoes can be serviced only within their prescribed time windows, and ships cannot transport incompatible cargoes. Furthermore, the routes must respect pairing and precedence constraints. The pairing constraint states that any pair must belong to the same route, and the precedence constraint states that any pickup must occur before its delivery .

Related literature. Early studies of ship routing and scheduling optimization date back to the 1970–80s. In a seminal study, Ronen (1983) discusses the differences between classical vehicle routing and ship routing and lists possible explanations for the scarcity of research at the time. The author also provides a comprehensive classification scheme for various types of ship routing and scheduling problems. Since this article, research on ship routing has flourished, as highlighted by a general survey of maritime transportation (Christiansen et al. 2007), and reviews focusing on routing and scheduling (Christiansen et al. 2013, Christiansen and Fagerholt 2014).

Many variations of ship routing and scheduling problems have been formulated, and these problems have grown in richness, complexity and accuracy over the years. To name a few, Brown et al. (1987) introduced an SP model to solve a full shipload routing and scheduling problem for a fleet of crude oil tankers. Fagerholt and Christiansen (2000b) proposed a DP algorithm to solve a traveling salesman problem with applications to ship scheduling subproblems. The same algorithm was later exploited by Fagerholt and Christiansen (2000a) to solve subproblems for a multi-ship PDPTW. A maritime PDPTW with split loads was studied by Andersson et al. (2011). The authors proposed two path-flow models and an exact algorithm that generates single ship schedules a priori. Stålhane et al. (2012) studied a ship routing and scheduling problem with split loads, and proposed a BCP algorithm. Vilhelmsen et al. (2014) presented a solution method based on CG to solve a full shipload tramp ship routing and scheduling problem with integrated bunker optimization.

Heuristics and metaheuristics have also been applied to solve several variants of ship routing problems. Some notable examples are the multi-start LS of Brønmo et al. (2007), the unified tabu search of Korsvik et al. (2009), and the large neighborhood searches of Korsvik et al. (2011) and Hemmati et al. (2014). Borthen et al. (2017) used a hybrid genetic search algorithm with great success to solve a multi-period supply vessel planning problem for offshore installations. Furthermore, the UHGS methodology of Vidal et al. (2012, 2014) has led to highly accurate solutions for a considerable number of VRP variants, including the classical capacitated VRP, the VRPTW (Vidal et al. 2013), and several prize-collecting VRP with profits and service selections (Vidal et al. 2016, Bulhões et al. 2018). However, this methodology has never been extended to heterogeneous fixed fleet or pickup-and-delivery problems, which require structurally different neighborhood searches and proper precedence and pairing between the pickups and deliveries in the crossover and split operators.

3 Branch-and-Price Algorithm for the ITSRSP

A simple SP formulation of the ITSRSP is given in Equations (1) to (5). Let be the set of all feasible routes for ship

. This formulation uses a binary variable

to indicate whether or not route of vehicle is used in the current solution for a cost of . Moreover,  is a binary constant that is equal to if and only if the route of ship transports cargo , and otherwise. Each variable is equal to if and only if cargo is transported by a charter instead of being included in a route.


Objective (1) minimizes the routing and charter costs. Constraints (2) ensure that each ship is used at most once, and Constraints (3) guarantee that each cargo is either transported or chartered.

3.1 Column generation

Formulation (15) clearly contains an exponential number of variables, and therefore we will use a CG algorithm to solve its linear relaxation Moreover, each ship in the ITSRSP has a different starting location, cargo compatibility, capacity, travel cost, and time matrix. For this reason, we must solve a collection of pricing subproblems (one for each ship) rather than a single one.

Let and be the dual variables associated with Constraints (2) and (3). The reduced cost of a route, defined in Equation (6), can be distributed into reduced costs for each arc, as shown in Equation (7).


Since the dual costs are originally associated with cargoes (p–d pairs), we opted to associate all the dual costs with the out-arcs from the pickup nodes and depot, and none with those emerging from delivery nodes. This methodological choice led to better performance in our initial experiments, and it is in agreement with previous studies on the PDPTW (e.g.  Ropke and Cordeau 2009).

Pricing. The pricing subproblem is an ESPPRC, which is -hard and often difficult to solve for large instances. Various studies present ways to solve it more efficiently by using route relaxation techniques, at the cost of a slightly weaker linear relaxation (Christofides et al. 1981, Baldacci et al. 2011). Again based on preliminary experiments, we decided to maintain the pairing and precedence constraints, since relaxing these constraints leads to a strong deterioration in the linear relaxation bound. In contrast, elementarity tends to be “naturally” satisfied in most situations without any specific measure in the labeling and dominance. This occurs because after performing a pickup and delivery, the ship is usually far from the pickup point from a spatial and temporal viewpoint, and a new visit to the pickup may be impossible because of the time-window constraints. A similar property has been exploited in Bertsimas et al. (2018) to optimize taxi-fleets services.

We take advantage of this observation by initially relaxing the elementarity and reintroducing it using DSSR, to accelerate the solution of the pricing subproblems (Righini and Salani 2008). This is done by defining a set of pickups that cannot be opened again. The set is initialized as at the start of the process, and it is augmented each time that a repeated service is identified.

The ESPPRC is solved using a forward DP algorithm. For each path , we define a label containing, respectively, the last vertex of the path, the accumulated reduced cost, the total load, the total time, the set of opened p–d pairs, and the set of unreachable pairs. As in Dumas et al. (1991) and Ropke and Cordeau (2009), the set of opened pairs contains the visited pickup nodes for which the corresponding delivery node has not been visited. A pair is unreachable if the pickup node has already been visited. Finally, a valid route is a feasible path such that .

Extending a path to a vertex is allowed only if , , and:


If an extension is allowed, it generates the new label presented in Equation (9):


To reduce the number of labels during the DP algorithm, we use the following dominance rule: a path dominates a path if Condition (10) holds.


Note that implies that . However, as discussed in Ropke and Cordeau (2009), a subset-based dominance between and is valid only if the reduced costs satisfy the DTI: . When we define the reduced costs as in Equation (7), the DTI is valid if the original distances satisfy it. However, it may become violated by branching constraints that introduce new dual costs. Correction techniques need to be applied in these situations, as discussed in Section 3.2.

To the best of our knowledge, the DSSR approach has not been tested for PDPTW problems, but it was suggested as a promising research avenue in Ropke and Cordeau (2009). The same authors also suggested relaxing the sets. We tested this option, but observed that it led to a much slower convergence.

Ship ordering. Since the ITSRSP includes a heterogeneous fleet of ships, it becomes necessary to solve the pricing subproblems associated with each ship type. To reduce the number of subproblems, we tested various approaches based on ship grouping and different orderings. We opted to simply include all the ships in a circular list, and we systematically call the pricing subproblem for the last ship with which a negative-cost route was last obtained. When the pricing algorithm fails to generate a negative reduced-cost route, the procedure selects the next ship in the list. The CG terminates when a full round has been performed without generating any new routes.

Initialization and heuristic pricing. Because of the charter variables in the SP formulation, we simply start with empty sets. To reduce the computational effort, the CG initially uses a fast heuristic pricing in which only the label with the minimum reduced cost for each vertex and time value is kept. The CG starts using the exact pricing when a full round on the ship list fails to generate a new route with the heuristic pricing.

Preprocessing. Finally, our CG exploits various preprocessing techniques to eliminate arcs from the pricing subproblem. A simple version of these strategies was used in Dumas et al. (1991). These procedures are extended to consider the different attributes of the ITSRSP and quickly determine which requests cannot be closed from a given (customer, time) pair, in such a way that it is possible to filter label extensions in the pricing algorithms using bitwise operations in .

3.2 Branch-and-bound

The CG presented in the previous section produces strong lower bounds for the ITSRSP. To obtain integer optimal solutions, we embed it into a branch-and-bound algorithm to form a Branch-and-Price (B&P) method.

Branching rules. The BnP uses three branching rules, giving priority to the most fractional element, as explained below.

  • [nosep]

  • Branching on charters. If any variable is fractional, generate two branches with  and .

  • Branching on ships. If is fractional for any ship , generate two branches with and .

  • Branching on edges for all ships. Given , a binary constant indicating whether or not route of ship  traverses arc , let the number of times any ship traverses or  be . If is fractional, generate two branches with  and .

Delivery triangle inequality. Rule B1 has no impact on the pricing subproblems. In contrast, each new constraint generated by rules B2 and B3 introduces a new dual variable that is included in the reduced costs. The dual variables associated with rule B2 cannot lead to a DTI violation, since the first node of every route must be a pickup node. In contrast, the constraints resulting from rule B3 introduce a dual variable that will be subtracted from the right-hand side of Equation (7), and can lead to violations of the DTI. To circumvent this issue, we use a method similar to that of Ropke and Cordeau (2009) to fix the DTI. To reduce the computational complexity of this approach, we check and fix violations in an incremental manner, focusing on the newly generated dual variables.

Artificial variables. The branching rules presented in this section may make the solution of a child node infeasible. However, it is not immediately possible to be sure about the infeasibility because the solution may simply be missing some columns. For this reason, at every node of the BnP that results in an infeasible solution, the algorithm uses an approach like that of the two-phase simplex method. It introduces an artificial variable on each branching constraint and changes the objective function to minimize their sum, thus minimizing the infeasibility. When the solution becomes feasible again, the artificial variables are removed and the original objective function is restored. If the CG terminates before reaching this state, then the solution is confirmed to be infeasible.

Heuristic strong branching. We apply strong branching to predict which element will result in better solutions, thus reducing the size of the BnP tree. After solving the CG of each node, we build a set of branching candidates from the most fractional elements found by the branching rules, and we simulate the branching for each element by solving both child nodes. Since it is prohibitively expensive to solve the exact pricing several times, we perform heuristic strong branching by executing the CG with the heuristic pricing. Even the non-optimal linear solution gives a good prediction for the quality of each child node and can be used to compare the candidates. The method then retains the best one, i.e.  the branching with the best worst child node. Moreover, a branching with one infeasible child node (from the heuristic pricing viewpoint) is always considered to be better than one with no infeasible child node, and a branching with two infeasible child nodes is immediately chosen. When the branching candidate is chosen, the exact pricing is executed on both child nodes.

3.3 Route enumeration

The techniques discussed to this point lead to an efficient method, the results of which will be discussed in Section 5. Since a good upper bound is known from the heuristic presented in Section 4, we decided to test additional route enumeration techniques that, when applicable, can allow us to solve larger problem instances. Given the known upper bound, the algorithm attempts to enumerate all feasible routes within the integrality gap, and it aborts if more than million routes are created. The route enumeration is done by a DP algorithm similar to that for the exact pricing procedure, using dominance rule (11):


This rule is weaker since it cannot discard a route unless there is another with the exact same set of opened and unreachable customers, and therefore it leads to a much larger number of labels during the DP algorithm. To deal with this issue, we first execute a backward pricing using completion bounds (the best reduced cost of a path ending at a given customer and time) calculated from the last run of the forward exact pricing. From the results of the backward pricing we then calculate backward completion bounds (the best reduced cost of a path starting at a given customer and time). As the name suggests, the backward pricing is the exact pricing algorithm executed from the end of the time horizon to the start, changing dominance rule (10) into () to avoid enforcing the pickup triangle inequality (see Gschwind et al. 2018). At first the completion bounds from the forward pricing subproblems seem to be incompatible with the backward pricing because of the DTI fix. However, we observe that they are an underestimate of the correct completion bounds and therefore can be used to fathom labels. Finally, we generate completion bounds from the backward pricing and use them to fathom labels during the route enumeration procedure.

Upon success, the enumerated routes may be fed into the SP formulation to obtain an optimal solution. However, the number of routes is often prohibitively large to be directly tackled with a mixed integer programming (MIP) solver. Therefore, we continue the search with the same BnP approach and rely on 3-SRC to improve the value of the linear relaxations (Jepsen et al. 2008). For the ITSRSP, the 3-SRC can be obtained from Constraints (3), resulting in the valid inequality presented in (12), where represents the number of times route from vehicle visits the customers of the 3-SRC.


This leads to a BnP algorithm where the separation and pricing is done by simple route inspection, as in Contardo and Martinelli (2014). Note, in addition, that the 3-SRC can be first separated at the root node to improve the linear relaxation and reduce the number of routes resulting from the enumeration. Moreover, 3-SRC are not separated when evaluating candidates during strong branching; they are instead used on the two chosen branches.

4 Hybrid Genetic Search

As demonstrated in Section 5

, the proposed exact approaches lead to remarkable results for industrial-size ship routing instances, but the variance of CPU times can be a drawback for industrial applications. To provide heuristic solutions in a more consistent manner, as well as initial upper bounds for the exact method, we now introduce an HGS, a sophisticated extension of the UHGS of 

Vidal et al. (2012, 2014) which includes problem-tailored search operators and an SP-based intensification procedure.

As Algorithm 1 indicates, our method follows the same general scheme as UHGS with the addition of the SP procedure. It jointly evolves a feasible and an infeasible subpopulation of individuals representing solutions. At each iteration, two parents are selected from the union of the subpopulations. A crossover operator is applied to generate an offspring, which is improved by LS and inserted into the adequate subpopulation according to its feasibility (Lines 3–6). If the offspring is infeasible, a repair procedure is called in an attempt to generate a feasible solution (Lines 7–10).

1 Initialize population;
2 while number of iterations without improvement and  do
3       Select parent solutions and ;
4       Apply the crossover operator on and to generate an offspring ;
5       Educate offspring by local search;
6       Insert into respective subpopulation;
7       if  is infeasible then

    With probability

, repair (local search) and
9             insert it into respective subpopulation;
11      if maximum subpopulation size reached then
12             Select survivors;
14      if best solution not improved for iterations then
15             Diversify population;
17      if best solution not improved for iterations then
18             Run set partitioning;
20      Adjust penalty coefficients for infeasibility;
22Return best feasible solution;
Algorithm 1 Hybrid Genetic Search (HGS)

Whenever a subpopulation reaches a maximum size, a survivor selection procedure is applied to remove individuals based on their quality and contribution to the population diversity (Lines 11–13). If no improving solution is found after successive iterations, new individuals are added to the population in order to diversify it (Lines 14–16). Similarly, if no improving solution is found after successive iterations, an intensification procedure is triggered, in the form of an SP model that aims to construct better solutions from high-quality routes identified in the search history. Finally, the penalty coefficients are periodically adjusted to control the proportion of feasible individuals in the search.

Hybrid genetic searches with a similar structure have been used with great success for a variety of VRP (see, e.g., Vidal et al. 2014, 2016, Borthen et al. 2017). Still, the ITSRSP includes such a diversity of constraints that most components of the method had to be carefully tailored in order to obtain an effective algorithm. The following subsections describe each component in more detail.

4.1 Search space

Previous studies have demonstrated that the controlled use of penalized infeasible solutions can help converging towards high-quality feasible solutions (Glover and Hao 2009, Vidal et al. 2015a). This is especially relevant for the ITSRSP, since this problem includes time windows, capacity constraints and incompatibility constraints, as well as precedence and pairing restrictions for p–d pairs. In the proposed HGS, we allow the exploration of infeasible solutions in which:

  • [nosep]

  • ship capacity constraints may be exceeded;

  • some cargoes may not be picked up or delivered within their respective time windows;

  • each ship may transport incompatible cargoes; but

  • no component of the HGS creates solutions that violate precedence or pairing constraints.

The load infeasibility is proportional to the difference between the peak load (largest load over the trip) and the ship capacity. To relax the time-window constraints, we use the “time-warp” approach of Nagata et al. (2010) and Vidal et al. (2013), which relaxes the precedence constraints between visits by allowing penalized “returns in time” upon a late arrival to a customer. Finally, the penalty associated with ship-cargo incompatibilities is proportional to the number of incompatible cargoes carried by each ship.

Let be a route for ship , starting from the initial position () and servicing a sequence of (pickup or delivery) nodes . The start-of-service time at the th node can be defined as:


Route can be characterized by the following quantities:

Travel cost: (14)
Peak load: (15)
Time warp use: (16)
Incompatibilities: (17)

Finally, we define the penalized cost of route for ship as:


where , , and are the respective penalty coefficients for peak-load, time-window, and incompatibility-constraint violations. The penalty coefficients will be adjusted during the search as described in Section 4.5. The penalized cost of a solution is the sum of the penalized costs of all its routes, that is, .

4.2 Solution representation and evaluation

A solution is represented in HGS as a giant tour that holds a permutation of nodes in and satisfies the precedence constraints between the p–d pairs. Such a representation greatly facilitates the design of an effective crossover operator. Moreover, segmenting this giant tour into different routes can be done efficiently using a variant of the Split algorithm (Prins 2004, Vidal 2016) to form a complete solution after crossover.

Split is a DP algorithm that was originally designed for the capacitated VRP but is flexible enough to be adapted to a variety of constraints and objectives (see, e.g., Prins et al. 2009, Velasco et al. 2009). When dealing with VRP with heterogeneous fleets, previous authors have assumed that Split should jointly optimize the giant tour segmentation and the choice of ship for each route (Duhamel et al. 2011). This extension, unfortunately, leads to a special case of the shortest path problem with resource constraints, for which only pseudo-polynomial algorithms are currently available. To avoid this issue, we opted to fix the sequence of ships and restrict the Split algorithm to the segmentation of the tours, so that the ships are considered one by one in their order of appearance. To avoid any possible bias from the instance representation, we shuffle the order of the ships when reading the data and keep this order fixed during the solution process. The possible use of spot charters for optional cargoes is modeled via a dummy ship of index with zero distance cost and a service cost equal to the charter price for each cargo.

The Split graph is defined as follows. Let be a directed acyclic graph with nodes and arcs . Each arc represents a route


associated with ship . If , then , representing an empty route. The cost of arc is set to when the route satisfies the pairing constraints (no open pickup or delivery), and infinity otherwise. With these definitions, an optimal segmentation of the giant tour into routes assigned to the ships to corresponds to a shortest path between nodes and  in . This shortest path can be obtained via Bellman’s algorithm in topological order, with a time complexity of and space complexity of . Moreover, note that there always exists at least one feasible path without necessarily relying on spot charters: a single route for ship containing all visits naturally satisfies the pairing and precedence restrictions, despite its high cost related to the penalized violation of all the other (load, time, and incompatibility) constraints.

Individual evaluation. As in UHGS, the quality of an individual is not based solely on its cost but also on its contribution to the subpopulation diversity. The combination of these two metrics is referred to as the biased fitness of in its subpopulation , and it is defined as


where is the penalized cost rank of in , and is the diversity contribution rank of in . Both ranks are relative to the subpopulation size, and parameter balances the weight of each rank. The diversity contribution of in is defined as the average distance to its closest individuals. We use the broken pairs distance, measuring the proportion of different edges between two solutions.

4.3 Parent selection and crossover

At each iteration, the algorithm selects two parents and by binary tournament based on their biased fitness. To produce a child , these parents are submitted to a specialized one-point crossover operator (Velasco et al. 2009), designed to enforce the pairing and precedence constraints between the pickups and deliveries:

  • [nosep,leftmargin=0cm]

  • Step 1) A cutting point is randomly selected with uniform probability in the giant tour of the first parent, and the sequence of visits is copied into .

  • Step 2) The second parent is swept from beginning to end, and any pending delivery ( such that ) is inserted at the end of .

  • Step 3) The second parent is swept a second time, and any missing node is inserted at the end of .

4.4 Education and repair

Each individual resulting from the crossover operator is decoded with the Split algorithm to obtain a complete solution, and then improved (i.e., educated) by an efficient LS based on a variety of neighborhoods tailored for pickup-and-delivery problems. The moves are evaluated in random order, and any improving move is directly applied (first-improvement strategy). The LS terminates as soon as no improving move exists. After thorough computational analyses, we selected five neighborhoods:

  • [nosep]

  • Relocate Pickup: Relocate a pickup in the same route after a node located before the corresponding delivery .

  • Relocate Delivery: Relocate a delivery in the same route after a node located after the corresponding pickup .

  • Relocate Pair: Relocate a p–d pair , placing after a node and placing no more than nodes after .

  • Swap Pair: Given two pairs and , swap with and with .

  • Swap Ships: Exchange the ships assigned to two different routes.

All these neighborhoods preserve the p–d pairing and precedence constraints. Neighborhoods and specialize in intra-route modifications, while and may be used for both intra- and inter-route modifications, and they can therefore change cargo-ship allocations. Finally, to maintain a low complexity, neighborhood is limited to .

Efficient move evaluations. All the move evaluations are performed in amortized time thanks to concatenation strategies (Vidal et al. 2014, 2015b). These strategies are based on the fact that all moves in to create routes that correspond to the concatenation of a constant number of route subsequences of the current solution. Therefore, preprocessing meaningful information on subsequences of consecutive visits prior to move evaluations (as well as after each route update) can facilitate the evaluation of complex constraints and objectives.

Name Base case Induction step – Concatenation
Travel cost
Peak load
Time warp use
Earliest possible completion time
Latest feasible starting time
Auxiliary computations
Table 1: Preprocessing and move evaluations by concatenation.

The information preprocessed on subsequences for the ITSRSP is listed in Table 1. This is done by induction on the operation of concatenation of two visit sequences, starting from the base case of a single node in the sequence. Note that, in contrast with other problems, this information depends on the ship type , requiring preprocessing time and space if a brute force approach is employed, since a solution contains subsequences of consecutive visits and  ships. Fortunately, this complexity can be reduced to time and space by observing that the following information is sufficient to evaluate all the moves:

  • [nosep,leftmargin=*]

  • The information on all subsequences of consecutive nodes in the incumbent solution, for their current ship type (for neighborhoods to ) ;

  • The information on each single node for all ship types, in (for and );

  • The information on each sequence representing a complete route for all ship types, which can be computed in time and stored in (for ).

Finally, to avoid redundant move evaluations, the HGS uses a simple memory scheme that registers the last-modified time of a route and the last-evaluated time for each move. By comparing these values, one can decide whether or not to re-evaluate a move. This strategy is as efficient as and much simpler than the “static move descriptors” discussed in Zachariadis and Kiranoudis (2010). Moreover, time stands for any non decreasing counter, e.g., the number of moves applied or tested in the method.

Vehicle-dependent neighborhood restrictions. Our LS uses static neighborhood restrictions similar to those of Vidal et al. (2013), by limiting move evaluations to those that create at least one directed arc with vehicle such that and or and . The set contains the most promising successors of vertex for ship , ranked according to a metric of spatial and temporal proximity between nodes and :


The first term of the equation represents the spatial proximity (distance), scaled by the ratio between travel time and distance to ensure that all terms have the same unit. The next two terms measure the temporal proximity, based on the unavoidable amount of waiting time and time warp when servicing and consecutively, with weights and .

Repair. The routes and solutions explored in the LS can include penalized violations of time-window, capacity, and incompatibility constraints. Therefore, this procedure may lead to an infeasible solution that will be stored in the infeasible subpopulation. In this event, a Repair phase is additionally called on this solution with probability . Repair temporarily multiplies all the penalty coefficients by and runs the LS. If the resulting solution remains infeasible, then the coefficients are again multiplied, this time by , and the LS is run again. If it is successful, the resulting feasible solution will be added to the feasible subpopulation.

4.5 Population management

As in UHGS, we rely on survivor selection, population diversification, and adaptive penalty mechanisms to find a good balance between population diversity and elitism. We also incorporate an additional intensification phase, in the form of an SP procedure that aims to build a better solution from existing routes from the search history.

Initialization. To initialize the population, the HGS generates random individuals. Random individuals are generated as giant tours where a sequence of pickups is shuffled and then deliveries are placed immediately after the respective pickups. These individuals are educated, possibly repaired, and inserted into their respective subpopulations.

Survivor selection and diversification. A survivor selection mechanism occurs whenever a subpopulation reaches the maximum size of individuals. As in Vidal et al. (2012), the individuals with maximum biased fitness are discarded, prioritizing individuals that have a clone. This selection procedure preserves the best individuals with respect to the penalized cost and finds a good balance between elitism and diversity. To explore an even wider diversity of solution characteristics, a diversification procedure is called whenever no improving solution was found during the last iterations. It discards all but the individuals with the smallest biased fitness in each subpopulation, and then generates new random individuals that are educated and possibly repaired.

Set partitioning. Because of the many constraints of the ITSRSP, even generating promising feasible routes can be a challenging task. In this context, it is natural to exploit as far as possible high-quality routes from the search history. Thus, in a similar way to Muter et al. (2010) and Subramanian et al. (2013), HGS triggers an intensification procedure whenever no improving solution has been found over consecutive iterations. This procedure formulates an SP model (Equations 15) that considers all feasible routes from local minima in the past search, and solves it with a MIP solver subject to a time limit of . Any improved solution obtained is inserted into the population. This intensification procedure complements the other operators well: instead of performing local improvements, it seeks good combinations of previously found routes and can be viewed as a form of large neighborhood search.

Adaptive penalties. Finally, the penalty coefficients are adjusted during the search to maintain a target proportion of feasible individuals with respect to the relaxed constraints (load, time windows, and incompatibilities). For each constraint , HGS records the proportion of feasible solutions w.r.t. constraint over the last iterations. This calculation is performed after the LS, before any possible repair. The penalty coefficients are then adjusted every iterations: if , is multiplied by , and if , is multiplied by .

5 Experimental analysis

This section reports our computational experiments with the two BnP algorithms (with and without enumeration) and the HGS. Our aim is threefold:

  • [nosep,leftmargin=*]

  • We compare the proposed exact algorithms, and evaluate their ability to solve practical-size ship routing problems to optimality in limited time;

  • In situations where a faster response is sought, we evaluate the quality of the solutions produced by the HGS metaheuristic;

  • Finally, we evaluate the impact of some of our most important methodological choices and new components: e.g., vehicle-dependent neighborhood restrictions and set-partitioning problem parameters for the HGS; advanced preprocessing techniques, DSSR, and completion bounds for the BnP algorithms.

We implemented all algorithms in C++ with double precision numbers, using CPLEX 12.7 to solve the BnP master problem and the integer SP inside the HGS. We conducted all experiments on a computer with an i7-3960X CPU and 64 GB of RAM.

We rely on the benchmark instances suite for the ITSRSP based on real-life scenarios, presented in Hemmati et al. (2014) and currently available at These instances are divided into four groups of , according to problem topology and cargo type: short sea mixed load (SS_MUN), short sea full load (SS_FUN), deep sea mixed load (DS_MUN), and deep sea full load (DS_FUN). Each group contains five instances for 12 different problem size values. The mixed-load instances have up to cargoes and ships, whereas the full-load instances have up to cargoes and ships. So far, instances remain open. In the full-load instances, each delivery should be visited immediately after its associated pickup due to the absence of residual capacity for other loads. This property is no longer valid for mixed-load instances. Furthermore, short and deep sea instances consider different geographical regions. The short sea instances represent shipments among European ports, whereas the deep sea instances involve long-distance shipments between different continents.

5.1 Exact solutions

To our knowledge, the study of Hemmati et al. (2014) is the only one to report lower bounds and optimal solutions for these benchmark instances, obtained by solving a MIP formulation. We establish a comparison with our branch-and-price without route enumeration (B&P1), using heuristic strong branching with 50 candidates at each iteration, as well as the same algorithm with route enumeration and 3-SRC using inspection pricing (B&P2). As an initial upper bound, we set , where is the best objective value found by our HGS metaheuristic. One hour of computation was allowed for each instance.

Table 2 reports the experimental results. Each line gives the average results for five instances with the same characteristics. The first group of columns presents the results of Hemmati et al. (2014): the time in minutes, the gap between the integer solution and the best bound found in the MIP formulation, and the number of instances solved to optimality. The second group presents the results for B&P1: “Gap0” and “T0” represent the percentage gap and the time in minutes for the root node. “GapF” and “TF” are the final percentage gap and time, “NF” is the number of nodes in the search tree, and “Opt” is the number of instances solved to optimality. The last group of columns presents the results for B&P2 (with route enumeration and inspection pricing): the columns “TE” and “RE” are the time in minutes for the route enumeration and the number of routes found. “Gap0”, “T0”, and “Cuts0” are the percentage gap, the overall time and the number of 3-SRC separated at the root node. “RF”, “GapF”, “TF”, “CutsF”, and “NF” are the final number of routes, the percentage gap, the overall time, the 3-SRC separated, and the number of nodes in the search tree. Finally, “Opt” is the number of instances solved to optimality in the group.

Hemmati et al. (2014) Branch-and-Price (B&P1) Branch-and-Price + Enumeration + SRCs (B&P2)
Instances T Gap Opt Gap0 T0 GapF TF NF Opt TE RE Gap0 T0 Cuts0 RF GapF TF CutsF NF Opt


C7-V3 0.0 0.00 5 0.00 0.00 0.00 0.00 1.0 5 0.00 9 0.00 0.00 0.0 9 0.00 0.00 0.0 1.0 5
C10-V3 0.0 0.00 5 0.00 0.00 0.00 0.00 1.0 5 0.00 13 0.00 0.00 0.0 13 0.00 0.00 0.0 1.0 5
C15-V4 1.4 0.00 5 1.07 0.00 0.00 0.00 3.4 5 0.00 75 0.80 0.00 3.2 64 0.00 0.00 4.6 1.8 5
C18-V5 42.5 3.12 4 0.56 0.00 0.00 0.00 5.4 5 0.00 83 0.16 0.00 2.2 40 0.00 0.00 2.2 1.8 5
C22-V6 50.6 15.47 1 1.40 0.00 0.00 0.01 6.6 5 0.00 468 0.61 0.00 13.8 242 0.00 0.00 17.0 4.2 5
C23-V13 60.0 26.88 0 0.41 0.01 0.00 0.05 7.8 5 0.00 180 0.17 0.01 5.8 74 0.00 0.01 5.8 3.0 5
C30-V6 60.0 79.46 0 0.73 0.01 0.00 0.11 15.0 5 0.00 913 0.36 0.01 25.8 290 0.00 0.01 31.4 2.6 5
C35-V7 60.0 82.66 0 0.66 0.03 0.00 0.36 23.0 5 0.00 2K 0.38 0.03 42.4 850 0.00 0.07 78.6 7.0 5
C60-V13 60.0 85.48 0 0.43 0.49 0.00 11.82 96.6 5 0.07 28K 0.29 0.48 57.2 13K 0.00 0.91 176.0 20.6 5
C80-V20 60.0 86.63 0 0.19 1.20 0.00 9.06 17.4 5 0.19 35K 0.08 1.15 51.2 12K 0.00 1.23 79.0 5.4 5
C100-V30 60.2 97.59 0 0.15 2.88 0.02 36.02 40.0 3 0.47 23K 0.10 2.81 34.4 9K 0.00 3.27 70.0 19.0 5
C130-V40 61.5 99.95 0 0.13 12.68 0.04 60.00 18.4 0 2.82 111K 0.07 13.44 41.8 36K 0.00 15.08 115.6 23.8 5
Overall 43.0 48.10 20 0.48 1.44 0.01 9.79 19.6 53 0.30 17K 0.25 1.49 23.2 6K 0.00 1.71 48.4 7.6 60


C8-V3 0.0 0.00 5 0.00 0.00 0.00 0.00 1.0 5 0.00 10 0.00 0.00 0.0 10 0.00 0.00 0.0 1.0 5
C11-V4 0.0 0.00 5 0.00 0.00 0.00 0.00 1.0 5 0.00 19 0.00 0.00 0.0 19 0.00 0.00 0.0 1.0 5
C13-V5 0.0 0.00 5 0.00 0.00 0.00 0.00 1.0 5 0.00 21 0.00 0.00 0.0 21 0.00 0.00 0.0 1.0 5
C16-V6 0.0 0.00 5 0.00 0.00 0.00 0.00 1.0 5 0.00 34 0.00 0.00 0.0 34 0.00 0.00 0.0 1.0 5
C17-V13 0.0 0.00 5 0.01 0.00 0.00 0.00 1.8 5 0.00 49 0.01 0.00 0.0 49 0.00 0.00 0.0 1.8 5
C20-V6 0.8 0.00 5 0.03 0.00 0.00 0.00 1.8 5 0.00 81 0.00 0.00 2.6 76 0.00 0.00 5.0 1.4 5
C25-V7 40.5 0.83 3 0.00 0.00 0.00 0.00 1.4 5 0.00 67 0.00 0.00 0.8 62 0.00 0.00 0.8 1.0 5
C35-V13 60.0 8.85 0 0.00 0.00 0.00 0.00 1.4 5 0.00 131 0.00 0.00 1.4 130 0.00 0.00 1.4 1.4 5
C50-V20 60.0 13.99 0 0.03 0.00 0.00 0.01 1.8 5 0.00 671 0.03 0.00 0.0 671 0.00 0.01 0.0 1.8 5
C70-V30 60.1 60.04 0 0.12 0.01 0.00 0.08 4.6 5 0.00 7K 0.12 0.01 0.0 7K 0.00 0.04 3.0 4.6 5
C90-V40 60.3 78.32 0 0.00 0.03 0.00 0.11 2.2 5 0.00 1K 0.00 0.03 1.0 1K 0.00 0.04 1.0 1.8 5
C100-V50 60.9 79.11 0 0.01 0.04 0.00 0.24 3.8 5 0.00 3K 0.01 0.04 2.0 2K 0.00 0.10 2.0 3.4 5
Overall 28.5 20.10 33 0.02 0.01 0.00 0.04 1.9 60 0.00 1K 0.02 0.01 0.7 978 0.00 0.02 1.1 1.8 60


C7-V3 0.0 0.00 5 1.15 0.00 0.00 0.00 1.4 5 0.00 12 0.00 0.00 0.6 8 0.00 0.00 0.6 1.0 5
C10-V3 0.0 0.00 5 1.92 0.00 0.00 0.00 3.4 5 0.00 30 1.81 0.00 1.0 28 0.00 0.00 1.0 1.8 5
C15-V4 0.4 0.00 5 1.02 0.00 0.00 0.00 1.8 5 0.00 80 0.00 0.00 1.4 22 0.00 0.00 1.4 1.0 5
C18-V5 16.3 2.18 4 0.51 0.00 0.00 0.00 2.6 5 0.00 97 0.25 0.00 5.0 78 0.00 0.00 5.0 1.4 5
C22-V6 26.2 3.28 4 1.45 0.00 0.00 0.00 5.8 5 0.00 154 1.01 0.00 2.8 108 0.00 0.00 3.4 2.6 5
C23-V13 26.4 4.62 3 0.27 0.00 0.00 0.00 2.6 5 0.00 79 0.13 0.00 0.4 55 0.00 0.00 0.4 1.8 5
C30-V6 60.0 52.25 0 1.56 0.00 0.00 0.04 15.0 5 0.00 4K 0.84 0.00 16.0 1K 0.00 0.01 30.6 3.0 5
C35-V7 60.0 54.25 0 1.13 0.01 0.00 0.06 14.6 5 0.00 8K 0.74 0.01 20.6 5K 0.00 0.01 25.6 3.4 5
C60-V13 60.0 89.22 0 0.58 0.10 0.00 3.62 45.4 5 0.02 33K 0.14 0.10 32.8 5K 0.00 0.70 88.2 9.4 5
C80-V20 60.1 91.56 0 0.43 0.37 0.00 12.57 49.4 5 0.10 231K 0.11 0.42 43.4 30K 0.00 0.59 64.0 7.8 5
C100-V30 60.3 99.03 0 0.52 0.57 0.09 20.72 63.4 4 0.21 1M 0.24 0.70 32.0 395K 0.00 2.13 99.4 22.2 5
C130-V40 61.3 100.00 0 0.41 3.81 0.16 60.00 30.0 0 4.86 13M 0.28 9.07 52.6 5M 0.01 17.64 168.0 19.2 4
Overall 35.9 41.37 26 0.91 0.40 0.02 8.08 19.6 54 0.43 1M 0.46 0.86 17.4 478M 0.00 1.76 40.6 6.2 59


C8-V3 0.0 0.00 5 0.00 0.00 0.00 0.00 1.0 5 0.00 10 0.00 0.00 0.0 10 0.00 0.00 0.0 1.0 5
C11-V4 0.0 0.00 5 0.00 0.00 0.00 0.00 1.0 5 0.00 14 0.00 0.00 0.0 14 0.00 0.00 0.0 1.0 5
C13-V5 0.0 0.00 5 0.00 0.00 0.00 0.00 1.0 5 0.00 18 0.00 0.00 0.0 18 0.00 0.00 0.0 1.0 5
C16-V6 0.0 0.00 5 0.01 0.00 0.00 0.00 1.4 5 0.00 28 0.01 0.00 0.0 28 0.00 0.00 0.0 1.4 5
C17-V13 0.0 0.00 5 0.02 0.00 0.00 0.00 1.4 5 0.00 34 0.02 0.00 0.0 34 0.00 0.00 0.0 1.4 5
C20-V6 0.0 0.00 5 0.11 0.00 0.00 0.00 1.4 5 0.00 67 0.11 0.00 0.0 67 0.00 0.00 0.0 1.4 5
C25-V7 0.1 0.00 5 0.00 0.00 0.00 0.00 1.0 5 0.00 80 0.00 0.00 0.0 80 0.00 0.00 0.0 1.0 5
C35-V13 45.6 4.37 2 0.03 0.00 0.00 0.00 1.8 5 0.00 145 0.03 0.00 0.0 145 0.00 0.00 0.0 1.8 5
C50-V20 56.3 8.81 1 0.00 0.00 0.00 0.00 1.0 5 0.00 182 0.00 0.00 0.0 182 0.00 0.00 0.0 1.0 5
C70-V30 60.1 11.37 0 0.00 0.01 0.00 0.01 1.4 5 0.00 682 0.00 0.01 0.4 682 0.00 0.01 0.4 1.4 5
C90-V40 60.2 48.44 0 0.00 0.02 0.00 0.04 1.8 5 0.01 1K 0.00 0.02 1.4 1K 0.00 0.03 5.2 1.8 5
C100-V50 60.5 52.36 0 0.00 0.02 0.00 0.07 1.8 5 0.01 682 0.00 0.02 0.6 651 0.00 0.04 1.0 1.4 5
Overall 23.6 10.45 38 0.01 0.00 0.00 0.01 1.3 60 0.00 267 0.01 0.00 0.2 265 0.00 0.01 0.6 1.3 60
Table 2: Performance comparison – Exact approaches

As observed in Table 2, the MIP formulation already fails to produce optimal solutions on some small instances with 20 to 30 cargoes. In contrast, B&P1 solves all the full cargo instances in less than one minute, as well as 107 out of the 120 mixed cargo instances. The column generation produces good lower bounds, with an average gap of at the root node. However, for the mixed-cargo instances, the average time needed to solve each pricing subproblem (ratio TF/NF) increases with the number of cargoes, whereas the lower bounds tend to deteriorate, leading to larger branch-and-bound trees. To go further, one can concentrate on improving the pricing problem solution or the lower bounds. For the largest open instances, the root-node gap seems small enough to allow a complete route enumeration via sophisticated DP algorithms (Section 3.3). We therefore derived B&P2 from this premise: the enumeration of the routes allows subsequent pricing by inspection and permits to introduce SRCs without significant consequences on CPU time.

The remaining columns of Table 2 report the performance of B&P2. This approach outperforms the standard B&P1 in terms of number of instances solved and average solution time. Route enumeration at the root node takes a maximum of 17.8 minutes. The number of routes may, however, rise up to 60 millions (on instance DS-MUN-C130-V40-HE-2). Having enumerated the routes allows to efficiently separate the SRCs and decrease the root-node gap (from 0.36% to 0.19%). As a consequence, the route set can be reduced further (from 311K to 121K on average) as well as the size of the branch-and-bound tree (from 10.6 to 4.2 on average). Using this method, all available instances but one (DS-MUN-C130-V40-HE-2) are solved within a time limit of one hour. A larger time limit allows to solve the last remaining instance in 4 hours and 23 minutes.

Figure 1: Number of instances solved over time by B&P1 and B&P2

Figure 1 displays the number of instances solved by B&P1 and B&P2 as a function of the CPU time limit. B&P2 visibly produces superior results, but this method is also less flexible: its performance depends on the ability to do a complete route enumeration at the root node within the optimality gap. This process may possibly fail (due to time or memory limitations, or an unusually large gap), such that we recommend to first consider B&P1 in exploratory analyses without any prior knowledge of the structure of the ITSRSP instances.

5.2 Heuristic solutions

As seen in the previous section, our exact methods can solve the majority of the instances, but their CPU time can widely vary, even for instances with similar characteristics. Therefore, fast heuristic solutions remain essential for applications requiring a response in a guaranteed short time. This section compares the performance of our HGS with that of the existing ITSRSP heuristics from Hemmati et al. (2014) and Hemmati and Hvattum (2016). To highlight the impact of the SP component, we evaluated two versions of our algorithm: without SP (HGS1) and with SP (HGS2).

We used the same parameters as Vidal et al. (2013) to avoid any problem-specific overfit. Therefore, , , , , , , and . The penalty coefficients take as a starting value, where and . To compare with previous literature in similar CPU time, we set , and . Finally, the full-load instances (FUN) require directs visits from pickup to delivery points, such that neighborhoods and can be disabled and neighborhood is restricted to . In the neighborhood restrictions, is the only possible candidate successor of , and only a pickup node can follow a delivery .

Table 3 compares the results of the ALNS of Hemmati et al. (2014) with those of HGS1 and HGS2. Each line gives average results over 10 runs for five instances with the same characteristics. Each column “Gap” reports the percentage gap of one method relative to the optimal solutions found in Section 5.1, column “T” gives the total CPU time in minutes, and column “T*” gives the attainment time (in minutes) needed to reach the final solution. For each instance class, the best gap is highlighted in boldface. The detailed results of Hemmati et al. (2014) were kindly provided by the authors, and our complete detailed results (per instance) are also available in the electronic companion of this paper, located at

Hemmati et al. (2014) HGS1 HGS2
Instances Gap T Gap T T* Gap T T*


C7-V3 0.00 0.03 0.00 0.02 0.00 0.00 0.02 0.00
C10-V3 0.00 0.04 0.00 0.04 0.00 0.00 0.04 0.00
C15-V4 0.58 0.09 0.00 0.08 0.00 0.00 0.08 0.00
C18-V5 0.51 0.13 0.07 0.12 0.02 0.02 0.13 0.03
C22-V6 1.82 0.19 0.05 0.19 0.05 0.00 0.18 0.04
C23-V13 0.58 0.25 0.13 0.25 0.09 0.00 0.23 0.07
C30-V6 1.60 0.38 0.17 0.36 0.17 0.00 0.33 0.12
C35-V7 1.92 0.54 0.81 0.51 0.26 0.03 0.54 0.27
C60-V13 1.69 2.01 1.03 2.44 1.81 0.14 2.25 1.37
C80-V20 2.45 4.13 1.79 5.21 4.09 0.03 3.84 2.44
C100-V30 2.68 7.77 1.67 9.97 8.01 0.01 5.31 3.02
C130-V40 2.57 16.95 2.28 13.84 11.87 0.02 12.43 8.39
Overall 1.37 2.71 0.67 2.75 2.20 0.02 2.12 1.31


C8-V3 0.00 0.03 0.00 0.01 0.00 0.00 0.01 0.00
C11-V4 0.13 0.05 0.00 0.02 0.00 0.00 0.02 0.00
C13-V5 0.07 0.07 0.00 0.02 0.00 0.00 0.02 0.00
C16-V6 0.05 0.10 0.00 0.03 0.00 0.00 0.03 0.00
C17-V13 0.01 0.14 0.00 0.04 0.00 0.00 0.05 0.00
C20-V6 0.14 0.18 0.00 0.04 0.00 0.00 0.05 0.00
C25-V7 0.22 0.27 0.00 0.07 0.02 0.00 0.07 0.01
C35-V13 0.29 0.60 0.01 0.19 0.08 0.00 0.17 0.05
C50-V20 0.52 1.38 0.22 0.67 0.43 0.00 0.43 0.15
C70-V30 1.29 3.51 0.68 1.86 1.28 0.00 1.09 0.42
C90-V40 1.45 6.98 0.70 3.90 2.83 0.00 1.93 0.76
C100-V50 0.96 9.79 0.50 5.98 4.46 0.00 2.83 1.20
Overall 0.43 1.93 0.18 1.07 0.76 0.00 0.56 0.22


C7-V3 0.00 0.03 0.00 0.02 0.00 0.00 0.02 0.00
C10-V3 0.01 0.04 0.00 0.04 0.00 0.00 0.04 0.00
C15-V4 1.26 0.08 0.00 0.08 0.01 0.00 0.08 0.01
C18-V5 0.47 0.13 0.00 0.12 0.01 0.00 0.12 0.01
C22-V6 2.18 0.19 0.01 0.18 0.05 0.00 0.18 0.05
C23-V13 0.12 0.24 0.02 0.19 0.05 0.00 0.20 0.04
C30-V6 1.04 0.37 0.22 0.36 0.17 0.03 0.32 0.12
C35-V7 1.08 0.51 0.14 0.53 0.28 0.00 0.44 0.17
C60-V13 3.74 1.92 1.41 2.58 1.96 0.03 1.81 1.06
C80-V20 3.10 4.26 1.61 5.69 4.58 0.04 3.56 2.20
C100-V30 3.69 8.00 3.27 9.22 7.27 0.01 9.87 6.84
C130-V40 5.18 17.47 5.08 13.66 11.64 0.09 14.81 12.88
Overall 1.82 2.77 0.98 2.72 2.17 0.02 2.62 1.95


C8-V3 0.00 0.03 0.00 0.01 0.00 0.00 0.01 0.00
C11-V4 0.00 0.05 0.00 0.02 0.00 0.00 0.02 0.00
C13-V5 0.00 0.06 0.00 0.02 0.00 0.00 0.04 0.00
C16-V6 0.03 0.10 0.00 0.03 0.00 0.00 0.04 0.00
C17-V13 0.00 0.13 0.00 0.05 0.00 0.00 0.05 0.00
C20-V6 0.01 0.16 0.00 0.04 0.00 0.00 0.05 0.00
C25-V7 0.41 0.26 0.00 0.07 0.01 0.00 0.07 0.01
C35-V13 1.03 0.59 0.01 0.26 0.14 0.00 0.21 0.08
C50-V20 0.61 1.41 0.23 0.71 0.46 0.00 0.45 0.17
C70-V30 0.59 3.55 0.31 2.02 1.43 0.00 1.02 0.38
C90-V40 1.10 7.01 0.44 4.27 3.16 0.00 2.20 1.00
C100-V50 1.07 9.85 0.41 6.37 4.80 0.00 3.13 1.44
Overall 0.41 1.93 0.12 1.16 0.83 0.00 0.61 0.26
Table 3: Performance comparison – Metaheuristics

As visible in the results of Table 3, both versions of the HGS largely outperform previous algorithms. HGS2 obtains near-optimal solutions within an average gap of , compared to for the ALNS. The full-load instances are generally easier to solve than the mixed-load instances since the problem simplifies. Moreover, the SP-based procedure largely contributes to the performance of the algorithm: contrary to intuition, it does not increase the overall CPU time, but even decreases it from  min (HGS1) to  min (HGS2) on average. Indeed, the SP helps to reach optimal solutions more quickly, such that the method only performs additional iterations before reaching the termination criteria instead of gradually improving over a longer time. The effect is manifest when comparing the average attainment time of HGS1 (T* ) with that of HGS2 (). In a follow-up work, Hemmati and Hvattum (2016) investigated the impact of some design decisions and parameters of their ALNS and reported the results of six variants of their algorithm on a smaller subset of instances. Drawing a comparison with these methods leads to similar conclusions: HGS2 largely outperforms these methods. For the sake of brevity, this comparison is presented in the electronic companion.

5.3 Sensitivity analyses

Finally, to highlight the role of each main strategy and parameter, we have conducted extensive sensitivity analyzes with the BnP and HGS algorithms. We started with the standard configurations of each method (B&P1, B&P2 and HGS2) and generated a number of alternative configurations by deactivating a component or modifying a single parameter to study its impact. The results of these analyzes are presented in Tables 45. In Table 4, Column “Gap” represent the average gap, Column “T” reports the average CPU time in minutes, Column “Root” counts the number of instances for which the root node solution was completed, and Column “Opt” counts the number of optimal solutions found. By convention, failing to solve the root node gives a Gap of 100%. In Table 5, Column “Gap” refers to the average gap, and Columns “T” and “T*” represent the average CPU time and attainment time.

Gap T Root Opt Gap T Root Opt
Standard (B&P1) 0.00 0.02 120 120 0.01 8.94 120 107
A. No Heuristic Pricing 0.00 0.02 120 120 0.02 11.20 120 103
B. No Strong Branching 0.00 1.09 120 118 0.06 17.54 120 86
C. No Preprocessing 0.00 0.03 120 120 6.76 15.10 112 95
D. No DSSR 0.00 0.04 120 120 25.02 19.27 90 84
Standard (B&P2) 0.00 0.01 120 120 0.00 1.74 120 119
E. No Completion Bounds 0.00 0.01 120 120 13.33 9.97 104 104
F. No Subset-Row Cuts 0.00 0.01 120 120 0.00 2.47 120 118
Table 4: Impact of some of the key components of the B&P
Gap T T* Gap T T* 0.00 0.58 0.24 0.9450 0.02 2.37 1.63 0.7908
G. No SP Intensification (HGS1) 0.15 1.11 0.80 0.82 2.74 2.18 0.00 0.58 0.24 0.9450 0.02 2.29 1.54 0.7867
I. Longer SP Intensification () 0.00 0.58 0.24 0.02 2.37 1.59 0.00 0.72 0.29 0.9675 0.05 3.49 2.35 0.8400
K. No diversity management () 0.00 0.58 0.25 0.03 2.39 1.71
Table 5: Impact of some of the key components and parameters of the HGS

In these experiments, again, the instances of class FUN are solved more easily. Since most methods achieve the same solution quality for this class, we will primarily rely on the harder MUN instances in our analyzes.

Table 4 analyzes the components of the B&P algorithms. As illustrated by the results of Configuration A, heuristic pricing significantly reduces the overall pricing time. Without this component, four additional instances remain open and the CPU time increases by 25%, though this effect is generally less marked than in other VRP variants (see, e.g. Desaulniers et al. 2008, Martinelli et al. 2014).

Deactivating strong branching (Configuration B) has a larger impact. Strong branching helps predicting good branching decisions and reducing the search tree. Without it, the number of solved instances decreases from 107 to 86. The method can still nearly close the optimality gap (0.06% on average), but it often fails to complete the optimality proof.

Configurations C and D evaluate the impact of our advanced preprocessing strategies and DSSR (Section 3.1). Both components focus on enhancing the speed of the DP pricing algorithm. These components play a decisive role. Deactivating just one of these components makes it impossible to compute the root-node relaxation in a reasonable time for many instances. As an immediate consequence, the number of optimal solutions dramatically decreases (down to 84/120) and the optimality gap soars (up to 25.02%) due to the number of incomplete root node calculations.

The remaining analyzes of Table 4 concern the B&P2 algorithm, based on route enumeration. In addition to the preprocessing strategies and DSSR, the DP algorithm used for route enumeration exploits a sophisticated succession of completion bounds (Section 3.3). Deactivating these bounds, as in Configuration E, hinders the performance of the route enumeration algorithm, which fails on 16/120 largest instances. Remark that any instance which is successfully enumerated is subsequently solved. Finally, deactivating the SRC separation is moderately detrimental: it leads to an 42% increase of CPU time and one additional open instance.

The results of the sensitivity analysis of the HGS, in Table 5, essentially highlight the importance of the SP-based intensification procedure. As visible in Configuration F, the solution quality of HGS significantly decreases (up to 0.82% average gap on the MUN instances) when this strategy is deactivated. HGS identifies the routes (columns) belonging to the optimal solutions within the available number of iterations, but due to the large number of ships, ship types, and the cargo-ship incompatibility matrix, combining these routes into a complete solution can be a challenging task. In contrast, the SP component is perfectly suited for this role, such that the combination of both techniques leads to a particularly effective matheuristic.

The impact of other components and parameter settings is less marked: increasing or decreasing the time limit of each SP model (Configurations H and I) does not impact the solution method, due to the fact that most SP models are solved within a few seconds. Moreover, deactivating our vehicle-dependent neighborhood restrictions (Configuration J) and population-diversity management strategies (Configuration K) leads to a small but significant decrease of solution quality.

6 Conclusions and Perspectives

As demonstrated in this paper, the literature on maritime logistics has attained a turning point where state-of-the-art exact algorithms can solve industrial and tramp ship routing optimization problems of a realistic scale. The B&P algorithm that we designed for this purpose capitalizes on multiple methodological elements to find a good balance between relaxation strength and pricing speed. As demonstrated by our sensitivity analyzes, its most critical method components concern the efficiency of the DP pricing and enumeration algorithms: our sophisticated preprocessing techniques and filters, the use of DSSR to reintroduce elementarity and the successive completion bounds are essential to solve large ITSRSP instances. These observations are in line with the works of Pecin et al. (2017), Sadykov et al. (2017) and the general research on MIP for VRPs which, for a large part, focuses on finding tighter relaxations and faster DP algorithms.

From the heuristic viewpoint, our experiments with a problem-tailored HGS show that the routes of optimal solutions can usually be quickly identified, but that crossover and local search methods are easily tricked into suboptimal route selections. Hybridizing the HGS with an SP solver fixes this issue and allows to attain near-optimal solutions within minutes.

Overall, the algorithms presented in this paper have contributed to push the limits of performance and problem tractability, but multiple avenues of research remain open concerning model accuracy. Indeed, despite its relevance for maritime transportation companies, the ITSRSP remains a mere simplification of reality. As highlighted in Christiansen and Fagerholt (2014), it does not consider load-dependent fuel consumption, possible load splitting and flexible cargo quantities, or the possibility of slow-steaming on selected route segments. Emission control areas and sea conditions (e.g. depth and currents) are also largely ignored. Last but not the least, considerable reductions of turnaround time may be achieved by jointly optimizing ship routing and port operations within integrated supply chains. Adapting state-of-the-art exact and heuristic algorithms to handle these complex attributes remain a significant challenge for the future.


  • Andersson et al. (2011) Andersson H, Christiansen M, Fagerholt K (2011) The maritime pickup and delivery problem with time windows and split loads. INFOR: Information Systems and Operational Research 49(2):79–91.
  • Baldacci et al. (2011) Baldacci R, Mingozzi A, Roberti R (2011) New route relaxation and pricing strategies for the vehicle routing problem. Operations Research 59(5):1269–1283.
  • Bertsimas et al. (2018) Bertsimas D, Jaillet P, Martin S (2018) Online vehicle routing: The edge of optimization in large-scale applications. Operations Research, Forthcoming .
  • Borthen et al. (2017) Borthen T, Loennechen H, Wang X, Fagerholt K, Vidal T (2017) A genetic search-based heuristic for a fleet size and periodic routing problem with application to offshore supply planning. EURO Journal on Transportation and Logistics .
  • Brønmo et al. (2007) Brønmo G, Christiansen M, Fagerholt K, Nygreen B (2007) A multi-start local search heuristic for ship scheduling—a computational study. Computers & Operations Research 34(3):900–917.
  • Brown et al. (1987) Brown GG, Graves GW, Ronen D (1987) Scheduling ocean transportation of crude oil. Management Science 33(3):335–346.
  • Bulhões et al. (2018) Bulhões T, Hà M, Martinelli R, Vidal T (2018) The vehicle routing problem with service level constraints. European Journal of Operational Research 265(2):544–558.
  • Christiansen and Fagerholt (2014) Christiansen M, Fagerholt K (2014) Ship routing and scheduling in industrial and tramp shipping. Toth P, Vigo D, eds., Vehicle Routing, 381–408 (Society for Industrial and Applied Mathematics).
  • Christiansen et al. (2007) Christiansen M, Fagerholt K, Nygreen B, Ronen D (2007) Maritime transportation. Barnhart C, Laporte G, eds., Transportation, 189–284 (Elsevier).
  • Christiansen et al. (2013) Christiansen M, Fagerholt K, Nygreen B, Ronen D (2013) Ship routing and scheduling in the new millennium. European Journal of Operational Research 228(3):467–483.
  • Christofides et al. (1981) Christofides N, Mingozzi A, Toth P (1981) Exact algorithms for the vehicle routing problem, based on spanning tree and shortest path relaxations. Mathematical Programming 20:255–282.
  • Contardo and Martinelli (2014) Contardo C, Martinelli R (2014) A new exact algorithm for the multi-depot vehicle routing problem under capacity and route length constraints. Discrete Optimization 12:129–146.
  • Desaulniers et al. (2008) Desaulniers G, Lessard F, Hadjar A (2008) Tabu search, partial elementarity, and generalized k-path inequalities for the vehicle routing problem with time windows. Transportation Science 42(3):387–404.
  • Duhamel et al. (2011) Duhamel C, Lacomme P, Prodhon C (2011) Efficient frameworks for greedy split and new depth first search split procedures for routing problems. Computers & Operations Research 38(4):723–739.
  • Dumas et al. (1991) Dumas Y, Desrosiers J, Soumis F (1991) The pickup and delivery problem with time windows. European Journal of Operational Research 54(1):7–22.
  • Fagerholt and Christiansen (2000a) Fagerholt K, Christiansen M (2000a) A combined ship scheduling and allocation problem. Journal of the Operational Research Society 51(7):834–842.
  • Fagerholt and Christiansen (2000b) Fagerholt K, Christiansen M (2000b) A travelling salesman problem with allocation, time window and precedence constraints — an application to ship scheduling. International Transactions in Operational Research 7(3):231–244.
  • Glover and Hao (2009) Glover F, Hao JK (2009) The case for strategic oscillation. Annals of Operations Research 183(1):163–173.
  • Gschwind et al. (2018) Gschwind T, Irnich S, Rothenbächer AK, Tilk C (2018) Bidirectional labeling in column-generation algorithms for pickup-and-delivery problems. European Journal of Operational Research 266(2):521–530.
  • Hemmati and Hvattum (2016) Hemmati A, Hvattum LM (2016) Evaluating the importance of randomization in adaptive large neighborhood search. International Transactions in Operational Research 24(5):929–942.
  • Hemmati et al. (2014) Hemmati A, Hvattum LM, Fagerholt K, Norstad I (2014) Benchmark suite for industrial and tramp ship routing and scheduling problems. INFOR: Information Systems and Operational Research 52(1):28–38.
  • Jepsen et al. (2008) Jepsen M, Petersen B, Spoorendonk S, Pisinger D (2008) Subset-row inequalities applied to the vehicle-routing problem with time windows. Operations Research 56(2):497–511.
  • Korsvik et al. (2009) Korsvik JE, Fagerholt K, Laporte G (2009) A tabu search heuristic for ship routing and scheduling. Journal of the Operational Research Society 61(4):594–603.
  • Korsvik et al. (2011) Korsvik JE, Fagerholt K, Laporte G (2011) A large neighbourhood search heuristic for ship routing and scheduling with split loads. Computers & Operations Research 38(2):474–483.
  • Martinelli et al. (2014) Martinelli R, Pecin D, Poggi M (2014) Efficient elementary and restricted non-elementary route pricing. European Journal of Operational Research 239(1):102–111.
  • Muter et al. (2010) Muter I, Birbil S, Sahin G (2010) Combination of metaheuristic and exact algorithms for solving set covering-type optimization problems. INFORMS Journal on Computing 22(4):603–619.
  • Nagata et al. (2010) Nagata Y, Bräysy O, Dullaert W (2010) A penalty-based edge assembly memetic algorithm for the vehicle routing problem with time windows. Computers & Operations Research 37(4):724–737.
  • Pecin et al. (2017) Pecin D, Contardo C, Desaulniers G, Uchoa E (2017) New enhancements for the exact solution of the vehicle routing problem with time windows. INFORMS Journal on Computing 29(3):489–502.
  • Prins (2004)

    Prins C (2004) A simple and effective evolutionary algorithm for the vehicle routing problem.

    Computers & Operations Research 31(12):1985–2002.
  • Prins et al. (2009) Prins C, Labadi N, Reghioui M (2009) Tour splitting algorithms for vehicle routing problems. International Journal of Production Research 47(2):507–535.
  • Righini and Salani (2008) Righini G, Salani M (2008) New dynamic programming algorithms for the resource constrained elementary shortest path problem. Networks 51(3):155–170.
  • Ronen (1983) Ronen D (1983) Cargo ships routing and scheduling: survey of models and problems. European Journal of Operational Research 12(2):119–126.
  • Ropke and Cordeau (2009) Ropke S, Cordeau JF (2009) Branch and cut and price for the pickup and delivery problem with time windows. Transportation Science 43(3):267–286.
  • Sadykov et al. (2017) Sadykov R, Uchoa E, Pessoa A (2017) A bucket graph based labeling algorithm with application to vehicle routing. Technical Report L-2017-7, Cadernos do LOGIS-UFF, Niterói, Brazil.
  • Stålhane et al. (2012) Stålhane M, Andersson H, Christiansen M, Cordeau JF, Desaulniers G (2012) A branch-price-and-cut method for a ship routing and scheduling problem with split loads. Computers & Operations Research 39(12):3361–3375.
  • Subramanian et al. (2013) Subramanian A, Uchoa E, Ochi LS (2013) A hybrid algorithm for a class of vehicle routing problems. Computers & Operations Research 40(10):2519–2531.
  • UNCTAD (2017) UNCTAD (2017) Review of maritime transport URL, accessed: 2018-09-06.
  • Velasco et al. (2009) Velasco N, Castagliola P, Dejax P, Guéret C, Prins C (2009) A memetic algorithm for a pick-up and delivery problem by helicopter, 173–190 (Springer Berlin Heidelberg).
  • Vidal (2016) Vidal T (2016) Technical note: Split algorithm in O(n) for the capacitated vehicle routing problem. Computers & Operations Research 69:40–47.
  • Vidal et al. (2012)

    Vidal T, Crainic TG, Gendreau M, Lahrichi N, Rei W (2012) A hybrid genetic algorithm for multidepot and periodic vehicle routing problems.

    Operations Research 60(3):611–624.
  • Vidal et al. (2013) Vidal T, Crainic TG, Gendreau M, Prins C (2013) A hybrid genetic algorithm with adaptive diversity management for a large class of vehicle routing problems with time-windows. Computers & Operations Research 40(1):475–489.
  • Vidal et al. (2014) Vidal T, Crainic TG, Gendreau M, Prins C (2014) A unified solution framework for multi-attribute vehicle routing problems. European Journal of Operational Research 234(3):658–673.
  • Vidal et al. (2015a) Vidal T, Crainic TG, Gendreau M, Prins C (2015a) Time-window relaxations in vehicle routing heuristics. Journal of Heuristics 21(3):329–358.
  • Vidal et al. (2015b) Vidal T, Crainic TG, Gendreau M, Prins C (2015b) Timing problems and algorithms: time decisions for sequences of activities. Networks 65(2):102–128.
  • Vidal et al. (2016) Vidal T, Maculan N, Ochi L, Penna P (2016) Large neighborhoods with implicit customer selection for vehicle routing problems with profits. Transportation Science 50(2):720–734.
  • Vilhelmsen et al. (2014) Vilhelmsen C, Lusby R, Larsen J (2014) Tramp ship routing and scheduling with integrated bunker optimization. EURO Journal on Transportation and Logistics 3(2):143–175.
  • Wilson (2018) Wilson (2018) URL, accessed: 2018-09-06.
  • Zachariadis and Kiranoudis (2010) Zachariadis E, Kiranoudis C (2010) A strategy for reducing the computational complexity of local search-based methods for the vehicle routing problem. Computers & Operations Research 37(12):2089–2105.

Appendix – Detailed Results

Tables 6 to 9 give the detailed results of B&P2 and HGS2 (with the set partitioning component) on each individual instance. For the B&P2, Columns “Opt” and “T” represent the optimal value and the total CPU time (in minutes). For the HGS2, Columns “Best”, “Avg” and “Worst” represent the best, average and worst cost found over runs, column “T” gives the total CPU time in minutes, and column “T*” gives the attainment time (in minutes) needed to reach the final solution.

Instance B&P2 HGS2
Opt T Best Avg Worst T T*
SHORTSEA_MUN_C7_V3_HE_1 1476444 0.00 1476444 1476444 1476444 0.02 0.00
SHORTSEA_MUN_C7_V3_HE_2 1134176 0.00 1134176 1134176 1134176 0.02 0.00
SHORTSEA_MUN_C7_V3_HE_3 1196466 0.00 1196466 1196466 1196466 0.02 0.00
SHORTSEA_MUN_C7_V3_HE_4 1256139 0.00 1256139 1256139 1256139 0.02 0.00
SHORTSEA_MUN_C7_V3_HE_5 1160394 0.00 1160394 1160394 1160394 0.02 0.00
SHORTSEA_MUN_C10_V3_HE_1 2083965 0.00 2083965 2083965 2083965 0.04 0.00
SHORTSEA_MUN_C10_V3_HE_2 2012364 0.00 2012364 2012364 2012364 0.04