Exact separation of the rounded capacity inequalities for the capacitated vehicle routing problem

The family of Rounded Capacity (RC) inequalities is one of the most important sets of valid inequalities for the Capacitated Vehicle Routing Problem (CVRP). This paper considers the problem of separation of violated RC inequalities and develops an exact procedure employing mixed integer linear programming. The developed routine is demonstrated to be very efficient for small and medium-sized problem instances. For larger-scale problem instances, an iterative approach for exact separation of RC inequalities is developed, based upon a selective variable pricing strategy. The approach combines column and row generation and allows us to introduce variables only when they are needed, which is essential when dealing with large-scale problem instances. A computational study demonstrates scalability of the proposed separation routines and provides exact RC-based lower bounds to some of the publicly available unsolved CVRP instances. The same computational study provides RC-based lower bounds for very large-scale CVRP instances with more than 3000 locations obtained within appropriate computational time limits.


INTRODUCTION
The Vehicle Routing Problem is a well-known problem in the operations research domain, which goes back to [4].The problem seeks to find a set of routes for a number of vehicles (originally located at a depot location denoted by 0) to serve a set of customers denoted by .A route is defined as a cycle in the complete undirected graph G = ( ∪ {0}, E) passing through the depot, where E denotes the set of edges of the graph: E = {{i, j} | i, j ∈  ∪ {0}, i < j}.It is assumed that every customer has some deterministic demand for a commodity  i > 0, i ∈ , and every vehicle has identical capacity to serve that commodity Q > 0. Hence, the routes must satisfy the vehicle capacity constraint.In addition, this paper considers the simplest version of the problem, where delivery cannot be split between two (or more) vehicles, which subsequently implies that Q ≥ max i∈  i .The objective function of the problem is defined using the distance (or cost) matrix C = {c ij ≥ 0 | {i, j} ∈ E}, such that any feasible routing solution is evaluated in terms of its total distance using matrix C. Thus, the problem aims to find a solution with the minimum total travel cost.Finally, the focus of this paper is on the version of the problem with a symmetric matrix C, which is known in the literature as the symmetric Capacitated Vehicle Routing Problem.Typically and with very little loss of generality (at least from the practical standpoint), nonnegative parameters c ij ,  i and Q are further assumed to take integer values.This is the assumption used in this paper as well.
One of the early mathematical programming formulations of the symmetric CVRP with an unconstrained number of vehicles was proposed by [11]: x(S) x 0i ∈ {0, 1, 2}, {0, i} ∈ E, (1.4) where x ij is a decision variable denoting whether a vehicle traverses edge {i, j} in any direction.Such a decision variable is binary for {i, j} ∈ E, i ∈ , and may take, in addition to possible values 0 and 1, the value of 2 for an edge connecting the depot location to a customer: this case corresponds to the simplest one-customer route 0 − i − 0. The set of constraints (1.2) ensures that every customer is visited exactly once.The set of constraints (1.3) ensures that no integer solution contains subtours not connected to the depot, and that the vehicle capacity constraint is satisfied.Notation x(S) represents the summation of all edge variables x ij within subset S, x(S) = ∑ {i, j}∈E(S) x ij , (1.6) where E(S) = {{i, j} | i, j ∈ S, i < j} for S ⊂ .Hence, constraints (1.3) generalize the well-known subtour elimination constraints of [5] for the Traveling Salesman Problem to the CVRP.Constraints (1.3) are often used in its equivalent but different form in terms of constraints on the capacities of cut-sets: where S =  ∪ {0} ⧵ S, (S, S) is a cut-set associated with the set of customers S ⊂ , and x(S, S) = j∈S, j<i x ji .The equivalence of (1.3) and (1.7) is easy to see.Consider the aggregation of constraints (1.2) over i ∈ S: which can be represented as follows: From this representation, it is easy to see that a single constraint in the set (1.3) is equivalent to the corresponding constraint in (1.7) when each of the constraints is considered together with the corresponding subset of degree constraints (1.2).Note, however, that the family of RC inequalities is in fact the relaxation of the stronger set of constraints called weak capacity inequalities, where r(S) is the optimal objective of the corresponding bin-packing problem [8] for set S. Since the bin-packing problem itself has the same NP-hard complexity as the CVRP, the value of r(S) is often relaxed to ⌈ ∑ i∈S  i ∕Q ⌉ , the round-up of the fraction ∑ i∈S  i ∕Q.This is why the corresponding constraint is called the Rounded Capacity (RC) inequality, a simple and easy-to-compute lower bound for the number of required vehicles to serve the subset of customers S. Such inequalities constitute the main focus of this paper.A more detailed discussion about constraints of this type can be found in [14].
Given the exponential size of the constraint pool (1.7), such constraints are imposed in a cutting planes optimization framework.Using a cutting planes framework requires a separation routine to find RC inequalities violated by a given solution, which can be either heuristic or exact.An efficient heuristic algorithm for separation of RC inequalities was proposed in [15].One of the most efficient publicly available codes CVRPSEP [13] used heuristic approaches for separation of the RC inequalities, together with other families of valid inequalities.The question of using an exact versus a heuristic separation routine is typically influenced by the computational complexity of the separation problem.While NP-hardness of the RC inequalities separation problem had been assumed since 1990s [1,15], the first formal study proving this claim appeared only recently, [6], which justifies the wide use of heuristic approaches for separation of RC inequalities.At the same time, given the significant progress in development of modern mixed integer linear programming solvers, exact separation of RC inequalities using mixed integer programming is a reasonable idea that has received limited attention in the literature.This paper argues that mixed integer linear programming is indeed a promising approach for the solution of the RC cuts separation problem.
The outline of the paper is as follows.Section 2 of this paper develops a new mixed integer linear programming formulation for exact separation of RC inequalities.Section 3 develops an iterative approach for exact separation of RC inequalities, based upon a selective variable pricing strategy, combining simultaneous cut and column generation and applicable for large-scale instances.Section 4 presents additional implementation details of the iterative approach.Section 5 evaluates the performance of the proposed approaches and demonstrates the possibility to find all violated RC inequalities for problem instances with up to 500 customers within reasonable CPU times.The iterative separation approach appears to be scalable for problem instances with more than 1000 locations.For example, Section 5 provides results of computational experiments, where lower bounds for instances with up to 7000 locations within 3 or 5 h of computer runtime are obtained.Section 6 concludes.

