Stabilized Branch-Price-and-Cut for the Commodity-constrained Split Delivery Vehicle Routing Problem

In the commodity-constrained split delivery vehicle routing problem (C-SDVRP), customer demands are composed of sets of different commodities. The C-SDVRP asks for a minimum-distance set of vehicle routes such that all customer demands are met and vehicle capacities are respected. Moreover, whenever a commodity is delivered by a vehicle to a customer, the entire amount requested by this customer must be provided. Different commodities demanded by one customer, however, can be delivered by different vehicles. Thus, the C-SDVRP is a relaxation of the capacitated vehicle routing problem and a restriction of the split delivery vehicle routing problem. For its exact solution, we propose a branch-price-and-cut algorithm that employs and tailors stabilization techniques that have been successfully applied to several cutting and packing problems. More precisely, we make use of (deep) dual-optimal inequalities which are particularly suited to reduce the negative effects caused by the inherent symmetry of C-SDVRP instances. One main issues here is the interaction of branching and cutting decisions and the different classes of dual inequalities. Extensive computational tests on existing and extended benchmark instances show that all stabilized variants of our branch-price-and-cut are clearly superior to the non-stabilized version. On the existing benchmark, we are significantly faster than the state-of-the-art algorithm and provide several new optima for instances with up to 60 customers and 180 tasks. Lower bounds are reported for all tested instances with up to 80 customers and 480 tasks, improving the bounds for all unsolved instances and providing first lower bounds for several instances.


Introduction
The commodity-constrained split delivery vehicle routing problem (C-SDVRP) has been introduced by Archetti et al. (2016) to implement one of the possible distribution policies applicable when different commodities have to be distributed to customers. The underlying idea was to study the impact on variable routing costs from adopting different distribution policies, i.e., from using vehicles dedicated to a single commodity compared with using flexible vehicles capable of carrying any set of commodities, and from allowing or forbidding split deliveries of individual commodities. The policy associated with the C-SDVRP allows the homogeneous vehicles to visit a customer several times, but when a commodity is delivered by a vehicle to a customer, the entire amount requested by the customer has to be provided. No compatibility restrictions exist so that vehicles can carry any subset of commodities. Thus, the number of commodies requested by a customer is an upper bound on the number of visits to this customer. Overall, the C-SDVRP is a relaxation of the classical capacitated vehicle routing problem (CVRP) and a restriction of the split delivery vehicle routing problem, see the book by Toth and Vigo (2014), in particular chapters (Semet et al., 2014;Poggi and Uchoa, 2014;Irnich et al., 2014).
In order to evaluate the proposed distribution policies, Archetti et al. (2016) carried out worst-case as well as experimental analyses. In particular, they devised exact and heuristic algorithms. The tests were run on 64 small instances, with 15 customers and up to three commodities, 80 mid-size instances with 20, 40, 60, or 80 customers and up to three commodities, and large instances with 100 customers. One of the main conclusions drawn is that the delivery policy associated with the C-SDVRP is almost the best option w.r.t. to costs and, at the same time, it is likely more acceptable to customers than allowing splitting all the deliveries.
A straightforward approach for modeling and solving the C-SDVRP is to reduce it to a standard CVRP, as done in (Archetti et al., 2016) to solve the C-SDVRP heuristically. The authors duplicated each customer vertex as many times as the number of commodities requested by the customer, associating with each duplicated vertex the customer's demand for the corresponding commodity. Nevertheless, in order to limit the size of the model and mitigate symmetry issues, the branch-and-cut proposed by the authors is based on a direct model for the C-SDVRP. Archetti et al. (2015) focused on the exact solution of the problem. While the branch-and-cut algorithm proposed in Archetti et al. (2016) was able to solve only 25 out of the 64 small instances, the branch-priceand-cut (BPC) algorithm proposed in Archetti et al. (2015) was able to solve all the small instances within the same time limit, and other instances with up to 40 customers and three commodities.
Under the name discrete split delivery vehicle routing problem, the C-SDVRP had already been introduced in the literature by Nakao and Nagamochi (2007) to model a specific practical routing problem as a variant of the SDVRP. The authors proposed a fast heuristic based on dynamic programming to solve real-world instances.
Another related problem is the discrete split delivery vehicle routing problem with time windows (DSD-VRPTW) proposed by Salani and Vacca (2011), where the demand of a customer consists of several items. The authors assumed that demand can be split in orders, i.e., feasible combinations of items, that each vehicle can serve at most one order per customer and that service time at a customer's location depends on the delivered combination of items. The problem is solved by means of a branch-and-price algorithm. The DSDVRPTW can indirectly model the C-SDVRP by defining the set of orders to associate with a customer as the set of all possibles subsets of commodities to deliver to him.
In this paper, we focus on the exact solution of the C-SDVRP via column-generation techniques (Lübbecke and Desrosiers, 2005;Desaulniers et al., 2005). We enhance the above-mentioned BPC algorithm of Archetti et al. (2015) in several aspects. First, we integrate powerful techniques such as implicit bidirectional labeling for solving the column-generation subproblem (Righini and Salani, 2006;Bode and Irnich, 2012) and additional cuts to strengthen the linear relaxation of the master program (Jepsen et al., 2008). Second and more importantly, we tailor recent stabilization techniques originally suggested for bin-packing and vertex-coloring problems to the C-SDVRP (Ben Amor et al., 2006;Gschwind and Irnich, 2016). Here, the stabilization of the column-generation process is achieved by the addition of (deep) dual-optimal inequalities, i.e., inequalities known to hold for (some) optimal dual solutions. We show that stabilization with dual inequalities is particularly suited to reduce the negative effects caused by the inherent symmetry of C-SDVRP instances. Indeed, when the number of tasks, i.e., the number of commodities to deliver to all customers, is kept constant, the most symmetric instances are those that have the largest number of commodities per customer, and these instances become less difficult to solve when dual inequalities are added. From a methodological point of view, the challenge is to correctly handle the impact that branching and cutting have on the solution of the subproblem and the validity of different classes of dual inequalities used for stabilization. The paper puts specific emphasis on these interdependencies.
The remainder of the paper is structured as follows. The C-SDVRP is formally defined in Section 2. Details on the new BPC algorithm including the generation of route variables, stabilization techniques, strengthening of the linear relaxation by adding valid inequalities, and branching are provided in Section 3. The results of our computational studies are provided in Section 4, before the paper closes with conclusions drawn in Section 5.

