VRP Multi-Period Model infeasible
Answered
Hi Gurobi Team,
I'm working on my undergraduate thesis about Multi Period-Capacitated-Time Windows-VRP to minimize time travel. I'm working on scheduling and pickup routing for goods. It's said my model is infeasible. I have tried using IIS, and it's shows several decision variable that error, but I still can't figure out why the value doesn't want to come out. Also, I can't get attributes for decision variable e.
For this model, the decision variables are:
- Schedule in period t (k index) with vehicles i (y[i,j,k]
- Routing in period t (k index) with vehicles i (z[i,j,l,k]
- Pickup quantity to satisfy the demand (e[i, j, k])
- Service Start Time for each scheduling. (s[i,j,k])
Here are the decision variables:
def decisionVariables(M, gb, m, n, t, e, I, y, z, s):
'''decision variable y'''
for i in range(m):
for j in range(n + 1): # +1 is for depot
for k in range(t):
y[i, j, k] = M.addVar(vtype=gb.GRB.BINARY, name='y%d,%d,%d' % (i, j, k))
M.update()
'''variable s'''
for i in range(m):
for j in range(n + 1):
for k in range(t):
s[i, j, k] = M.addVar(vtype=gb.GRB.CONTINUOUS, name="s%d,%d,%d" % (i, j, k))
M.update()
'''variable z'''
for i in range(m):
for j in range(n + 1):
for l in range(n + 1):
for k in range(t):
z[i, j, l, k] = M.addVar(vtype=gb.GRB.BINARY, name="z%d,%d,%d,%d" % (i, j, l, k))
M.update()
"variable e"
for i in range (m):
for j in range (n): #nanti cek dulu n or n+1
for k in range (t):
e[i, j, k] = M.addVar(vtype=gb.GRB.INTEGER, name="e%d,%d,%d" % (i, j, k))
M.update()
"variable I"
for k in range (t):
I[k] = M.addVar(vtype=gb.GRB.INTEGER, name="I%d" % (k))
M.update()
Here are my constraints:
def Constraints(M, gb, e, I, y, z, s, m, n, t, et, lt, sd, grid, bigM, d, P, initial_inventory, ss, df_V, tao):
for j in range(n):
'''constraint 5.5: total number of visit'''
M.addConstr((gb.quicksum(gb.quicksum(y[i, j + 1, k] for i in range(m)) for k in range(t)) == 1), name = "C.55")
M.update()
'''constraint 5.9: every vehicle start from depot'''
for i in range(m):
for k in range(t):
M.addConstr((y[i, 0, k] == 1), name="C.9") # j = 0 for depot location
M.update()
for i in range(m):
for j in range(n + 1):
for k in range(t):
if j != 0:
'''constraint 5.11 and 5.12: every node visited once'''
M.addConstr((gb.quicksum(z[i, j, l, k] for l in range(n + 1)) == y[i, j, k]), name="C.5.11")
M.addConstr((gb.quicksum(z[i, l, j, k] for l in range(n + 1)) == y[i, j, k]), name = "C.5.12")
if j != 0:
'''constraint 5.13: ommiting a loop inside a node'''
M.addConstr((z[i, j, j, k] == 0), name="C.513")
M.update()
for i in range(m):
for j in range(n + 1):
for k in range(t):
'''constraint 5.14 and 5.15: time window'''
M.addConstr((s[i, j, k] + bigM *(1 - y[i, j, k]) >= et[j] + (grid[0, j] * z[i, 0, j, k])),name="C.14")
M.addConstr((s[i, j, k] + sd[j] - bigM *(1 - y[i, j, k]) + (grid[j, 0] * z[i, j, 0, k]) <= lt[j]),name="C.15")
for i in range(m):
for j in range(n+1):
for k in range(t):
'''constraint 5.16. service start time'''
M.addConstr((s[i, j, k] <= lt[j] * y[i, j, k]), name="5.16")
M.update()
for i in range(m):
for k in range(t):
for j in range(n):
for l in range(n):
'''constraint 5.17. scheduling constraint'''
M.addConstr((s[i, j+1, k] + (sd[j+1] + grid[j+1, l+1]) * z[i, j+1, l+1, k] <= s[i, l+1, k] + bigM * (1 - z[i, j+1, l+1, k])), name="5.17")
M.update()
"constraint pickup has to be equal as the nodes' produce"
for k in range (t):
for j in range (n):
M.addConstr(((gb.quicksum(y[i, j, k] * P[k] for i in range(m))) == gb.quicksum(e[i, j, k] for i in range (m))), name = "pickup-produce")
M.update
"dummy pickup minimal qty for checking the model"
for k in range (t):
M.addConstr((gb.quicksum(gb.quicksum(e[i, j, k] for i in range (m)) for j in range (n)) == d[k]), name="dummy pickup")
M.update()
'''constraint 5.18: vehicles' maximum capacity'''
for i in range(m):
for k in range(t):
M.addConstr((gb.quicksum(e[i, j, k] for j in range (n)) <= df_V["C"][i]), name="5.18")
M.update()
# '''constraint 5.19: vehicles' maximal total travelling time'''
for i in range(m):
M.addConstr((gb.quicksum(gb.quicksum((y[i, j + 1, k] * sd[j+1]) for j in range (n)) for k in range(t)) + gb.quicksum(gb.quicksum(gb.quicksum((grid[j,l] * z[i, j, l, k]) for j in range (n+1)) for l in range (n+1)) for k in range (t)) <= tao), name= "spoilage")
M.update()
# "constraint pickup minimal quantity - considering multi period inventory"
# for k in range (t+1):
# if k == 0:
# M.addConstr(gb.quicksum(gb.quicksum(e[i, j, k] for i in range (m)) for j in range (n)) >= initial_inventory - d[k] - ss) #apakah ga masalah d[k] aja?
# if k <= t:
# M.addConstr(gb.quicksum(gb.quicksum(e[i, j, k+1] for i in range (m)) for j in range (n)) >= I[k] - d[k+1] - ss) #not really sure penggunaannya
# M.update()
"constraint inventory"
for k in range (t):
if k == 0:
M.addConstr((I[k] == initial_inventory - d[k] - gb.quicksum(gb.quicksum(e[i, j, k] for i in range (m)) for j in range (n))),name="inventory 1")
else:
M.addConstr((I[k] == I[k-1] - d[k] - gb.quicksum(gb.quicksum(e[i, j, k] for i in range (m)) for j in range (n))), name= "inventory 2")
M.update()
Then it shows that my program is infeasible. Then, I used IIS for finding the error and I got this:
Computing Irreducible Inconsistent Subsystem (IIS)...
Constraints | Bounds | Runtime
Min Max Guess | Min Max Guess |
--------------------------------------------------------------------------
0 261 - 0 46 - 0s
4 4 - 1 1 - 0s
IIS computed: 4 constraints, 1 bounds
IIS runtime: 0.01 seconds
Gurobi Optimizer version 9.1.0 build v9.1.0rc0 (win64)
Thread count: 4 physical cores, 8 logical processors, using up to 8 threads
Optimize a model with 261 rows, 214 columns and 1033 nonzeros
Model fingerprint: 0x03bd8d7c
Variable types: 24 continuous, 190 integer (168 binary)
Coefficient statistics:
Matrix range [1e+00, 2e+04]
Objective range [4e+00, 1e+04]
Bounds range [1e+00, 1e+00]
RHS range [1e+00, 1e+04]
Presolve removed 4 rows and 5 columns
Presolve time: 0.00s
Explored 0 nodes (0 simplex iterations) in 0.01 seconds
Thread count was 1 (of 8 available processors)
Solution count 0
Model is infeasible
Best objective -, best bound -, gap -
When I check the iis result, it shows that there are errors in the decision variable for demand satisfaction. It doesn't want to create value of it.
\ Model problem_copy
\ LP format - for model browsing. Use MPS format to capture full model detail.
Minimize
Subject To
dummy_pickup: e0,0,0 + e0,1,0 + e0,2,0 + e0,3,0 + e0,4,0 + e1,0,0 + e1,1,0
+ e1,2,0 + e1,3,0 + e1,4,0 = 54
dummy_pickup: e0,0,1 + e0,1,1 + e0,2,1 + e0,3,1 + e0,4,1 + e1,0,1 + e1,1,1
+ e1,2,1 + e1,3,1 + e1,4,1 = 99
inventory_1: e0,0,0 + e0,1,0 + e0,2,0 + e0,3,0 + e0,4,0 + e1,0,0 + e1,1,0
+ e1,2,0 + e1,3,0 + e1,4,0 + I0 = 6
inventory_2: e0,0,1 + e0,1,1 + e0,2,1 + e0,3,1 + e0,4,1 + e1,0,1 + e1,1,1
+ e1,2,1 + e1,3,1 + e1,4,1 - I0 + I1 = -99
Bounds
e0,0,0 free
e0,0,1 free
e0,1,0 free
e0,1,1 free
e0,2,0 free
e0,2,1 free
e0,3,0 free
e0,3,1 free
e0,4,0 free
e0,4,1 free
e1,0,0 free
e1,0,1 free
e1,1,0 free
e1,1,1 free
e1,2,0 free
e1,2,1 free
e1,3,0 free
e1,3,1 free
e1,4,0 free
e1,4,1 free
I0 free
Generals
e0,0,0 e0,0,1 e0,1,0 e0,1,1 e0,2,0 e0,2,1 e0,3,0 e0,3,1 e0,4,0 e0,4,1
e1,0,0 e1,0,1 e1,1,0 e1,1,1 e1,2,0 e1,2,1 e1,3,0 e1,3,1 e1,4,0 e1,4,1 I0
I1
End
Could you please help to check where the possible error is?
Thank you
-
Hi Maura,
To simplify your IIS I will introduce the following
le1 = e0,0,0 + e0,1,0 + e0,2,0 + e0,3,0 + e0,4,0 + e1,0,0 + e1,1,0
+ e1,2,0 + e1,3,0 + e1,4,0
le2 = e0,0,1 + e0,1,1 + e0,2,1 + e0,3,1 + e0,4,1 + e1,0,1 + e1,1,1
+ e1,2,1 + e1,3,1 + e1,4,1The constraints in your IIS can then be written as
dummy_pickup: le1 = 54
dummy_pickup: le2 = 99
inventory_1: le1 + I0 = 6
inventory_2: le2 - I0 + I1 = -99From here it is easy to deduce that
I0 = 6 - 54 = -48
I1 = -99 + I0 - le2 = -246The problem is I1 looks to have a lower bound of 0 since it is not listed in the "free" variables (the default lower bound is 0). But if its lower bound is 0 then it cannot have a value of -246.
At quick glance I'm wondering if your inventory constraints should be more like this
M.addConstr((I[k] == I[k-1] - d[k] + gb.quicksum(e[i, j, k] for i in range (m) for j in range (n))), name= "inventory 2")
- Riley
0 -
Hi Riley,
Thank you for the heads up! Yes you are right, I did some mistake while writing the inventory constraints. After giving revision of it, it still giving infeasible result and showing this result with IIS:
\ Model problem_copy
\ LP format - for model browsing. Use MPS format to capture full model detail.
Minimize
Subject To
pickup-produce: - e0,0,0 - e1,0,0 = 0
pickup-produce: - e0,1,0 - e1,1,0 = 0
pickup-produce: - e0,2,0 - e1,2,0 = 0
pickup-produce: - e0,3,0 - e1,3,0 = 0
pickup-produce: - e0,4,0 - e1,4,0 = 0
dummy_pickup: e0,0,0 + e0,1,0 + e0,2,0 + e0,3,0 + e0,4,0 + e1,0,0 + e1,1,0
+ e1,2,0 + e1,3,0 + e1,4,0 = 54
Bounds
e0,0,0 free
e0,1,0 free
e0,2,0 free
e0,3,0 free
e0,4,0 free
e1,0,0 free
e1,1,0 free
e1,2,0 free
e1,3,0 free
e1,4,0 free
Generals
e0,0,0 e0,1,0 e0,2,0 e0,3,0 e0,4,0 e1,0,0 e1,1,0 e1,2,0 e1,3,0 e1,4,0
EndI do not understand why the pickup-produce is still zero. Why does the y decision variable doesn't multiply the production qty? Could you please give me clue on that?
Also, I still find the "e" attribute error:
AttributeError Traceback (most recent call last) <ipython-input-7-54227bb1f309> in <module> 12 13 # get attributes ---> 14 sol_e, sol_I, sol_y, sol_z, sol_s = M.getAttr('e', e), M.getAttr('I', I), M.getAttr('y', y), M.getAttr('z', z), M.getAttr('s', s) # get decision variables 15 16 # print_solution.printSolution(M, sol_x, sol_y, sol_z, sol_s, m, n, t, q, f, df, df_n) # print solution src\gurobipy\model.pxi in gurobipy.gurobipy.Model.getAttr() src\gurobipy\attrutil.pxi in gurobipy.gurobipy.__getattrinfo() AttributeError: 'gurobipy.Model' object has no attribute 'e'
Meanwhile, I have written the e variable in the decision variable function:
"variable e"
for i in range (m):
for j in range (n): #nanti cek dulu n or n+1
for k in range (t):
e[i, j, k] = M.addVar(vtype=gb.GRB.INTEGER, name="e%d,%d,%d" % (i, j, k))
M.update()And call it in the same syntax:
decisionVariables(M, gb, m, n, t, e, I, y, z, s) # add decision variables to model
objectiveFunction(M, gb, z, n, m, t, grid) # add objective function to model
Constraints(M, gb, e, I, y, z, s, m, n, t, et, lt, sd, grid, bigM, d, P, initial_inventory, ss, df_V, tao) # add constraints to model
M.computeIIS()
M.write('iismodel.ilp')
# get attributes
sol_e, sol_I, sol_y, sol_z, sol_s = M.getAttr('e', e), M.getAttr('I', I), M.getAttr('y', y), M.getAttr('z', z), M.getAttr('s', s) # get decision variablesCould you please help to check on it also where the possible error will be?
Thank you very much for the kind help
0 -
For giving more context,
In the pickup-produce constraint, I would like to create a scenario where everytime the vehicles visit that node in specific time (node can only be visited once for all period and by any vehicle), it will pickup all the quantities that being produced by the node.
0 -
Hi Maura,
I would suggest running the following code, based off the pickup-produce constraint, and seeing if it is what you expect it to be
for k in range (t):
for j in range (n):
print(f"k={k}", f"j={j}", f"P[k]={P[k]}", f"empty LHS expected={P[k]==0 or m==0}")- Riley
0 -
Hi Riley,
Thank you for your help, after changing the pickup-produce constraint as you suggested, the model is not infeasible anymore. But I still have this difficulty in the attribute e error. Could you help to check on this?
AttributeError Traceback (most recent call last) <ipython-input-20-df5f488c5933> in <module> 12 13 # get attributes ---> 14 sol_e, sol_I, sol_y, sol_z, sol_s = M.getAttr('e', e), M.getAttr('I', I), M.getAttr('y', y), M.getAttr('z', z), M.getAttr('s', s) # get decision variables 15 16 # print_solution.printSolution(M, sol_x, sol_y, sol_z, sol_s, m, n, t, q, f, df, df_n) # print solution src\gurobipy\model.pxi in gurobipy.gurobipy.Model.getAttr() src\gurobipy\attrutil.pxi in gurobipy.gurobipy.__getattrinfo() AttributeError: 'gurobipy.Model' object has no attribute 'e'
As I mentioned in the earlier message, I have define the attribute 'e' as the decision variable. Could you help to check where the possible error is?
Thank you very much. Looking forward for your response
0 -
Hi Maura,
Solution values are accessed from the variables, not the model object, so not sure what has led you down this path. If you have a variable v then you can obtain its value in the optimal solution using v.X
What is the "e" parameter that you are passing into your function? If it is a dictionary then you could create another dictionary of solution values like this
sol_e = {k:v.X for k,v in e.items()}or to create a dictionary of all variables:
sol = {v.VarName:v.X for v in M.getVars()}
- Riley1 -
Hi Riley,
If it's like that, so what's the '.X' define in the code?
Also, the constraint that you created previously for decision variable "e" is not creating the boundaries that need to be fulfilled for the model.
The expected coding is to fulfill this constraint:
So I code:
for j in range (n):
M.addConstr(((gb.quicksum(gb.quicksum(e[i, j, k] for i in range (m)) for k in range (t))) == P[j]), name = "pickup qty = production")
M.update
for i in range (m):
for j in range (n):
for k in range (t):
M.updateThen it still shows the result is infeasible, and this is that shows in the IIS result
Subject To
pickup_qty_=_production: e0,4,0 + e0,4,1 + e1,4,0 + e1,4,1 = 80
pickup_qty_definition: - e0,4,0 = 0
pickup_qty_definition: 30 y0,4,1 - e0,4,1 = 0
pickup_qty_definition: - e1,4,0 = 0
pickup_qty_definition: 30 y1,4,1 - e1,4,1 = 0
Bounds
e0,4,0 free
e0,4,1 free
e1,4,0 free
e1,4,1 free
Binaries
y0,4,1 y1,4,1
Generals
e0,4,0 e0,4,1 e1,4,0 e1,4,1
EndCould you please help in this issue again? But to make sure the constraint is fulfilled, as the previous constraint doesn't give the optimum result.
Thank you very much
0 -
Hi Maura,
The "X" attribute of a variable is its value in the solution. You can find a list of all attributes available here.
If you take the pickup_qty_=_production constraint in the IIS and use the 4 pickup_qty_definition constraints underneath it to perform substitutions you can arrive at the following equation:
30 y0,4,1 + 30 y1,4,1 = 80
We obviously can't satisfy this when the y variables are binary, which should give you a hint as to why the model is infeasible. I can't see a definition for the pickup_qty_definition constraints in your code and so it is hard to give more insight than this
Also, the constraint that you created previously for decision variable "e" is not creating the boundaries that need to be fulfilled for the model.
I'm not sure what this means, but let's address the infeasibility first then return to this.
- Riley
0 -
Hi Riley, noted. Thank you very much for your help on this, I can understand where the error within my code. After giving several revisions, the code looks like this. I also put several constraint that I didn't put before which related to variable "e".
for j in range (n):
M.addConstr(((gb.quicksum(gb.quicksum(e[i, j, k] for i in range (m)) for k in range (t))) == P[j]), name = "pickup qty = production")
M.update
for i in range (m):
for j in range (n):
for k in range (t):
M.addConstr((y[i, j, k] * P[k] == e[i, j, k]), name= "pickup qty definition")
M.update
"dummy pickup minimal qty for checking the model" #--> ada update
for k in range (t):
M.addConstr((gb.quicksum(gb.quicksum(e[i, j, k] for i in range (m)) for j in range (n)) >= d[k]), name="dummy pickup")
M.update()
'''constraint 5.18: vehicles' maximum capacity'''
for i in range(m):
for k in range(t):
M.addConstr((gb.quicksum(e[i, j, k] for j in range (n)) <= df_V["C"][i]), name="5.18")
M.update()Then I checked my solution and get this result in IIS.
Subject To
C.5.9: y1,0,0 = 1
C.5.9: y1,0,1 = 1
pickup_qty_=_production: e0,0,0 + e0,0,1 + e1,0,0 + e1,0,1 = 30
pickup_qty_definition: 30 y0,0,0 - e0,0,0 = 0
pickup_qty_definition: 30 y0,0,1 - e0,0,1 = 0
pickup_qty_definition: 30 y1,0,0 - e1,0,0 = 0
pickup_qty_definition: 30 y1,0,1 - e1,0,1 = 0
Bounds
e0,0,0 free
e0,0,1 free
e1,0,0 free
e1,0,1 free
Binaries
y0,0,0 y0,0,1 y1,0,0 y1,0,1
Generals
e0,0,0 e0,0,1 e1,0,0 e1,0,1
EndFor context, this is C.5.9. Each period, the vehicle needs to start from depot (node j=0) to assigned node (but still, each node can only be visited once in all period). Variable 'y' consider depot in the decision variable, hence 0 is depot. Variable 'e' doesn't consider depot, hence 0 is the 1st node that needs to be visited.
'''constraint 5.9: every vehicle start from depot'''
for i in range(m):
for k in range(t):
M.addConstr((y[i, 0, k] == 1), name="C.5.9") # j = 0 for depot location
M.update()For your information:
- i = index for vehicle
- j = index for nodes
- k = index for time
Could you help to assist with the possible error cause?
Thank you very much, very appreciate it.
0 -
Hi Maura,
Can you explain to me exactly what it means when the variable y[a,b,c] = 1? Something to do with scheduling vehicles?
- Riley
0 -
Yes, it shows the schedule for that nodes and vehicle in specific timing period. If it's giving y[i,j,k] = 1, it means that vehicle i is visiting node j in period k. Hence, the variable 'e' shows the pickup qty that needs to be taken in every visit, in this case, all the qty that is produced by node j. What limits the schedule will be the capacity of the vehicle (because it's a single visit) and the time window.
0 -
Hi Maura,
In that case these constraints:
'''constraint 5.9: every vehicle start from depot'''
for i in range(m):
for k in range(t):
M.addConstr((y[i, 0, k] == 1), name="C.5.9") # j = 0 for depot locationare saying that every vehicle should visit the depot in every period. This surely is incorrect, and likely one source of your infeasibility.
- Riley
0 -
Ah i see, i think that's the error comes from. So I have set of n indexed with j. The length is 6, with 0 is depot and 1-5 is the node that needs to be visited. Since generally we use lots of 1-5 data than 0-5, I define n as n-1.
In variable 'e' it uses
for j in range (n):
Hence it will use the data 1-5 only. So the e[0,0,0] with j=0 will refer to the 1st node that needs to be visited.
In variable 'y' it uses
for j in range (n+1):
The context in here is to also schedule the depot. So, the y[0,0,0] with j=0 will refer to the depot (the snippet that you just sent is meant to define that depot is always visited daily, but the other nodes are only visited once in all period).
I think I get the error from here where the code can't define which '0' should it refer to. Could you help to give me suggestion regarding this matter? To give distinction when trying to call 0-5 and 1-5?
Thank you!
0 -
Hi Maura,
To get an iterator for 0-5 you use range(6). To get an iterator from 1-5 you can use range(1,6). If you want the integer sequence a, a+1, ..., b-1, b then you can use range(a, b+1).
When the sequence is small you can also just explicitly write it, eg
for i in (1,2,3,4,5):
...- Riley
0
Please sign in to leave a comment.
Comments
14 comments