EXACT SEPARATION OF THE ROUNDED CAPACITY INEQUALITIES
Suppose there is a set of potentially fractional decision variables {x * ij , {i, j} ∈ E}, and we want to determine exactly whether a violated RC inequality, that is, a cut, exists or not.A violated RC inequality corresponds to a cut-set (S, S).Let  i ∈ {0, 1} indicate whether i ∈ S or not.Then, an edge {i, j} belongs to the cut-set (S, S) if and only if i and j are in different sets S vs S.This condition can be modelled as  i +  j = 1.Alternatively, condition  i +  j = 1 is equivalent to another more convenient condition  i +  j − 2 i  j = 1, which allows us to express the capacity of the cut-set (S, S) as follows: subject to ) ) ) where  ij =  i  j , which is linearized using one of the standard techniques due to [9].It follows that the mixed integer linear program for the identification of a RC cut with the maximal violation can be represented as follows: subject to ) ) ) Constraint (2.7) is the key to the separation routine: the optimal  +1 must equal ⌈ ∑ i∈S  i ∕Q ⌉ due to the objective function (2.6).Presence of an  1 ∈ (0, 1] in the left-hand side of (2.7) is crucial for the following reason.Let S be such that ∑ i∈S  i = kQ for an integer k.Then, constraint (2.7) would imply  = k, which would result in the incorrect target capacity in the objective (2.6) equal to 2(k + 1) instead of 2k.At the same time, when (k − 1)Q ≤ ∑ i∈S  i < kQ, then (2.7) implies  = k − 1 and corresponds to the correct target capacity of 2k for such S in the objective (2.6).Hence, an  1 ∈ (0, 1] "adjusts" the target  from k to k − 1 for cases when ∑ i∈S  i = kQ and makes RCC−Sep an exact separation routine.In addition, RCC−Sep omits previously introduced (2.4) constraints as being fulfilled at optimality for such an objective function that aims to maximize values of  ij variables.Note that a practical implementation of RCC−Sep should introduce a tolerance parameter  2 > 0 and consider edges {i, j} ∈ E() ∶ x ij ≥  2 only.Specific values of parameters used in our computational study are presented in Section 4.
The mixed integer linear programming-based approach is not totally new to the literature of the RC cuts separation problem.Fukasawa et al. [7] proposed a slightly different approach for finding a RC cut with the maximum violation based on mixed integer linear programming.However, the target capacity  in their approach was declared as an integer-valued parameter, which required executing the procedure multiple times with different values of the parameter.Instead, the RCC−Sep procedure is complete in this sense and searches for violated inequalities for every possible , which speeds up the computational run-time significantly.Moreover, computational performance of the new RCC−Sep procedure can be further enhanced using the following simple observation.When the ultimate aim is to find all violated RC inequalities, RCC−Sep needs to be executed multiple times, adding one violated RC inequality per iteration.The number of iterations can be decreased if we could add more than one RC cut per iteration.The following callback is introduced to the RCC−Sep procedure: add (S, S) to the set of violated RC constraints. ( A single iteration of solving RCC−Sep with Callback enabled results in the generation of multiple RC cuts before another iteration of RCC−Sep is launched and, in general, reduces the overall number of required iterations and hence the CPU time.
It is important to note that Callback is implementable using the standard tools modern optimization solvers provide.

EXACT SEPARATION OF THE ROUNDED CAPACITY INEQUALITIES FOR LARGE-SCALE PROBLEM INSTANCES
RCC−Sep is a mixed integer linear programming program that grows quadratically in size with the number of customers.Hence, solving the separation program to optimality might be a great challenge for larger problem instances.For this reason, a special technique might be desirable to develop in order to handle such problem instances.We adapt a technique described by [3] for finding an optimal solution to the linear programming relaxation of the [5] formulation of the symmetric Traveling Salesman Problem involving a small subset of edges only.The approach starts from solving the linear programming relaxation of the (1.1)-(1.5)model defined over a subset E ′ of x ij variables (referred to as the master problem), which makes solving the corresponding RC cuts separation problem over the reduced set of edges easier.Then, the initial set of edges is evaluated to see whether additional edges (variables) need to be added to the master problem, and the initial set of edges is expanded if the evaluation is positive.The procedure continues until it is proved that x ij variable must equal 0 for every other edge in any optimal solution of the linear programming relaxation of (1.1)-(1.5).Essentially, we have developed a procedure based on the column generation principles for solving linear programming problems and refer the reader to the remark in the end of this section for an alternative description of the nature of this procedure.
The linear programming relaxation of (1.1)-(1.5)we aim to solve is presented below together with associated dual variables: ) ) For a subset of edges E ′ ⊂ E with E ′ (S) = E(S) ∩ E ′ , consider the following restricted version of P over the subset of edges E ′ : ∑ {i, j}∈E ′ (S) ) ) As we cannot solve P ′ directly due to the size of constraint pool (3.9), we can use RCC−Sep to identify a (often) small subset of essential constraints  ⊂ 2  , sufficient to solve P ′ .This is done by dropping constraints (3.9) from P ′ and solving the resulting linear program first, then iteratively running RCC−Sep and incorporating identified RC cuts into P ′ .Linear program P ′ is solved when no more violated RC inequalities exist.Set  consists of all subsets S corresponding to identified RC cuts, and P ′′ is the resulting tractable linear program that can be solved directly: ∑ {i, j}∈E ′ (S) ) ) An optimal solution x ′′ to P ′′ is also optimal to P ′ and it is guaranteed that no constraint (3.9) for S ∉  is violated in P ′ .At the same time, the obtained optimal solution to P ′ for a given E ′ does not necessarily allow us to construct a corresponding optimal solution to P. Therefore, the aim of the procedure is to construct a sufficiently small set E ′ , such that solving P ′ is tractable and at the same time to ensure that an optimal solution to P ′ extends to an optimal solution to P.
To facilitate execution of the procedure, it is important to choose E ′ such that the restriction problem P ′ is feasible; we discuss how to establish that in the next section.Assuming this is the case, we extend an optimal solution x ′′ to P ′ to a feasible solution x * to P: The aim for x * is to be not only a feasible but also an optimal solution to P. The following paragraphs provide explicit theoretical arguments on how to determine whether a given set E ′ is sufficiently large, that is, that the restricted problem P ′ solves the original problem P. One can use the linear programming duality theory to check if this is the case.The linear programming duals of P, P ′ , P ′′ are available and can be expressed as follows: ) .27) ) ) Y S ≤ 0, S ∈  , (3.36) where  S (i) denotes the standard indicator function  S (i) = 1 if i ∈ S and 0, otherwise.Let ( y ′′ , Y ′′ , u ′′ ) denote an optimal solution to D ′′ that corresponds to an optimal solution x ′′ to P ′′ .This solution extends naturally to an optimal solution to D ′ : is indeed optimal to D ′ for the following two reasons: first, it is clearly feasible to D ′ and second, its objective value is equal to the optimal objective value of the dual problem P ′ .Subsequently, ( y ′ , Y ′ , u ′ ) extends to a potentially feasible solution (y * , Y * , u * ) to the problem D as follows: At the same time, by the weak duality of linear programming for the pair P ←→ D we have: Subsequently, taking into account the way solutions x * and (y * , Y * , u * ) were constructed, we obtain which, together with (3.44), implies that The last identity (3.47) demonstrates the primal-dual pair of feasible solutions with identical objective functions, which proves that x * is an optimal solution to P. 2. If (y * , Y * , u * ) is infeasible to D, some of the constraints (3.21), (3.22) must be violated by this solution.Violations are only possible for constraints corresponding to edges {i, j} ∈ E⧵E ′ .Thus, we need to find and subsequently add all edges {i, j} corresponding to violated constraints (3.21), (3.22) to E ′ .Notice that by construction of (y * , Y * , u * ), the following simplified inequality has to be checked: for every edge {i, j} ∈ E ⧵ E ′ , i ∈  in order to find violated constraints in (3.22), and inequality for every edge {0, i} ∈ E ⧵ E ′ in order to find violated constraints in (3.21) (since u * ij = 0 for every {i, j} ∈ E ⧵ E ′ ).
The complete procedure to solve P using an initial subset E ′ of edges via column generation can now be summarized as follows: RCC−Sep Col. Gen. ∶