Problem Definition
Let K = {1, . . . , κ} be the set of commodities that has to be distributed from the depot to the customers. Let N = {1, . . . , n} be the set of customers. The demand of commodity k ∈ K to deliver to customer i ∈ N is denoted by d ik . The set K i = {k ∈ K : d ik > 0} comprises the commodities to be delivered to customer i ∈ N . The number of delivery tasks is therefore m = i∈N |K i |. Customers are served by means of a fleet of F homogeneous vehicles, each of which with a capacity of Q. Vehicles are flexible and can deliver any subset of commodities. Each customer may be visited more than once. When a commodity is delivered by a vehicle to a customer, the entire amount of the commodity requested by the customer has to be provided.
Let G = (V, A) be a directed graph with vertex set V = N ∪ {0, n + 1}, and arc set A. Vertices 0 and n + 1 represent the depot where vehicle routes start and end, respectively. Each arc (i, j) ∈ A, with i = n + 1 and j = 0, represents the possibility for the vehicles to travel from the location of i to the location of j, and it is associated with nonnegative variable routing costs c ij . In the following, the term route is used for the combination of (i) an elementary 0-(n + 1)-path representing the sequence in which customers are visited by the vehicle, (ii) and a selection of commodities delivered to each visited customer i ∈ N , i.e., a non-empty subset S i ⊆ K i . Let Ω be the set of all feasible routes. For each r ∈ Ω, the routing costs are c r = (i,j)∈A(r) c ij , and route r is feasible if the total demand delivered does not exceed the vehicle capacity, i.e., i∈N (r) k∈Si d ik ≤ Q, where N (r) ⊆ N is the set of customers visited along route r. Similarly, we define A(r) ⊆ A as the set of arcs defining the path associated with route r.
The C-SDVRP is the problem of determining a set of least-cost feasible routes to be associated with the vehicles such that all customers' demands are met.

Branch-Price-and-Cut
In order to solve the C-SDVRP, we adopted the same set-partitioning formulation as proposed by Archetti et al. (2015). The formulation has, for each route r ∈ Ω, one binary variable λ r that assumes a value equal to 1 if route r is performed, and 0 otherwise. Moreover, the non-negative integer variable φ models the number of routes performed (i.e., number of vehicles used), and x ij and z i are non-negative integer variables representing the number of times that an arc (i, j) ∈ A is traversed and a customer i ∈ N is visited, respectively. The formulation is as follows: The objective function (1a) calls for the minimization of the total variable routing costs. Constraints (1b) ensure that all m customers' demands are met (i.e., tasks are fulfilled), where a r ik is a binary coefficient indicating that commodity k ∈ K is delivered to customer i ∈ N . Constraints (1c) and (1d) limit the number of vehicles to use, and (1e) are domain definition constraints for variables λ r .
The number of times a customer i ∈ N is visited and an arc (i, j) ∈ A is traversed is limited due to constraints (1f)-(1i). For this purpose, the binary coefficients e r i and b r ij are equal to 1 if customer i ∈ N is visited and arc (i, j) ∈ A is traversed by r, respectively, and 0 otherwise. Note that constraints (1f)-(1i) are redundant. They may only be added when violated by a fractional solution to the linear relaxation of (1) and for branching, see Section 3.4.
For many VRP variants, BPC is the most successful and leading exact solution approach. Hence, we solve formulation (1) with BPC. The starting point is a restricted master program (RMP) which is the linear relaxation of formulation (1) defined over a small subset Ω ⊂ Ω of the route variables λ r . Column generation alternates between the optimization of the RMP and the solution of the column-generation subproblem that generates negative reduced cost route variables λ r if they exist . Archetti et al. (2015) have shown that this subproblem can be formulated and solved as a variant of the shortest-path problem with resouce constraints (SPPRC, Irnich and Desaulniers, 2005). If no negative reduced cost routes exist, a solution to the linear relaxation of (1) is found. The corresponding lower bound can be strengthened by the addition of valid inequalities. Finally, branching is required to ensure integer solutions.

