The Two-index Subtour Formulation

Run LPL Code  ,  PDF Document

Problem

The capacitated vehicle routing problem (CVRP) is defined as follows. Given K > 0 vehicles with capacity C > 0, and a network (graph) G = (V,E,d) with the node set i,j ∈ V = {1…n}, an edge set E and a (time)-distance matrix d_i,j, let 1 ∈ V be the depot (starting node). Furthermore let c_i the amount of capacity to load (or unload) at location i. The K vehicles must start at the depot and return to the depot, making a tour by delivering (or collecting) goods to/from costumers located at the nodes on the graph G. Each customer must be visited by exactly one vehicle. Find the shortest round-trip by adding the round trips of all vehicles. Implement a mathematical model for this problem. For more information see [5], [1], or [3]).

Modeling Steps

At first view, this problem is similar to the multiple TSP (see tsp-m), that is, a TSP with multiple salesmen starting all from the same location (the depot). Unfortunately, there is an important difference: each vehicle has a limited capacity and, therefore, can only visit a reduced number of customers.

Interestingly, the problem can be formulated similarly to the cut set formulation of the TSP (here we define the problem on an undirected graph).

1. Let’s introduce a binary variable x_i,j which is 1 if an edge is traversed by a vehicle, otherwise it is 0.

2. We also define a value for all subsets S ⊆ V :

   

   b(S) expresses the minimal number of vehicles that are needed to serve all the customer in S from the capacity point of view.

3. The first constraint A (below) says that the number of vehicles leaving and arriving at the depot is 2K.

4. The second constraint B says that all customers (except the depot) must be visited exactly once (once arriving, once leaving)).

5. Constraint H excludes all subtours and all capacity exceeding subsets. Suppose for a moment that b(S) = 1 for all subsets S ⊂ V , then the constraint say that no subtour is allowed as in the TSP, that is we require that at least two edge leave (or enter) any subset S, it corresponds to the subtour elimination constraints of the TSP. Suppose now that for a particular S we have b(S) = 2, that is, at least two vehicles are needed to serve the customers in S. This means that at least 4 edges must leave or enter the subset S (because 2 vehicles are needed to serve all costumer in S. Hence, constraint H would be violated otherwise, because its capacity is exceeded. That is, constraint H prevents subtours and tours which exceeds the capacity of C (even if they start and end at the depot).

The model is as follows:

Further Comments

The model has an exponential number of constraints as already seen in the TSP formulation. We may try to solve the problem in a similar way: Starting without any subset constraints H, and then add repeatedly subsets of violating constraints. The model formulation is only slightly more complicated than the TSP: The constraint H not only prevents subtours that do not include the depot, but also prevents over-capacitated tours, as described. This is done by first finding the components of the graph as in the TSP approach and collect them, except the component that contains the depot. In a second step within the loop, we remove all edges with the depot as an endpoint and find all components with nodes S for which b(S) > 1. All subset are added to the subsets list of violating constraints and the problem is solved, and the process is repeat until no such component is found anymore.

model cvrp1 "The Two-index Subtour Formulation"; 
  set i,j,k   "nodes set"; 
  parameter c{i}; D{i,j}; K; C; 
  binary variable x{i,j|i<>j}; 
  constraint 
    A{i}: sum{j} x[i,j] =  if(i=1,K,1); 
    B{i}: sum{j} x[j,i] =  if(i=1,K,1); 
    S2{i,j|i<j}: x[i,j]+x[j,i] <= 1; 
    H{s|s<=cNr and bS}: sum{i,j|S[s,i] and ~S[s,j]} (x[i,j]+x[j,i]) >= 2*bS; 
  set s:=[1..2000]; 
  parameter cNr; cNr1:=2; cNr2; passed; CO{i}; size{i}; 
    Count{i}; maxS; OBJ; 
    S{s,i}; bS{s}; 
  while (cNr1>1 or cNr2) and passed<=20 do 
    minimize obj: sum{i,j} D*x; OBJ:=obj; 
    passed:=passed+1; 
    cNr1:=Graph.Components(x,CO); 
    {i|CO<>CO[1]} (S[CO+cNr,i]:=1, bS[CO+cNr]:=1); 
    cNr:=cNr+cNr1; 
    {i} (x[1,i]:=0, x[i,1]:=0); 
    cNr2:=Graph.Components(x,CO); 
    size{i}:=0; cNr2:=0; 
    {i} (size[CO]:=size[CO]+c); 
    {i|size>C} (cNr:=cNr+1, cNr2:=cNr2+1, 
                {j|CO[j]=i} (S[cNr,j]:=1), bS[cNr]:=Ceil(size/C) ); 
  end; 
  minimize obj1: sum{i,j} D*x; //optimize a last time 
end

Solution

This approach works well for the schoolbus problem. However, even for small problems from the VRP-library found here , it is difficult to find the optimal solution. The code lists three problems instances from the library: Problem 3 (A-n32-k5.vrp) has 32 nodes, 5 vehicles with capacity of 100. Gurobi 11.0 does not find the optimum (784) even after solving 121 mip models and adding 165 subtour constraints (H). A more appropriate approach for solving larger instances is the model cvrp1

The data model is defined as:

  model data; 
    string File:=['P-n19-k2.vrp'];  //opt=212 (problem 1) (K=2) 
    --string File:=['E-n22-k4.vrp'];//opt=375 (problem 2) (K=4) 
    --string File:=['A-n32-k5.vrp'];//opt=784 (problem 3) (K=5) 
    parameter X{i}; Y{i}; dum;; 
    Read(File&',%1;-1:CAPACITY',dum,dum,C); 
    Read{i}('%1:NODE_COORD_SECTION:DEMAND_SECTION',i,X,Y); 
    Read{i}('%1:DEMAND_SECTION:DEPOT_SECTION',i,c); 
    K:=Ceil(sum{i} c/C); 
    D{i,j} := Trunc(Sqrt((X[i]-X[j])^2+(Y[i]-Y[j])^2)+0.5); 
    parameter min0:=min{i,j|D} D; max0:=max{i,j|D} D; 
    set e{i,j}:=D<2*min0; 
    set e1{i,j}:=D>=max0/1.5; 
    set e2{i,j}:=D<=3.5*min0; 
  end

Larger problem cannot be solved in reasonable time using this approach. However, the method suggests to use this model as an lower bound generator, and together with a “repair” heuristic.

Questions

1. Replace the data and the output model with the data of the schoolbus problem and run. How many times the loop is passed and how many subtour elimination constraints are added?

Answers

1. The run loops 6 times and 6 subtour elimination constraints are added until it finds the optimal solution.

References