The Electric Vehicle Routing Problem (EVRP) Model is always infeasible
Answered
I am trying to implement the Electric Vehicle Routing Problem (EVRP) model described in the paper by David Woller, Victor Kozak and Miroslav Kulich (2020). Here is an image describing the constraints: 
Now, I implemented such model in Gurobipy but I keep getting an error saying it is infeasible (I created a super simple benchmark instance to use with the model):
from gurobipy import Model, GRB, quicksum
import numpy as np
import matplotlib.pyplot as plt
D = [0] # D = Depot
I = [1, 4] # I = Customers
V = [0, 1, 2, 3, 4] # V = Depot U Customers U Charging Stations
F = [2, 3] # F = Charging Stations
FuD = [0, 2, 3] # FuD = Depot U Charging Stations
demand = [0, 5, 0, 0, 3] # Demands (only customers have demand > 0)
b = {i: demand[i] for i in V} # Dictionary of Demands
A = [(i, j) for i in V for j in V if i != j] # Set of arcs
C = 10 # Vehicle capacity
Q = 15 # Battery capacity
w = {
(0, 1): 20, # Weights for each arc.
(0, 2): 12,
(0, 3): 13,
(0, 4): 30,
(1, 0): 20,
(1, 2): 2,
(1, 3): 1,
(1, 4): 3,
(2, 0): 12,
(2, 1): 2,
(2, 3): 5,
(2, 4): 11,
(3, 0): 13,
(3, 1): 1,
(3, 2): 5,
(3, 4): 3,
(4, 0): 30,
(4, 1): 3,
(4, 2): 11,
(4, 3): 3,
}
mdl = Model("EVRP")
x = mdl.addVars(A, vtype=GRB.BINARY, name="currentArc") # CURRENT ARC???
u = mdl.addVars(V, vtype=GRB.CONTINUOUS, name="currentCap") # CURRENT CAPACITY
y = mdl.addVars(V, vtype=GRB.CONTINUOUS, name="batteryLevel") # CURRENT BATTERY
# Eqtn 1 --> Objective Function
mdl.modelSense = GRB.MINIMIZE
mdl.setObjective(quicksum(w[i, j] * x[i, j] for i, j in A))
# Eqtn 2
mdl.addConstrs(quicksum(x[i, j] for j in V if j != i) == 1 for i in I)
# Eqtn 3
mdl.addConstrs(quicksum(x[i, j] for j in V if j != i) <= 1 for i in F)
# Eqtn 4
mdl.addConstrs(
(quicksum(x[j, i] for i in V if i != j) - (quicksum(x[i, j] for i in V if i != j)))
== 0
for j in V
)
# Eqtn 5
mdl.addConstrs(
(x[i, j] == 1) >> (u[j] <= (u[i] - (b[i] * x[i, j]) + C * (1 - x[i, j])))
for i in V
for j in V
if i != j
)
# Eqtn 6
mdl.addConstrs(u[i] >= 0 for i in V)
mdl.addConstrs(u[i] <= C for i in V)
# Eqtn 7
mdl.addConstrs(
(x[i, j] == 1) >> (y[j] <= (y[i] - (w[i, j] * x[i, j]) + Q * (1 - x[i, j])))
for i in I
for j in V
if i != j
)
# Eqtn 8
mdl.addConstrs(
(x[i, j] == 1) >> (y[j] <= (Q - (w[i, j] * x[i, j])))
for i in FuD
for j in V
if i != j
)
# Eqtn 9
mdl.addConstrs(y[i] >= 0 for i in V)
mdl.addConstrs(y[i] <= Q for i in V)
# Time to optimize is set really short since model is infeasible.
mdl.Params.TimeLimit = 10
mdl.optimize()
Now, when I run the above, I get the following error:
Set parameter TimeLimit to value 10
Gurobi Optimizer version 9.5.1 build v9.5.1rc2 (win64)
Thread count: 8 physical cores, 16 logical processors, using up to 16 threads
Optimize a model with 29 rows, 30 columns and 76 nonzeros
Model fingerprint: 0xb2ad0dc6
Model has 40 general constraints
Variable types: 10 continuous, 20 integer (20 binary)
Coefficient statistics:
Matrix range [1e+00, 1e+00]
Objective range [1e+00, 3e+01]
Bounds range [1e+00, 1e+00]
RHS range [1e+00, 2e+01]
GenCon rhs range [1e+01, 2e+01]
GenCon coe range [1e+00, 5e+01]
Presolve added 2 rows and 0 columns
Presolve removed 0 rows and 7 columns
Presolve time: 0.00s
Presolved: 31 rows, 23 columns, 110 nonzeros
Variable types: 7 continuous, 16 integer (16 binary)
Root relaxation: objective 6.800000e+00, 8 iterations, 0.00 seconds (0.00 work units)
Nodes | Current Node | Objective Bounds | Work
Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
0 0 6.80000 0 6 - 6.80000 - - 0s
0 0 7.00000 0 6 - 7.00000 - - 0s
Cutting planes:
Learned: 1
Gomory: 2
Implied bound: 4
MIR: 5
Explored 1 nodes (17 simplex iterations) in 0.01 seconds (0.00 work units)
Thread count was 16 (of 16 available processors)
Solution count 0
Model is infeasible
Best objective -, best bound -, gap -
Can someone please help me? I apologize if I made a rookie mistake, I have been doing a lot of research online but just can't put my finger on what the problem is.
Thank you in advance for taking the time to go through my entire post. I appreciate all the help I can get.
-
Hi Vanny!
I could not spot an obvious error in your formulation, except that you are using indicator constraints in (5), (7), and (8). I would also put the upper bounds on your capacity variables right into the variable declaration instead of adding them as constraints (6) and (9):
u = mdl.addVars(V, vtype=GRB.CONTINUOUS, ub=C, name="currentCap") # CURRENT CAPACITY
y = mdl.addVars(V, vtype=GRB.CONTINUOUS, ub=Q, name="batteryLevel") # CURRENT BATTERYI could not generate a feasible solution even when increasing the capacities and after removing the indicator restrictions.
You should compute the IIS of that model and inspect that to figure out what causes the infeasibility:
mdl.computeIIS()
mdl.write("model.ilp")I hope that helps!
Cheers,
Matthias0 -
Hello Matthias Miltenberger!
Thank you for your incredibly fast response. I have followed your suggestions, and I removed constraints 6 & 9 (and added the upper bounds as you explained). I have also removed the indicator constraints in (5), (7), (8).
I also ran the computeIIS() after running the mdl.optimize() command. I have been reading about IIS and as I understand, it tells us what constraints should be relaxed in order to make the model feasible. Am I correct? However, I have been trying to understand the output of it and I just can't wrap my head around it. Would you please let me know what this output means?:\ Model EVRP_copy
\ LP format - for model browsing. Use MPS format to capture full model detail.
Minimize
Subject To
R1: currentArc[4,0] + currentArc[4,1] + currentArc[4,2] + currentArc[4,3]
= 1
R4: currentArc[0,1] + currentArc[0,2] + currentArc[0,3] + currentArc[0,4]
- currentArc[1,0] - currentArc[2,0] - currentArc[3,0] - currentArc[4,0]
= 0
R5: - currentArc[0,1] + currentArc[1,0] + currentArc[1,2]
+ currentArc[1,3] + currentArc[1,4] - currentArc[2,1] - currentArc[3,1]
- currentArc[4,1] = 0
R6: - currentArc[0,2] - currentArc[1,2] + currentArc[2,0]
+ currentArc[2,1] + currentArc[2,3] + currentArc[2,4] - currentArc[3,2]
- currentArc[4,2] = 0
R8: - currentArc[0,4] - currentArc[1,4] - currentArc[2,4]
- currentArc[3,4] + currentArc[4,0] + currentArc[4,1] + currentArc[4,2]
+ currentArc[4,3] = 0
R10: 10 currentArc[0,2] - currentCap[0] + currentCap[2] <= 10
R11: 10 currentArc[0,3] - currentCap[0] + currentCap[3] <= 10
R13: 15 currentArc[1,0] + currentCap[0] - currentCap[1] <= 10
R14: 15 currentArc[1,2] - currentCap[1] + currentCap[2] <= 10
R15: 15 currentArc[1,3] - currentCap[1] + currentCap[3] <= 10
R16: 15 currentArc[1,4] - currentCap[1] + currentCap[4] <= 10
R17: 10 currentArc[2,0] + currentCap[0] - currentCap[2] <= 10
R18: 10 currentArc[2,1] + currentCap[1] - currentCap[2] <= 10
R19: 10 currentArc[2,3] - currentCap[2] + currentCap[3] <= 10
R20: 10 currentArc[2,4] - currentCap[2] + currentCap[4] <= 10
R21: 10 currentArc[3,0] + currentCap[0] - currentCap[3] <= 10
R22: 10 currentArc[3,1] + currentCap[1] - currentCap[3] <= 10
R23: 10 currentArc[3,2] + currentCap[2] - currentCap[3] <= 10
R24: 10 currentArc[3,4] - currentCap[3] + currentCap[4] <= 10
R26: 13 currentArc[4,1] + currentCap[1] - currentCap[4] <= 10
R27: 13 currentArc[4,2] + currentCap[2] - currentCap[4] <= 10
R28: 13 currentArc[4,3] + currentCap[3] - currentCap[4] <= 10
R33: 45 currentArc[4,0] + batteryLevel[0] - batteryLevel[4] <= 15
R37: 20 currentArc[0,1] + batteryLevel[1] <= 15
R40: 30 currentArc[0,4] + batteryLevel[4] <= 15
Bounds
currentCap[0] free
currentCap[1] free
currentCap[2] free
currentCap[3] free
currentCap[4] free
Binaries
currentArc[0,1] currentArc[0,2] currentArc[0,3] currentArc[0,4]
currentArc[1,0] currentArc[1,2] currentArc[1,3] currentArc[1,4]
currentArc[2,0] currentArc[2,1] currentArc[2,3] currentArc[2,4]
currentArc[3,0] currentArc[3,1] currentArc[3,2] currentArc[3,4]
currentArc[4,0] currentArc[4,1] currentArc[4,2] currentArc[4,3]
EndThank you once again! :)
0 -
Hi Vanny,
I think the infeasibility may be coming from the demand data and capacity constraint.
If we have positive demand for any node this is always infeasible unless one of the nodes has the same negative demand (supply). The original paper contains information about the initialisation of the variables, \(u_0=C, y_0=Q\), but I think this is still infeasible.Cheers,
David0 -
You can try replacing the capacity constraints with an equivalent constraint.
For example, you can use the \(\texttt{addVehicleLoadConstraintsByFlows}\) in the cvrp notebook:https://github.com/ruthmair/vrp/blob/main/cvrp.ipynb
I tried this and it worked as (remove your constraints 5 and 6 and replace them with)
z = mdl.addVars(A, lb=0, ub=C, name="z")
for i in I:
z[0, i].UB = 0
mdl.addConstrs(
(z.sum("*", j) + b[j] == z.sum(j, "*") for j in I),
name="flowConservation",
)
mdl.addConstrs(
(z[i, j] >= b[i] * x[i, j] for i in I for j in V if i != j),
name="loadLowerBound",
)
mdl.addConstrs(
(z[i, j] <= (C - b[j]) * x[i, j] for i in I for j in V if i != j),
name="loadUpperBound",
)Cheers,
David0 -
Hello Dear David,
First of all, thank you for taking the time to explain the code and suggest a solution.
I would like to let you know that I have implemented your suggestions, and even though the model stops being infeasible, the solutions are not correct because of the following reasons:
1. The vehicles should depart from and return to the depot (node 0), but they leave from and return to node 1 (which is a customer).
2. Taking a look at the variables of the solution obtained, the battery information when the vehicle is at the nodes it visited is incorrect (I know this because I generated this very small dataset to test the model). You can see it is incorrect as the optimal solution should be [0, 2, 1, 4, 3, 0] --> but gurobi says it is [1,3,4,1].
Here is an image showing what the actual values and solution should be:
But when I take a look at my variables I see the following:
<gurobi.Var batteryLevel[0] (value 0.0)> <gurobi.Var batteryLevel[1] (value 1.0)> <gurobi.Var batteryLevel[2] (value 0.0)> <gurobi.Var batteryLevel[3] (value 0.0)> <gurobi.Var batteryLevel[4] (value 4.0)>
Which from what I understand is completely wrong based on the formula for battery consumption.
Here is the code I ran to get such output:from gurobipy import Model, GRB, quicksum
mdl = Model('EVRP')
x = mdl.addVars(A, vtype=GRB.BINARY, name="currentArc") # CURRENT ARC???
u = mdl.addVars(V, vtype=GRB.CONTINUOUS, ub=C, name="currentCap") # CURRENT CAPACITY
y = mdl.addVars(V, vtype=GRB.CONTINUOUS, ub=Q, name="batteryLevel") # CURRENT BATTERY
z = mdl.addVars(A, lb=0, ub=C, name="z")
# Eqtn 5
for i in I:
z[0, i].UB = 0
# Eqtn 1 --> Objective Function
mdl.modelSense = GRB.MINIMIZE
mdl.setObjective(quicksum(w[i,j]*x[i,j] for i,j in A))
# Eqtn 2
mdl.addConstrs(quicksum(x[i,j] for j in V if j!=i)==1 for i in I)
# Eqtn 3
mdl.addConstrs(quicksum(x[i,j] for j in V if j!=i)<=1 for i in F)
# Eqtn 4
mdl.addConstrs((quicksum(x[j,i] for i in V if i!=j)-(quicksum(x[i,j] for i in V if i!=j)))==0 for j in V)
mdl.addConstrs(
(z.sum("*", j) + b[j] == z.sum(j, "*") for j in I),
name="flowConservation",
)
mdl.addConstrs(
(z[i, j] >= b[i] * x[i, j] for i in I for j in V if i != j),
name="loadLowerBound",
)
mdl.addConstrs(
(z[i, j] <= (C - b[j]) * x[i, j] for i in I for j in V if i != j),
name="loadUpperBound",
)
# Eqtn 7
mdl.addConstrs((y[j] <= (y[i]-(w[i,j]*x[i,j])+Q*(1-x[i,j]))) for i in I for j in V if i!=j)
# Eqtn 8
mdl.addConstrs((y[j] <= (Q-(w[i,j]*x[i,j]))) for i in FuD for j in V if i!=j)Again, thank you for all your time and patience. I would really appreciate if you can take a look at the code and suggest what could be going on so that my vehicle ignores the depot and the battery capacity is wrong.
Thank you!0
Please sign in to leave a comment.
Comments
5 comments