Column Generation
In this section, we sketch the SPPRC-based solution approach of Archetti et al. (2015) for the columngeneration subproblem. In addition, we refine parts of this approach by introducing another representation for the knapsack problem associated with each customer and further exploiting the inherent symmetry of the subproblem.
Recall that an instance of the column-generation subproblem is defined by the dual prices π ik , i ∈ N, k ∈ K i associated with constraints (1b), σ associated with (1c), µ i , i ∈ N associated with (1f), and ρ ij , (i, j) ∈ A associated with (1h). The reduced cost of a route r that delivers commodities Archetti et al. (2015) have shown that the combined problem of determining the path and the commodities to deliver can be formulated as an SPPRC over a multi-graph. Slightly differing from their work, we define the multi-graph as an undirected graph (but otherwise identical) as follows: The vertices of the multi-graph comprise two copies 0 and 0 of the depot as well as two copies i and i for each customer i ∈ N . Moreover, for all original arcs (i, j) ∈ A, there are two edges (i , j ) and (i , j ) in the multi-graph for modeling the movement of the vehicle. Finally, for all non-empty subsets S i ⊆ K i , i ∈ N , there are parallel edges between i and i , denoted as (i , i ) Si , that represent the respective deliveries made to customer i. An example of the multi-graph for an instance with three customers is depicted in Figure 1.  Any 0 -0 -path in the multi-graph that alternates between vertices in V := {0 } ∪ {n : n ∈ N } and V = {0 } ∪ {n : n ∈ N } corresponds with a route, and vice versa. Defining demands of zero on the edges (i , j ) and (i , j ) as well as demands on edges (i , i ) Si , such a 0 -0 -path represents a feasible route if the demands accumulated over its edges do not exceed Q. Obviously, edges (i , i ) Si with a demand exceeding Q can be eliminated from the multi-graph. With the following definitions, the multi-graph has a completely symmetric reduced-cost structure if the original instance is symmetric (c ij = c ji ): where we additionally define µ 0 := σ. Note that the benchmark instances for the C-SDVRP are all symmetric, see Section 4.1. The solution approach of Archetti et al. (2015) proceeds in two phases as follows: Phase 1: Pre-computation of all Pareto-optimal pairs (c Si i ,i , d Si i ,i ) for each customer i ∈ N ; Phase 2: Solution of an SPPRC defined over a reduced multi-graph containing only Pareto-optimal edges between i and i for all customer i ∈ N .
Pareto-optimal Deliveries. In principle, Pareto-optimal commodity combinations S i ⊆ K i could be determined via enumeration, since the number of commodities per customer usually does not exceed ten. However, we present an SPPRC-based technique that is also beneficial for separating violated dual-optimal inequalities (see Section 3.2), when non-robust cuts are added to (1) (see Section 3.3), and when branching constraints must be taken into account (see Section 3.4). For each customer i ∈ N , we first order the commodities to deliver at i arbitrarily so that K i = {k 1 , k 2 , . . . , k |Ki| } (for the sake of convenience, we omit an extra index i). We define the customer network as a digraph (N i , A i ) with |K i | + 1 vertices given by N i = {0, 1, 2, . . . , |K i |}. For each index j ∈ N i , j = 0, there are two arcs (j − 1, j) 0 and (j − 1, j) 1 between j − 1 and j. The arc (j − 1, j) 0 models that the jth commodity k j is not selected, while arc (j − 1, j) 1 models the selection of k j . Accordingly, arc (j − 1, j) 0 has zero weight and zero profit, and arc (j − 1, j) 1 has weight d ikj and profit π ikj . Figure 2 displays an example of such a customer network.
is in one-to-one correspondence to a (possibly empty) subset S i of K i . If the weight of the 0-|K i |-path is positive and does not exceed Q, the corresponding subset S i represents a feasible delivery. Hence, solving an SPPRC with resources profit (to maximize) and weight (constrained by Q) produces all Pareto-optimal sets S i . Note that dominance between SPPRC labels must be slightly modified: as S i = ∅ is not allowed (at least one commodity must be delivered), labels (0, 0) must not dominate and eliminate other labels.
SPPRC over the Multi-Graph. The reduced multi-graph changes from one column-generation iteration to the next because all edges (i , i ) Si whose subsets S i are not Pareto-optimal are temporarily removed. The actual SPPRC is then solved using the following resources: (i) accumulated reduced cost according to (3); (ii) accumulated demand according to (2); and (iii) visit indicators for each customer i ∈ N . In monodirectional forward labeling, labels at a vertex i are only propagated towards i (same i), while labels at i ∈ V are only propagated to vertices j ∈ V (with different j = i).
All resources are initially set to 0. While reduced costs are unconstrained, the accumulated demand is bounded from above by Q. The visit indicator of customer i ∈ N is increased when one of the arcs (i , i ) Si is traversed, and it is bounded from above by 1.
As this elementary SPPRC is NP-hard in the strong sense, relaxations are typically employed. We use the ng-path relaxation of Baldacci et al. (2011), in which a cycle over a task i ∈ N is possible if i is not in the neighborhood of a vertex of the cycle. Neighborhoods must be pre-defined, and in most implementations reported in the literature a global number controls the maximum size of all neighborhoods. The larger the neighborhoods, the less cycles are possible, but on average the computational difficulty increases. Good tradeoffs were often obtained with neighborhoods of size between 5 and 20.
Bidirectional labeling for SPPRCs was coined by Righini and Salani (2006) in order to mitigate the explosion of labels typically observed when partial paths grow longer. In bidirectional labeling, both forward partial paths and backward partial paths are created, but processed only up to a so-called half-way point. In VRP variants with only capacity constraints, the termination criterion is often that the accumulated demand has reached or exceeded half of the vehicle's capacity, i.e., the half-way point is Q/2. When forward and backward labeling terminates, suitable forward and backward labels must be merged to obtain a complete feasible route. Several subsequent works have shown that bounded bidirectional labeling algorithms are usually superior to their monodirectional counterparts.
If the SPPRC is completely symmetric, e.g., guaranteed by definitions (3), the computational effort can be further reduced by using an implicit bidirectional approach. Note that forward and backward propagation produce essentially identical partial paths, since forward partial paths start at 0 and always traverse edges (i , i ) Si in the given direction, while backward partial paths start at 0 and traverse edges (i , i ) Si in backward direction (from i to i ). Hence, swapping all i and i , respectively, allows to produce backward partial paths out of forward partial paths, and vice versa. Thus, label propagation in one direction can be omitted, so that, e.g., only forward partial paths are combined in the merge procedure. This implicit bidirectional techniques has already been successfully implemented and applied in (Bode and Irnich, 2012;Goeke et al., 2017).