3.
Solve the linear programming problem P ′ using the exact RC cuts separation approach described in the previous section and obtain an equivalent problem P ′′ with a reduced set  of RC inequalities.Let x ′′ be its (arbitrary) optimal solution and find the corresponding dual optimal solution (y ′′ , Y ′′ , u ′′ ).

Remark. The above presented procedure has an interpretation coming from the revised simplex method for linear programming. Notice that expression c ij − y
is in fact the reduced cost for variable x ij currently not in the model P ′ .Negativity of the reduced cost in the simplex method for solving linear programming programs implies that the objective function can be improved when that variable is added to the basis.This is exactly the same criterion for adding an edge to the E ′ set.Hence, it can be argued that RCC−Sep Col. Gen. is simply a modification of the revised simplex method, where all variables with negative reduced costs are added to the model instead of just one as it is in the revised simplex method.Moreover, no variable ever leaves the model, which is why the number of iterations of the procedure is polynomially bounded from the above.

RCC−SEP VIA COLUMN GENERATION: IMPLEMENTATION DETAILS
While the previous section contains a description (together with an underlying theoretical justification) of a fully functional procedure for finding RC-based bounds for the CVRP, a practical implementation of the procedure requires additional details.Our implementation of the procedure is available at the following link https://bitbucket.org/jonx0878/rcc_separation/src/master/.
This section provides additional details about our implementation.The implemented procedure differs slightly from the one developed in Section 3. The main difference is due to the existence of equivalent representations of RC inequalities.In addition to already discussed equivalent versions of RC inequalities (1.3) and (1.7), the CVRPSEP manual [12] provides yet another version of a RC inequality: where S =  ⧵ S. In practice, a given RC cut should be incorporated in the form that is most appropriate.The notion of appropriateness in the context of linear programming deals with the sparsity of the associated constraints.
∑ {i, j}∈E ′ (S) ) ) ) where  1 (E ′ ) contains sets S whenever |E ′ (S)| ≤ |E ′ (S)| and  2 (E ′ ) contains the remaining ones.The corresponding linear programming dual then takes the following form: subject to ) ) ) Thus, the following conditions represent the actual criteria used for inclusion of an edge into the set E ′ in our implementation: The initial set E ′ must be chosen in order to implement a working RCC−Sep Col. Gen. procedure.As mentioned above, the set E ′ should be such that corresponding P ′ is feasible.This requirement can be satisfied by including in E ′ all the edges that appear in the current best known solution (BKS) to the problem instance.In addition, we find it effective from the computational point of view to have additional edges included into initial set E ′ .Our choice of additional edges was to identify m shortest edges for every node and include them to E ′ .The value of parameter m used in our computational experiments was 15.
Finally, we would like to emphasize that the dual simplex method was used explicitly to solve (and resolve upon adding cuts) the corresponding linear programs in all implemented procedures.This choice of the solution method is standard in the setting when linear programs are changed and resolved iteratively as it preserves feasibility of the optimal basic feasible solution obtained in the previous iteration.Finally, we employed the following values of the parameters:  1 = 0.5 and  2 = 10 −4 while all other inequality conditions (such as (2.15), (3.50), (3.52), (4.15)) have been checked with the precision  3 = 10 −3 .

