Trouble with optional nodes in a VRP model using Gurobi for python

Question

I am currently working on a VRP problem where the vehicles are electric. This means that I have range restrictions, and charging stations as nodes which are optional to visit.

In my model I have the following:

O = [0] # depot
N = [1,2,3,4,5,6,7] # deliveries
V = [0,1,2,3,4,5,6,7] # all nodes which must be visited
F = [8,9] # charging stations.
G = [0,1,2,3,4,5,6,7,8,9] # all the nodes

In addition I have a distance matrix with distances between all nodes. I also have demand in weight at every delivery point.

My current model looks like this:

m = Model('E-VRP')
x = m.addVars(G, G, vtype=GRB.BINARY)
lw = m.addVars(G, vtype=GRB.CONTINUOUS) # for calculating load weight
d = m.addVars(G, vtype=GRB.CONTINUOUS) # for calculating distance driven

m.setObjective(quicksum(distance[i,j]*x[i,j] for i in G for j in G if i!=j), GRB.MINIMIZE)

# Restrictions

m.addConstrs(quicksum(x[i,j] for j in V if j!=i)== 1 for i in N)  #1
m.addConstrs(quicksum(x[i,j] for i in V if i!=j)== 1 for j in N)  #2

m.addConstrs ((x[i,j]==1) >> (lw[i] + Dw[i] == lw[j])  #3
              for i in G for j in G if i!=0 if j!=0)
m.addConstrs (lw[i] >= Dw[i] for i in G);  #4 (Dw being demand in weight at delivery)
m.addConstrs (lw[i] <= Qw for i in G);     #5 (Qw being weight capacity of vehicle)

m.addConstrs ((x[i,j]==1) >> (d[i] + distance[i,j] == d[j])
              for i in G for j in G if j>0);   #6
m.addConstrs (d[i] <= Qr for i in G);          #7 (Qr being range restriction of vehicle)

A combination of the first three restrictions make the model a working capacitated vehicle routing problem. It starts and ends at the depot (0), and if Qw is less than lw, the model will make multiple routes.

The problem comes when I impose range restrictions, where I have two problems:

First, the model will not make charging stations part of the route. I understand that restrictions #1 and #2 state that there must be a connection from list V to N, and vice versa, meaning that for node j, there must at be an i in list V (which does not include charging stations). However, whenever I try to "loose up" the two first restrictions, the model will suggest subtours that do not start nor end at the depot. Does anyone have any suggestions as to how I can solve this?

Second, my variable d, which keeps track of the distance travelled, will only add up the distance to the last node prior to returning to the depot. I realize that my code j>0 states that the depot cannot be a destination in my travel variable, however, this is the only way I have had it to work for all the other nodes.

Answer 1

Welcome to OR.SE.

It is not easy to enter into your problem and your model. You have modeled your problem in MILP, with each edge represented as a 0-1 variable. There are different possible formulations. All will be tedious to implement, even more if you have many constraints to implement. To make it easier and more scalable, we suggest you follow a heuristic column generation approach: generate some feasible "good" routes and then set a MILP to select the best subset of routes.

As explained above, Boolean modeling approaches are tedious to write and don't scale well for this kind of problem. Your problem can be modeled compactly by following a list-based modeling approach instead of the classical Boolean modeling approach, as you described in your question. This is a modeling approach offered by LocalSolver, which is different from traditional MILP solvers. Note that LocalSolver is commercial software. Nevertheless, it is free for faculty and students.

Below is the code snippet to model the Capacitated Vehicle Routing Problem (CVRP) by using the LocalSolver modeling language, namely LSP:

function model() {
  // Sequence of customers visited by each truck
  customersSequences[k in 1..nbTrucks] <- list(nbCustomers);

  // All customers must be visited by the trucks
  constraint partition[k in 1..nbTrucks](customersSequences[k]);

  for [k in 1..nbTrucks] {
      local sequence <- customersSequences[k];
      local c <- count(sequence);

      // A truck is used if it visits at least one customer
      trucksUsed[k] <- c > 0;

      // The quantity needed in each route must not exceed the truck capacity
      routeQuantity <- sum(0..c-1, i => demands[sequence[i]]);
      constraint routeQuantity <= truckCapacity;

      // Distance traveled by truck k
      routeDistances[k] <- sum(1..c-1, i => distanceMatrix[sequence[i - 1]][sequence[i]])
         + (c > 0 ? (distanceWarehouse[sequence[0]] + distanceWarehouse[sequence[c - 1]]) : 0);
  }

  nbTrucksUsed <- sum[k in 1..nbTrucks](trucksUsed[k]);

  // Total distance traveled
  totalDistance <- sum[k in 1..nbTrucks](routeDistances[k]);

  // Objective: minimize the number of trucks used, then minimize the distance traveled
  minimize nbTrucksUsed;
  minimize totalDistance;
}

You can write the same model in Python, Java, C#, or C++. Just have a look at the example here.

Having followed this modeling approach, LocalSolver finds quality solutions in minutes on a standard computer for instances with thousands of clients to serve. LocalSolver embeds and combines neighborhood search, constraint programming, and mixed-integer linear programming techniques under the hood to get such results.

Answer 2

In addition to @LocalSolver's answer, I believe you can solve your problem relatively easily with the or-tools routing library (free and open source). At its core, this library solves a TSP, over which you can add constraints with a resource based logic.

For example, the load on a vehicle is a resource, this resource is incremented when visiting a node, and the accumulated value of this resource cannot exceed a given value (the vehicle's capacity).

In your case, you could use a resource to keep track of the consumption of the electric charge of the battery, which would be incremented when using an edge with endpoint any node different from a charging station. Conversely, when passing through a charging station, the accumulated electric charge would go back to $0$.

Also, since visiting a charging station is optional, you could use the drop node option to allow skipping any charging station.

Or-tools is based on constraint programming, heuristics, and local search, so as stated by @LocalSolver, you can expect much better results than with a MILP when the problem size grows.