Stabilization and Dual Inequalities
In this section, we describe the key innovation of this research: the use of dual inequalities for the stabilization of the column-generation process. For a general introduction to the topic we refer to (Ben Amor et al., 2006;Gschwind and Irnich, 2016).
Dual Inequalities. We start with some basic definitions. For the C-SDVRP, any linear inequality in the dual variables is a dual inequality (DI). The DIs that we will consider in the following refer only to the dual variables (π ik ) for i ∈ N, k ∈ K i associated with the m partitioning constraints (1b). Let D * be the dualoptimal space, i.e., the set of optimal solutions of the dual model to the linear relaxation of (1). Following Ben Amor et al. (2006), a DI t π ≤ t (with t ∈ Z m and t ∈ Z) is a dual-optimal inequality (DOI) if D * ⊆ {π : t π ≤ t}. A set of DIs comprises deep dual-optimal inequalities (DDOIs) if at least one π * ∈ D * fulfills all inequalities. Hence, DOIs are always DDOIs, but the reverse is not necessarily true.
A trivial but already effective stabilization happens when the partitioning constraints (1b) are replaced by covering constraints. This is a valid replacement (preserving the linear relaxation bound) if the cost matrix fulfills the triangle inequality. We assume that the triangle inequality is fulfilled for the C-SDVRP instances. Note that this replacement is equivalent to the DIs π ik ≥ 0 for all i ∈ N, k ∈ K i (see also Gschwind and Irnich, 2016, Proposition 4).
The constraint matrix of (1) does also possess the row replacement property (see Gschwind and Irnich, 2016, Section 3.3, Definition 2, and Property 5). Indeed, for a fixed customer i ∈ N and two different commodities k, k ∈ K i with demands fulfilling d ik ≤ d ik , a route serving customer i and delivering k and not k can always be feasibly replaced by an otherwise identical route that delivers k and not k. For the sake of simplicity, we assume that all commodities delivered to customer i ∈ N are sorted by non-decreasing demand, i.e., holds for all i ∈ N . With this sorting, π i,kp − π i,kq ≤ 0 for p ≤ q are DDOIs, in the following denoted as pair inequalities. They result from the subset of |K i | − 1 so-called ranking inequalities which are as powerful as the pair inequalities. Hence, the |K i | − 1 ranking inequalities can be used instead of all |Ki| 2 pair inequalities. Note also that if two commodities k, k of a customer i ∈ N have identical demand, the two pair inequalities result in the equality π ikp = π ikq , which is equivalent to replacing the two corresponding covering constraints (1b) by an aggregated ≥ 2 constraint. For more details on the relationship between constraint aggregation and dual equalities we refer to (Gschwind and Irnich, 2016, Property 6). We do not apply aggregation here, because typical C-SDVRP instances have different demands for all commodities delivered to a customer, see Section 4.
As known from bin packing, the pair inequalities can be generalized in the following way: consider a customer i ∈ N , one of its commodities k ∈ K i and a subset S ⊆ K i with i.e., the amount to deliver of k is not smaller than the amount for the commodities k ∈ S. Then, is a so-called subset inequality (SI). In general, SIs are no DDOIs for the C-SDVRP. However, from Gschwind and Irnich (2016, Theorem 2) it follows that, if d i,k + min k∈S d i,k > Q, then pair or subset inequalities (5) are DOIs. (5), defined by the triple (i, k , S) with i ∈ N , k ∈ K i , and S ⊂ K i . Let Θ be the set of all active DIs (at some point in time). For the primal model (1), each DI for (i, k , S) ∈ Θ is implemented as a DI column with entry −1 for the row (i, k ) and entries +1 for all rows (i, k), k ∈ S of the constraints (1b). With non-negative variables y i,k,S ≥ 0 for (i, k, S) ∈ Θ, constraints (1b) are replaced by

Static versus Dynamic Addition of Dual Inequalities. The most general form of DIs we use is the SIs
Moreover, there are several possible strategies to add DIs. We describe three straightforward strategies now: Static: Right from the initialization of the master program, the DI columns are added to (1). The advantage is that the stabilization effect happens from the first iteration on. However, more columns in the RMP make the re-optimization slower. Dynamic: Only DIs that are violated may be added to the master problem (1).
The advantage is a generally smaller RMP with a faster re-optimization, but a less stabilized columngeneration process. Moreover, the identification of violated DIs, i.e., their separation, may impose additional effort. Mixed: Obviously, both strategies can be mixed if a subset of DIs is added a priori and other violated DIs are separated and added dynamically. In all three cases, we do not remove non-binding DI columns even though such a clearance step would be possible. Different strategies for the addition of DIs are presented and analyzed in Section 4.2.
Separation of Violated Dual Inequalities. The separation of violated SIs is a by-product of Phase 1 of the two-phase pricing approach described above. Indeed, forward labeling in the customer network (N i , A i ) (see Figure 2) produces labels at each vertex j ∈ N i , j = 0. Recall that each such vertex j is associated with the jth commodity k in K i . One of the labels at j is L = (π i,k , d i,k ), which results from the path that only includes the jth commodity but no other commodity. Any other non-zero label at j must have included one other or several commodities S ⊆ K i so that it is identical to L = ( k∈S π ik , k∈S d ik ). If L dominates L , in the SPPRC sense of Phase 1, this implies k∈S π ik ≥ π i,k and k∈S d ik ≤ d i,k . If the first inequality is strict, the violated SI for (i, k , S) has been identified. This procedure finds a most violated SI (i, k , S) because commodities are sorted by non-decreasing demand so that any subset S ⊂ K i with k∈S d ik ≤ d i,k comprises only commodities with an index smaller than j.
Primal Solutions. If additional variables y i,k ,S for SIs are added to the RMP, the solution to the master program may consist not only of columns of routes but of a mixture of these and DI columns. Such a solution can however be transformed into a pure route-columns solution if the DIs are DDOIs. Both statements are in fact equivalent as proven in (Gschwind and Irnich, 2016, Proposition 1). An iterative transformation procedure is described in detail in (Gschwind and Irnich, 2016, Algorithm 1). The basic idea is as follows: In each iteration, the algorithm chooses from the current (partly transformed) solution a route column with coefficients a r for a route r ∈ Ω and a compatible SI column (i, k , S), i.e., one that fulfills a r ik = 1 and a r ik = 0 for all k ∈ S. Letλ r andȳ i,k ,S be the solution values of the chosen columns, respectively. Then, a new solution is constructed from the current one by decreasing the solution values of the chosen route and SI columns by min{λ r ,ȳ i,k ,S } and increasing the solution value of the column corresponding to a new route with coefficient a + e iS − e ik , where the latter are the incidence vectors of {i} × S and (i, k ), respectively. The procedure continues until no more compatible routes and SI columns exist. If all SIs that have been added are DDOIs, it is guaranteed that none of them has positive value in the last transformed solution so that a pure route-columns solution is finally found.
Over Stabilization and Recovery Procedure. Recall that, in general, SIs are neither DOIs nor DDOIs for the C-SDVRP. However, they may in fact be DDOIs for many C-SDVRP instances. Even if not, their addition does have a stabilizing effect on the column-generation process. Moreover, this possible overstabilization can be remedied typically without significant additional computational effort. Indeed, if the added inequalities are no DDOIs for the given instance, i.e., all dual-optimal solutions are cut-off, a recovery procedure described in detail by Gschwind and Irnich (2016) can be used to detect and handle overstabilization. For the C-SDVRP, the recovery procedure first tries to build from the RMP solution a pure route-columns solution. If this is not possible the RMP has been overstabilized. In this case, there exists an SI column with positive value, but no compatible route column. Let (i, k ) be the task whose commodity k is replaced by other commodities in the SI column. The recovery procedure then eliminates all SIs that replace this task (i, k ) from the RMP, forbids their re-generation in the separation routine, and restarts the column-generation process to optimize the RMP. The process iterates until a pure route-columns solution is found.

