Skip to main content

The Electric Vehicle Routing Problem (EVRP) Model is always infeasible

Answered

Comments

5 comments

  • Matthias Miltenberger
    Gurobi Staff 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 CAPACITY
    y = 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

    0
  • vanny minanda
    First Comment
    First Question

    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]
    End

    Thank you once again! :)

    0
  • David Torres Sanchez
    Gurobi Staff 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

    0
  • David Torres Sanchez
    Gurobi Staff 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 = 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, 
    David

    0
  • vanny minanda
    First Comment
    First Question

    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.