1 Introduction
International trade depends heavily on ship transportation, as it is the only costeffective 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 (pickupanddelivery pairs) has been made available to the academic community. The authors also presented a compact mathematical formulation, and solved it with a branchandcut 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 subsetrow 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 highquality 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 pickupanddelivery problems as these problems require structurallydifferent localsearch 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 problemtailored crossover, LS operators and vehicledependent 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 stateoftheart 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 multiship PDPTW. A maritime PDPTW with split loads was studied by Andersson et al. (2011). The authors proposed two pathflow 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 multistart 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 multiperiod 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 prizecollecting 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 pickupanddelivery problems, which require structurally different neighborhood searches and proper precedence and pairing between the pickups and deliveries in the crossover and split operators.
3 BranchandPrice 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.(1)  
(2)  
(3)  
(4)  
(5) 
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 (1–5) 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).
(6)  
(7) 
Since the dual costs are originally associated with cargoes (p–d pairs), we opted to associate all the dual costs with the outarcs 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 timewindow constraints. A similar property has been exploited in Bertsimas et al. (2018) to optimize taxifleets 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:
(8) 
If an extension is allowed, it generates the new label presented in Equation (9):
(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.
(10)  
Note that implies that . However, as discussed in Ropke and Cordeau (2009), a subsetbased 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 negativecost route was last obtained. When the pricing algorithm fails to generate a negative reducedcost 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 Branchandbound
The CG presented in the previous section produces strong lower bounds for the ITSRSP. To obtain integer optimal solutions, we embed it into a branchandbound algorithm to form a BranchandPrice (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 righthand 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 twophase 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 nonoptimal 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):
(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 3SRC to improve the value of the linear relaxations (Jepsen et al. 2008). For the ITSRSP, the 3SRC 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 3SRC.
(12) 
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 3SRC can be first separated at the root node to improve the linear relaxation and reduce the number of routes resulting from the enumeration. Moreover, 3SRC 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 industrialsize 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 problemtailored search operators and an SPbased 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).
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 highquality 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 highquality 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 timewindow constraints, we use the “timewarp” 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 shipcargo 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 startofservice time at the th node can be defined as:
(13) 
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:
(18) 
where , , and are the respective penalty coefficients for peakload, timewindow, and incompatibilityconstraint 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 pseudopolynomial 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
(19) 
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
(20) 
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 onepoint 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 pickupanddelivery problems. The moves are evaluated in random order, and any improving move is directly applied (firstimprovement 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 intraroute modifications, while and may be used for both intra and interroute modifications, and they can therefore change cargoship 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  
Load  
Peak load  
Time warp use  
Earliest possible completion time  
Latest feasible starting time  
Duration  
Incompatibilities  
Auxiliary computations  
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 lastmodified time of a route and the lastevaluated time for each move. By comparing these values, one can decide whether or not to reevaluate 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.
Vehicledependent 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 :
(21)  
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 timewindow, 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 highquality 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 1–5) 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 practicalsize 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., vehicledependent neighborhood restrictions and setpartitioning 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 i73960X CPU and 64 GB of RAM.
We rely on the benchmark instances suite for the ITSRSP based on reallife scenarios, presented in Hemmati et al. (2014) and currently available at http://home.himolde.no/~hvattum/benchmarks/. 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 mixedload instances have up to cargoes and ships, whereas the fullload instances have up to cargoes and ships. So far, instances remain open. In the fullload 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 mixedload 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 longdistance 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 branchandprice without route enumeration (B&P_{1}), using heuristic strong branching with 50 candidates at each iteration, as well as the same algorithm with route enumeration and 3SRC using inspection pricing (B&P_{2}). 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&P_{1}: “Gap_{0}” and “T_{0}” represent the percentage gap and the time in minutes for the root node. “Gap_{F}” and “T_{F}” are the final percentage gap and time, “N_{F}” 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&P_{2} (with route enumeration and inspection pricing): the columns “T_{E}” and “R_{E}” are the time in minutes for the route enumeration and the number of routes found. “Gap_{0}”, “T_{0}”, and “Cuts_{0}” are the percentage gap, the overall time and the number of 3SRC separated at the root node. “R_{F}”, “Gap_{F}”, “T_{F}”, “Cuts_{F}”, and “N_{F}” are the final number of routes, the percentage gap, the overall time, the 3SRC 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)  BranchandPrice (B&P_{1})  BranchandPrice + Enumeration + SRCs (B&P_{2})  
Instances  T  Gap  Opt  Gap_{0}  T_{0}  Gap_{F}  T_{F}  N_{F}  Opt  T_{E}  R_{E}  Gap_{0}  T_{0}  Cuts_{0}  R_{F}  Gap_{F}  T_{F}  Cuts_{F}  N_{F}  Opt  
SSMUN 
C7V3  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 
C10V3  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  
C15V4  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  
C18V5  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  
C22V6  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  
C23V13  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  
C30V6  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  
C35V7  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  
C60V13  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  
C80V20  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  
C100V30  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  
C130V40  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  
SSFUN 
C8V3  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 
C11V4  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  
C13V5  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  
C16V6  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  
C17V13  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  
C20V6  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  
C25V7  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  
C35V13  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  
C50V20  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  
C70V30  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  
C90V40  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  
C100V50  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  
DSMUN 
C7V3  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 
C10V3  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  
C15V4  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  
C18V5  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  
C22V6  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  
C23V13  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  
C30V6  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  
C35V7  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  
C60V13  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  
C80V20  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  
C100V30  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  
C130V40  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  
DSFUN 
C8V3  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 
C11V4  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  
C13V5  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  
C16V6  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  
C17V13  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  
C20V6  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  
C25V7  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  
C35V13  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  
C50V20  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  
C70V30  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  
C90V40  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  
C100V50  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 
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&P_{1} 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 mixedcargo instances, the average time needed to solve each pricing subproblem (ratio T_{F}/N_{F}) increases with the number of cargoes, whereas the lower bounds tend to deteriorate, leading to larger branchandbound trees. To go further, one can concentrate on improving the pricing problem solution or the lower bounds. For the largest open instances, the rootnode gap seems small enough to allow a complete route enumeration via sophisticated DP algorithms (Section 3.3). We therefore derived B&P_{2} 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&P_{2}. This approach outperforms the standard B&P_{1} 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 DSMUNC130V40HE2). Having enumerated the routes allows to efficiently separate the SRCs and decrease the rootnode 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 branchandbound tree (from 10.6 to 4.2 on average). Using this method, all available instances but one (DSMUNC130V40HE2) 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 displays the number of instances solved by B&P_{1} and B&P_{2} as a function of the CPU time limit. B&P_{2} 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&P_{1} 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 (HGS_{1}) and with SP (HGS_{2}).
We used the same parameters as Vidal et al. (2013) to avoid any problemspecific 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 fullload 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 HGS_{1} and HGS_{2}. 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 https://w1.cirrelt.ca/~vidalt/en/VRPresources.html.
Hemmati et al. (2014)  HGS_{1}  HGS_{2}  
Instances  Gap  T  Gap  T  T*  Gap  T  T*  
SSMUN 
C7V3  0.00  0.03  0.00  0.02  0.00  0.00  0.02  0.00 
C10V3  0.00  0.04  0.00  0.04  0.00  0.00  0.04  0.00  
C15V4  0.58  0.09  0.00  0.08  0.00  0.00  0.08  0.00  
C18V5  0.51  0.13  0.07  0.12  0.02  0.02  0.13  0.03  
C22V6  1.82  0.19  0.05  0.19  0.05  0.00  0.18  0.04  
C23V13  0.58  0.25  0.13  0.25  0.09  0.00  0.23  0.07  
C30V6  1.60  0.38  0.17  0.36  0.17  0.00  0.33  0.12  
C35V7  1.92  0.54  0.81  0.51  0.26  0.03  0.54  0.27  
C60V13  1.69  2.01  1.03  2.44  1.81  0.14  2.25  1.37  
C80V20  2.45  4.13  1.79  5.21  4.09  0.03  3.84  2.44  
C100V30  2.68  7.77  1.67  9.97  8.01  0.01  5.31  3.02  
C130V40  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  
SSFUN 
C8V3  0.00  0.03  0.00  0.01  0.00  0.00  0.01  0.00 
C11V4  0.13  0.05  0.00  0.02  0.00  0.00  0.02  0.00  
C13V5  0.07  0.07  0.00  0.02  0.00  0.00  0.02  0.00  
C16V6  0.05  0.10  0.00  0.03  0.00  0.00  0.03  0.00  
C17V13  0.01  0.14  0.00  0.04  0.00  0.00  0.05  0.00  
C20V6  0.14  0.18  0.00  0.04  0.00  0.00  0.05  0.00  
C25V7  0.22  0.27  0.00  0.07  0.02  0.00  0.07  0.01  
C35V13  0.29  0.60  0.01  0.19  0.08  0.00  0.17  0.05  
C50V20  0.52  1.38  0.22  0.67  0.43  0.00  0.43  0.15  
C70V30  1.29  3.51  0.68  1.86  1.28  0.00  1.09  0.42  
C90V40  1.45  6.98  0.70  3.90  2.83  0.00  1.93  0.76  
C100V50  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  
DSMUN 
C7V3  0.00  0.03  0.00  0.02  0.00  0.00  0.02  0.00 
C10V3  0.01  0.04  0.00  0.04  0.00  0.00  0.04  0.00  
C15V4  1.26  0.08  0.00  0.08  0.01  0.00  0.08  0.01  
C18V5  0.47  0.13  0.00  0.12  0.01  0.00  0.12  0.01  
C22V6  2.18  0.19  0.01  0.18  0.05  0.00  0.18  0.05  
C23V13  0.12  0.24  0.02  0.19  0.05  0.00  0.20  0.04  
C30V6  1.04  0.37  0.22  0.36  0.17  0.03  0.32  0.12  
C35V7  1.08  0.51  0.14  0.53  0.28  0.00  0.44  0.17  
C60V13  3.74  1.92  1.41  2.58  1.96  0.03  1.81  1.06  
C80V20  3.10  4.26  1.61  5.69  4.58  0.04  3.56  2.20  
C100V30  3.69  8.00  3.27  9.22  7.27  0.01  9.87  6.84  
C130V40  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  
DSFUN 
C8V3  0.00  0.03  0.00  0.01  0.00  0.00  0.01  0.00 
C11V4  0.00  0.05  0.00  0.02  0.00  0.00  0.02  0.00  
C13V5  0.00  0.06  0.00  0.02  0.00  0.00  0.04  0.00  
C16V6  0.03  0.10  0.00  0.03  0.00  0.00  0.04  0.00  
C17V13  0.00  0.13  0.00  0.05  0.00  0.00  0.05  0.00  
C20V6  0.01  0.16  0.00  0.04  0.00  0.00  0.05  0.00  
C25V7  0.41  0.26  0.00  0.07  0.01  0.00  0.07  0.01  
C35V13  1.03  0.59  0.01  0.26  0.14  0.00  0.21  0.08  
C50V20  0.61  1.41  0.23  0.71  0.46  0.00  0.45  0.17  
C70V30  0.59  3.55  0.31  2.02  1.43  0.00  1.02  0.38  
C90V40  1.10  7.01  0.44  4.27  3.16  0.00  2.20  1.00  
C100V50  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 
As visible in the results of Table 3, both versions of the HGS largely outperform previous algorithms. HGS_{2} obtains nearoptimal solutions within an average gap of , compared to for the ALNS. The fullload instances are generally easier to solve than the mixedload instances since the problem simplifies. Moreover, the SPbased 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 (HGS_{1}) to min (HGS_{2}) 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 HGS_{1} (T* ) with that of HGS_{2} (). In a followup 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: HGS_{2} 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&P_{1}, B&P_{2} and HGS_{2}) 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 4–5. 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.
FUN  MUN  
Gap  T  Root  Opt  Gap  T  Root  Opt  
Standard (B&P_{1})  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&P_{2})  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 SubsetRow Cuts  0.00  0.01  120  120  0.00  2.47  120  118 
FUN  MUN  
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 (HGS_{1})  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 
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 rootnode 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&P_{2} 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 SPbased 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 cargoship 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 vehicledependent neighborhood restrictions (Configuration J) and populationdiversity 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 stateoftheart 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 problemtailored 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 nearoptimal 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 loaddependent fuel consumption, possible load splitting and flexible cargo quantities, or the possibility of slowsteaming 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 stateoftheart exact and heuristic algorithms to handle these complex attributes remain a significant challenge for the future.
References
 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 largescale applications. Operations Research, Forthcoming .
 Borthen et al. (2017) Borthen T, Loennechen H, Wang X, Fagerholt K, Vidal T (2017) A genetic searchbased 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 multistart 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 multidepot 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 kpath 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 columngeneration algorithms for pickupanddelivery 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) Subsetrow inequalities applied to the vehiclerouting 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 nonelementary 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 coveringtype optimization problems. INFORMS Journal on Computing 22(4):603–619.
 Nagata et al. (2010) Nagata Y, Bräysy O, Dullaert W (2010) A penaltybased 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 L20177, Cadernos do LOGISUFF, Niterói, Brazil.
 Stålhane et al. (2012) Stålhane M, Andersson H, Christiansen M, Cordeau JF, Desaulniers G (2012) A branchpriceandcut 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 http://unctad.org/en/PublicationsLibrary/rmt2017_en.pdf, accessed: 20180906.
 Velasco et al. (2009) Velasco N, Castagliola P, Dejax P, Guéret C, Prins C (2009) A memetic algorithm for a pickup 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 timewindows. 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 multiattribute 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) Timewindow 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 https://www.wilsonship.no/en, accessed: 20180906.
 Zachariadis and Kiranoudis (2010) Zachariadis E, Kiranoudis C (2010) A strategy for reducing the computational complexity of local searchbased 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&P_{2} and HGS_{2} (with the set partitioning component) on each individual instance. For the B&P_{2}, Columns “Opt” and “T” represent the optimal value and the total CPU time (in minutes). For the HGS_{2}, 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&P_{2}  HGS_{2}  

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 