Valid Inequalities and Cutting Strategy
We use three types of valid inequalities, namely capacity cuts, subset-row inequalities, and strong-degree constraints.
For any subset of customers For any CC added to the RMP of (1), let γ C be the corresponding dual price. Then, the value γ C /2 can be distributed symmetrically on the edges (i , j ) and (i , j ) for all (i, j) ∈ δ − (C) of the undirected SPPRC pricing network. Hence, CCs are robust cuts (Fukasawa et al., 2006). The second family comprises subset-row inequalities (SR inequalities), first introduced for the VRPTW by Jepsen et al. (2008). According to the definition of rows for partitioning constraints in (1b), the variant of the SR inequalities for the C-SDVRP is defined for subsets of tasks, i.e., As in several other works, we restrict ourselves to subset R of cardinality three. The associated SR inequality for R is For elementary routes, the SR inequality for R ensures that at most one route that fulfills at least two of three tasks is part of a feasible solution. We assume that the active SR inequalities are given by the sets R ∈ R. For any R ∈ R, let β R be the corresponding dual price. To correctly incorporate the dual price β R , a binary resource has to be added to the labels of both the customer networks and the SPPRC multi-graph.
The labeling algorithms are modified as explained by Jepsen et al. (2008). The third family of strong-degree constraints (SD constraints) was introduced by Contardo et al. (2014). For the C-SDVRP, there is one SD constraint per task (i, k) with i ∈ N and k ∈ K i ensuring that (i, k) is served by at least one route, may this route be elementary or not. Formally, the SD constraint is where the coefficient ξ r ik indicates whether or not route r delivers commodity k to customer i. Hence, ξ r ik ≤ a r ik with equality for elementary routes. Strictly less occurs, when a non-elementary route r delivers (i, k) twice or even more often. For any SD constraint added to the RMP of (1), let η ik be the corresponding dual price. As explained by Contardo et al. (2014), one binary resource per active SD constraint must be added to the SPPRC multi-graph.
Impact of Cutting on Dual Inequalities. SR inequalities and SD constraints have an impact on the SI. For a specific SI defined by (i, k , S), we define for each set R ∈ R the coefficient where δ (i,k ),R = 1 if (i, k ) ∈ R and zero otherwise. This coefficient is only positive if some tasks that define the SR inequality also occur in the SI. More precisely, the inner non-negative part in (6) counts the surplus of tasks added and removed by the SI that occur in R. An even inner part exactly gives the difference in reduced cost of a route when the respective replacement is performed, while for odd values the difference is possibly overestimated. Therefore, the SI for (i, k , S) is now: where the last term on the left-hand side is the correction for the SD constraint for (i, k ).
Overall Cutting Strategy. The overall cut-generation strategy is the following: For branch-and-bound nodes up to level 10, we separate violated CCs, SR inequalities, and SD constraints in this order. To avoid the introduction of too many resources, we limit the number of SR inequalities to 30 in total and ten in each round of separation. Moreover, for each task (i, k) we limit the number of SR inequalities defined by R with (i, k) ∈ R to a maximum of two per round. For branch-and-bound nodes deeper in the tree, we only separate violated SD constraints as they are required to guarantee the completeness of the branching rule described in the following section.

