• Gurobi Staff

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 CAPACITYy = mdl.addVars(V, vtype=GRB.CONTINUOUS, ub=Q, name="batteryLevel")  # CURRENT BATTERY

I 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,
Matthias

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] <= 15Bounds currentCap[0] free currentCap[1] free currentCap[2] free currentCap[3] free currentCap[4] freeBinaries 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]End

Thank you once again! :)

• Gurobi Staff

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,
David

• Gurobi Staff

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 = 0mdl.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,
David

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, quicksummdl = 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 CAPACITYy = mdl.addVars(V, vtype=GRB.CONTINUOUS, ub=Q, name="batteryLevel")    # CURRENT BATTERYz = mdl.addVars(A, lb=0, ub=C, name="z")# Eqtn 5for i in I:    z[0, i].UB = 0# Eqtn 1 --> Objective Functionmdl.modelSense = GRB.MINIMIZEmdl.setObjective(quicksum(w[i,j]*x[i,j] for i,j in A))# Eqtn 2mdl.addConstrs(quicksum(x[i,j] for j in V if j!=i)==1 for i in I)# Eqtn 3mdl.addConstrs(quicksum(x[i,j] for j in V if j!=i)<=1 for i in F)# Eqtn 4mdl.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 7mdl.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 8mdl.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!