COMPUTATIONAL EXPERIMENTS
This section presents results of computational experiments with the presented RC cuts separation routines.The corresponding linear and mixed integer linear programming programs for RCC−Sep were implemented and solved using Anaconda 3 64-bit Python 3.8.5 in connection with Gurobi 10.0.0 [10].Computational experiments with RCC−Sep were performed using an AMD EPYC 7601 2.20 GHz, 16GB RAM with four virtual cores (out of 32) machine equipped with Windows Server 2019.Separate code was developed to run computational experiments with CVRPSEP; it is available at https://github.com/Jonx0878/CVRPSEP-for-RCC.Computational experiments with CVRPSEP were performed on a PC with AMD Ryzen 7 5700U processor and 16GB RAM operating under Windows 11.
To begin with, a standard set of small CVRP instances considered by [2] was tested.Table 1 provides results of computational experiments with the RCC−Sep and RCC−Sep+Callback procedures using the aforementioned set of instances.Table 1 suggests that exact RC-based bounds can be easily obtained for problem instances with up to 100 locations.Moreover, Callback is a valuable addition to the separation routine that decreases the CPU times, roughly speaking, by the factor of 2. In addition, the table provides lower bounds obtained using the CVRPSEP package for separation of the RC inequalities for comparison, demonstrating that RCC−Sep with Callback provides the same or slightly better results compared to CVRPSEP while taking less computational time.
In order to test scalability of the newly proposed RC cuts separation routines further, larger problem instances have been considered.Specifically, we considered some of the problem instances introduced in [16].Experiments indicate that such instances can not be handled using the CVRPSEP package.
Table 2 provides results of computational experiments for that set of larger CVRP instances using RCC−Sep, with and without Callback.We observe that at the scale of over 300 locations, exact RC cuts separation might become intractable, especially if the number of required vehicles is large.Yet, Callback improves the CPU times by the same approximate factor of 2. Some of the problem instances considered in Table 2 are not currently solved to optimality, therefore the lower bounds presented are expressed in % of the objective values of the currently best known solutions (BKS).The current BKS for those instances are presented in Table 3.
Table 2 reveals that significant computational effort might be required to find an exact lower bound using the RCC−Sep procedure.However, we would like to note that this excessive potential CPU time should not be a limitation to applicability of the proposed separation procedures.Both the RCC−Sep and RCC−Sep+Callback procedures add new RC cuts over the time of execution, and that is why at any given time provide lower bounds to the corresponding CVRP instance, which improve over time upon addition of new cuts.Figure 1 provides a specific example of how the RCC−Sep lower bound evolves over time for a given problem instance.The figure suggests that a reasonably good bound can be obtained prior to the formal termination of the procedure when no other violated RC inequality can be obtained.
Table 3, subsequently, provides results of experiments with RCC−Sep+Callback solved via the column generation approach described in this paper.It can be observed that RCC−Sep+Callback Col. Gen. improves the CPU times for all the problem instances, sometimes substantially.
In order to test scalability of the RCC−Sep+Callback Col. Gen. approach, we have attempted to solve P for larger problem instances.Specifically, we considered the set of DIMACS (2021) problem instances.These problem instances are publicly available in the repository http://vrp.atd-lab.inf.puc-rio.br/index.php/en/.Table 4 provides results of computational experiments running RCC−Sep+Callback via column generation.The table demonstrates that it is possible to find exact RC-based bounds for instances with up to 400 and even 500 locations, especially if the number of required vehicles is low.However, finding exact RC-based lower bounds for larger problem instances from this set appears to be computationally challenging.For this reason, we have introduced the time limit parameter to our cuts separation procedure.Specifically, incorporation of the time limit works as follows: the cut separation procedure over a reduced set of edges E ′ operates until either the exact RC-based 43) is justified as follows.Notice that x * ij = 0 for any {i, j} ∈ E⧵E ′ by definition(3.19).In this case, clearly, the constraint (3.5) (or constraint(3.4)when i = 0) is nonbinding.Therefore, by the complementary slackness property of linear programming, the associated dual variable u * ij must be equal to 0. Two possibilities for (y * , Y * , u * ) exist: (1) (y * , Y * , u * ) is also feasible to D and (2) (y * , Y * , u * ) is infeasible to D. These outcomes lead to the different observations 1 and 2.

1 .
Assume that (y * , Y * , u * ) is a feasible solution to problem D. Then by the strong duality of linear programming for the pair P ′ ←→ D ′ we have: The form (1.3) of a RC inequality is sparser when |S| ≤ n(n + 1)∕(2n − 2), while (4.1) is sparser otherwise.Therefore, our implementation of the basic RCC−Sep procedure adds (1.3) cut when |S| ≤ n(n+1)∕(2n−2) and uses form (4.1), otherwise.In case of the RCC−Sep Col. Gen. procedure, we used a simpler logic to determine which form of a RC cut to use.Specifically, we compared |E ′ (S)| versus |E ′ (S)|, and (1.3) form was used whenever |E ′ (S)| ≤ |E ′ (S)|, and (4.1) was used otherwise.The adjusted version of P ′′ is as follows:

TABLE 1
Results of computational experiments with RCC−Sep, RCC−Sep+Callback procedures and CVRPSEP applied to a set of small CVRP instances.
Note: Lower bounds are expressed in % of the optimal objective values for corresponding instances.