Branching
We apply a similar five-level branching strategy as already applied in (Archetti et al., 2015). We briefly summarize this strategy. At the first level, we branch on the number of vehicles. At the second level, we branch on the number of visits to each customer. Here, the branching decision can make some subsets S i of commodities of customer i infeasible. In this case, we eliminate them from the customer network. For example, a single visit allows only the subset S i = K i . When forcing at least two visits, the subset S i = K i becomes infeasible. At the third level, we branch on the flow on each edge. All these branching rules are implemented by adding inequalities to the RMP. The only exception is the zero-flow decision on edges which is implemented by removing the edge.
Integer flows on all edges are however not sufficient to guarantee integer routes. Therefore, at the fourth and fifth level, we use the Ryan and Foster branching on the tasks. For two tasks (i, k) and (i , k ), let f i,k,i ,k = r∈Ω a r ik a r i k λ r be the information if (i, k) and (i , k ) are served by the same route. At the fourth level, we branch on pairs (i, k) and (i, k ) of tasks of the same customer i for which f i,k,i,k is fractional. For f i,k,i,k = 0 (separate branch), we modify the pre-computation of the Pareto-optimal deliveries S i ⊆ K i by introducing an additional binary resource into the SPPRC. Once that k or k is reached in the customer network of i the resource is set to one and the visit to the other commodity becomes impossible. The result is that all Pareto-optimal sets S i do not contain both k and k . For f i,k,i,k = 1 (together branch), the vertices k and k are merged together into one vertex of the customer network.
At the fifth level, we consider pairs (i, k) and (i , k ) for two different customers i and i . Both branches (separate and together) are implemented by introducing binary resources into the SPPRC multi-graph. The separate branch uses a single binary resource as described for the separate branch before. The together branch requires two binary resources indicating the service of (i, k) and (i , k ), respectively. The extension to the destination depot 0 and the merge of two labels are only allowed if both resources have identical value.
Completeness of Branching Scheme. The validity of the branching scheme relies on the use of the SD constraints that guarantee elementary routes. Otherwise, a minimum example with two customers 1 and 2 and tasks (1, k 1 ), (1, k 2 ), and (2, k) uses the following two non-elementary routes (0, (1, k 1 ), (2, k), (1, k 1 ), 0) and (0, (1, k 2 ), (2, k), (1, k 2 ), 0) exactly 0.5 times each. The result is that overall one vehicle is used, customer 1 is visited twice and customer 2 is visited once, and all edge flows and all f -values are integer. Therefore, none of the branching rules is applicable, but the introduction of the SD constraints for the tasks (1, k 1 ) and (1, k 2 ) cuts off this fractional solution.
Impact of Branching on Dual Inequalities. The branching decisions at the levels two, four, and five can have an impact on the validity of the DIs.
Branching decisions of level two can make some subsets S i of commodities of a customer i infeasible. In this case, we eliminate all SIs that can produce such infeasible S i . For example, in a branch that requires at least three visits to a customer i who has three commodities, any SI that replaces one commodity by two others is infeasible.
Branching decision using the Ryan and Foster branching at levels four and five impose the elimination of some SIs. In the together branch, the only valid SIs are those that either refer to none of the two tasks or SIs where the set S includes both. In the separate branch, again SIs that do not refer to either of the two tasks are valid. Moreover, SIs that replace one task by the other are also valid. All other SIs must be eliminated.
Tree Exploration and Variable Selection Strategy. We use a best-bound first tree exploration strategy in order to improve the dual bound as fast as possible. On all five levels, the specific branching variable is selected as one where the fractional part is closest to 0.5.

Acceleration Techniques
To accelerate the column-generation process, we use two types of reduced networks to solve the SPPRC pricing problem heuristically. The first one includes only one Pareto-optimal commodity set S i for each customer i in the SPPRC multi-graph. More precisely, we use the commodity combination S * i with the best reduced cost, i.e., S * i = arg min Si⊆Ki c Si i ,i . The second one uses only a subset of the inter-customer edges of the SPPRC multi-graph, namely we limit the number of edges adjacent to a customer to a given value D. In our implementation, we use D ∈ {2, 5, 10} to obtain a sequence of pricing heuristics with increasing search space and computational effort. The overall pricing strategy is to first consider the SPPRC multi-graphs with only one Pareto-optimal delivery per customer and D = 2, 5, 10, and |n| inter-customer edges, followed by the SPPRC multi-graphs with all Pareto-optimal deliveries and the same values for D. Whenever the labeling algorithm finds negative-reduced paths in one of the (reduced) networks they are returned to the RMP and the networks later in the sequence are not solved in this iteration.

Computational Results
We implemented the BPC algorithm in C++ and compiled the code in release mode under MS Visual Studio 2015 (64-bit version). At each column-generation iteration, the RMP was solved by means of CPLEX 12.6.2. The experiments were carried out on a standard PC with an Intel(R) Core(TM) i7-5930k, clocked at 3.5 GHz, and 64 GB of RAM, by allowing a single thread for each run. The time limit TL 1 for each run was set to one hour.

Benchmark Instances
As done in Archetti et al. (2015), we considered the small (n = 15) and mid-size (n = 20, 40, 60, 80) instances introduced by Archetti et al. (2016). In each instance, the customer locations are taken from the R101 or the C101 instance of the Solomon (1987) benchmark. Then, as for the small instances, the number κ of commodities is set to 2 or 3, and the probability p with which the single customer requires each specific commodity is set to 1.0 or 0.6. When the customer requires the commodity, the corresponding demand is randomly selected in the interval [1,100] or [40,60]. The vehicle capacity is finally set to β · max i∈N { k∈Ki d ik } with β ∈ {1.1, 1.5, 2.0, 2.5}. One instance is created for each combination of the defining parameters, for a total of 64 small instances.
The defining parameters for the mid-size instances differ as follows: κ and β are fixed to 3 and 1.5, respectively, and the demand for each commodity required by the customers can be randomly selected only in the interval [1, 100]. Five instances are generated for each combination of the parameters, for a total of 80 mid-size instances.
Recall that the C-SDVRP is inherently symmetric. The larger the number of commodities per customer, the larger the number of possible cost equivalent solutions to the problem that differ only in the delivery patterns associated with the same set of elementary 0-(n + 1)-path. We claimed that the use of DIs should have positive effects especially for such highly symmetric instances. In order to support this statement, we enlarged the set of small and mid-size instances considering additional values for the number κ of commodities, namely 4, 5, and 6. In particular, following the same procedure as described in Archetti et al. (2016), we generated 336 additional instances with more commodities (96 small and 240 mid-size).

Impact of Dual Inequalities
For the analysis of the impact of DIs, we restrict the test set to the 160 mid-size instances with 20 and 40 customers, because these instances pose a challenge due to the large number of up to 40 · 6 = 240 tasks. Moreover, comparisons on the basis of results obtained for smaller or bigger instances may lead to erroneous conclusions as they are generally either very easy or prohibitively hard to solve.
We define Plain as the version of the BPC algorithm in which DIs are not used. We analyze how much the DIs contribute to the performance of the BPC algorithm by comparing Plain against the following variants of the BPC algorithm in which: S-Rank: the DI columns associated with ranking inequalities (4) are statically added to the RMP; S-Rank-S-Subs: the DI columns associated with ranking inequalities (4) and a selection of the subset inequalities (5) are statically included into the RMP; D-Subs: the DI columns associated with violated subset inequalities (5) are dynamically added to the RMP; S-Rank-D-Subs: the DI columns corresponding to ranking inequalities (4) are statically added to the RMP, whereas DI columns associated with subset inequalities (5) are dynamically added when the corresponding inequalities are violated; S-Rank-S/D-Subs: as S-Rank-S-Subs but with the additional possibility to dynamically add to the master problem further DI columns associated with violated subset inequalities; We analyze the six different BPC algorithms with the help of performance profiles that allow the comparison of a set A of algorithms that are all applied to the same benchmark set (Dolan and Moré, 2002). The performance profile ρ A (τ ) of an algorithm A ∈ A is the fraction of instances that algorithm A can solve within a factor τ of the fastest algorithm (unsolved instances are taken into account with infinite run time). In particular, ρ A (1) is the percentage of instances on which A is a fastest algorithm compared to all other algorithms B ∈ A \ {A}. The value 1 − ρ A (∞) is the percentage of instances not solved by A. We start with the comparison for the solution of the linear relaxation of (1), in the following denoted by Root, for A = {Plain, S-Rank, S-Rank-S-Subs, D-Subs, S-Rank-D-Subs, S-Rank-S/D-Subs}. The performance profiles are depicted in Figure 3. The most striking result is that the Plain BPC is clearly inferior to all other variants that use DIs for a stabilization of the RMP. Without doubt, DIs significantly accelerate the column-generation process. The variant S-Rank that only uses ranking inequalities and adds them during the initialization of the RMP seems to be the clear winner. The four other variants perform better than Plain and worse than S-Rank. Among themselves, S-Rank-S-Subs seems slightly better than the three others which perform almost identically.
We expected rather similar results for integer solutions produced with the fully-fledged BPC (in the following denoted by Tree). Figure 4 depicts the performance profiles of all six BPC algorithms. Also here Plain is not at all competitive and the use of DIs can be strongly recommended. Moreover, the simple variant S-Rank is no longer superior, instead the variant S-Rank-D-Subs seems to be the overall best strategy.
In order to understand the different behavior at the root node and for the branch-and-bound tree, we deeper analyze the pricing iterations. Table 1     All five strategies reduce the total number of pricing iterations. This causes the observed reduction in computation times compared to Plain. On the other hand, all variants except for S-Rank need more exact pricing iterations. This effect is more prominent at the root node (ratios between 1.41 and 1.55) than for the complete branch-and-bound tree (ratios between 1.05 and 1.09). Only S-Rank requires less exact pricing iterations compared to Plain, at the root as well as in the tree. Therefore, the different performance shown in the profiles of Figures 3 and 4 can be explained with two concurring effects. Overall, the use of more than just the ranking inequalities (variant S-Rank) better stabilizes the RMP leading to a smaller total number of iterations, but more exact pricing iterations have to be performed. These additional exact pricing iterations occur more at the root node than in the branch-and-bound tree.
As stated above, on the basis of the performance profiles depicted in Figure 4, strategy S-Rank-S/D-Subs seems to be superior to the other strategies. This is also confirmed by the results shown in Table 1 where S-Rank-D-Subs is best w.r.t. both the number of optima and computation time (95 optima and 0.66 time ratio w.r.t. Plain). Therefore, the following experiments choose S-Rank-S/D-Subs as the stabilized version of the BPC.
Next, we show more detailed results comparing Plain and the winning variant S-Rank-S/D-Subs. Table 2 summarizes the results obtained on all 160 benchmark instances, grouped by the number κ of commodities and the number n of customers. Each group comprises 20 instances with the exception of the group for n = 40 and κ = 6, for which Plain was not able to solve the linear relaxation at the root node in one out of 20 cases. Again, Plain serves as the baseline algorithm and we provide absolute values for it, while values for S-Rank-S/D-Subs, apart from those reported in column opt, are given as geometric means of the ratios relative to Plain. The columns with header Root refer to the solution of the linear relaxation and those with header Tree to the full BPC. The columns of the table have the following meaning: #: Number of instances for which the linear relaxation is solved by both algorithms, all averages in section "Root" are taken only over those instances; iter: Number of column-generation iterations; cols: Number of columns (of routes as well as of DIs) in the RMP at termination; PP: Accumulated time (in seconds) for the solution of the pricing subproblems; RMP: Accumulated time (in seconds) for the reoptimization of the RMP; root: Total time (in seconds) for solving the linear relaxation; opt: Number of instances solved to proven optimality; time: Total time (in seconds) for solving branch-and-bound tree, averages are taken only over those instances solved optimally by both.
We start with observations regarding the solution of the root node: The larger the number n of customer and the number κ of commodities, the more iterations are needed and more columns are generated with the Plain BPC. The same is true for the pricing and re-optimization times and, therefore, for the overall root computation times. For a constant number of tasks, e.g., 120 tasks, a comparison of the rows for κ = 3 and n = 40 and for κ = 6 and n = 20 reveals that the latter instances are significantly more difficult for Plain. This is due to the higher symmetry in the latter group of instances with κ = 6 commodities. The columns for the S-Rank-S/D-Subs show that the number of iterations and generated columns as well as the computation time for pricing and re-optimization of the RMP are clearly reduced by adding DI columns. Larger instances benefit more from stabilization than smaller instances. Indeed, for the largest instances with κ = 6 and n = 40, i.e., 240 tasks, the average time ratio for the overall root computation goes down to 0.43. Comparing PP and RMP times, the RMP times decrease more. The point is probably that PP times are dominated by the calls to the exact labeling algorithm over the full network. These exact iterations cannot be avoided with stabilization, but the number of iterations and generated columns is reduced so that RMP re-optimization times are reduced by a larger extend.
The comparison of the Plain and the S-Rank-S/D-Subs regarding the branch-and-bound tree shows the following: In total, computation times of the stabilized BPC are reduced to 66 % on average over all instances. The faster computation times allow the variant S-Rank-S/D-Subs to compute four more optimal solutions (95 compared to 91 computed with Plain).

Comparison with Results of Archetti et al. (2015)
We briefly compare our chosen variant S-Rank-S/D-Subs of the BPC with the implementation of Archetti et al. (2015) which does not use any stabilization technique. All instances used in the computational study of Archetti et al. (2015) have two or three commodities (κ = 2 or 3). Since we focus on the stabilization aspect, only the instances with κ = 3 are interesting for a detailed comparison. In contrast to the previous Section 4.2, however, we use all the available instances with n = 20, 40, 60, and 80 customers.
Tables 3 and 4 present instance-by-instance results for the 80 instances, grouped by the original Solomon instance, number of customers n, and commodity probability p. The other columns of the tables have the following meaning: z * : Lower bound obtained by BPC at time of termination; z * : Upper bound obtained by BPC at time of termination; gap: The percentage gap 100 · (z * − z * )/z * if not solved to proven optimality (opt); time: Total computation time in seconds.
Note that Archetti et al. (2015) have set the computation time limit TL 2 to two hours (7200 seconds), while we set the time limit TL 1 for variant S-Rank-S/D-Subs to one hour (3600 second) to account for the better computer used in our experiments.
Over the benchmark, our stabilized BPC algorithm solves 32 out of 80 instances optimally, while the unstabilized algorithm of Archetti et al. (2015) solves 24. In particular, our algorithm provides 4 out of 20 optima for the 60-customer instances on which the unstabilized algorithm was not able to prove any optimality. However, all 80-customer remain unsolved, but we can at least provide lower bounds in all cases.

Results for New Instances with More Commodities
As mentioned at the beginning of the section, we created 96 small (n = 15) and 240 mid-size (n = 20, 40, 60, and 80) new C-SDVRP instance with four to six commodities (κ ∈ {4, 5, 6}). The new benchmark is available at http://logistik.bwl.uni-mainz.de/benchmarks.php. Table 5 shows aggregated results whereas the Appendix gives instance-by-instance results, i.e., Tables 6-8 for the small instances and Tables 9-16 for the mid-size instances. The additional columns of the tables have the following meaning: z LP : Linear-relaxation lower bound of the RMP; Instance (Archetti et al., 2015) S-Rank-S/D-Subs   Archetti et al. (2015), mid-size instances with κ = 3 (Part 1).
gap LP : The percentage gap 100 · (z * − z LP )/z * of the linear-relaxation lower bound; nodes: Number of solved nodes by BPC at time of termination; CC: Number of generated CC by BPC at time of termination; SR: Number of generated SR inequalities by BPC at time of termination; SD: Number of generated SD constraints by BPC at time of termination; rec: Number of times the recovery procedure found an overstabilized RMP solution in the BPC at time of termination. Overall, we are able to solve all but two of the 96 small instances with 15 customers and most (51 of 60) of the 20 customer instances. For the larger instances, the number of optima proven decreases significantly: only 16 and two optima for the 40 and 60 customer instances, respectively, while we cannot find any optimal solution for the instances with 80 customers. Lower bounds, however, are provided for all instances.

Conclusions
In this paper, we have developed a new BPC algorithm tailored to the commodity-constrained split delivery vehicle routing problem (C-SDVRP). The main novelty is the use of dual-optimal inequalities for the stabilization of the column-generation process. We have focused on the interaction of branching and cutting decisions and the different classes of dual inequalities. For the first time, non-robust VRP cuts such as subset-row inequalities and strong-degree constraints have been used in conjunction with dual subset inequalities that were originally introduced for packing and cutting applications. In order to keep as many dual inequalities as possible in the overall branch-and-bound tree, we have derived precise rules for their validity in the presence of branching and cutting constraints.
In an extensive computational analysis, we compared different strategies for the a-priori addition and dynamic generation of dual inequalities. A first result is that all strategies that employ stabilization based on dual inequalities clearly outperform the unstabilized version of the BPC algorithm. Interestingly, the BPC   algorithms with the analyzed strategies perform differently when solving the linear relaxation and within the branch-and-bound tree. A key observation in this context is that certain strategies require less overall column-generation iterations, but more calls to the exact pricing algorithm. This effect is, in the C-SDVRP, more pronounced at the root node resulting in the observed and different behaviors of the strategies. Our best strategy is one that a priori adds dual ranking inequalities and a few dual subset inequalities, but identifies the most violated dual subset inequalities in each column-generation iteration and adds them to the RMP. Comparing this strategy to the state-of-the-art algorithm on the existing benchmark, we are significantly faster. We have computed several new proven optima for instances with up to 60 customers and 180 tasks. In addition, we have determined lower bounds for all tested instances with up to 80 customers and 480 tasks. This improves the bounds for all unsolved instances and proves lower bounds for the first time for several instances. This appendix is supposed to become online supplementary material.

Instance
S-Rank-S/D-Subs