Adding indicator constraints giving 'model infeasible solution' output
AnsweredHello Everyone,
I am a new user of Gurobi and I am getting used to modeling in it. Right now I have a multi-trip (reloads up to vehicle's capacity from the depot whenever the vehicle needs to) capacitated VRP and I am modeling it in Gurobi. I am able to model the necessary constraints required for the multi-trip CVRP problem in the Gurobi and I am able to get some output optimal solution. But when I add additional indicator constraint, the model throws an error saying that the model is infeasible. The indicator constraint that I am planning to add is, when the vehicle's capacity is more than some threshold limit, then it can drop the reloading depot. Only when the vehicle capacity is less or equal to that threshold limit, then it can go to the reloading depot else, vehicle can drop it. I will send you the model code that I made for this problem. I will highlight/bold the portion where I used the indicator constraints below. Can anyone help me find what is the issue that I am facing? Looking forward to hearing from you. Thank you.
import pandas as pd
import numpy as np
from gurobipy import *
Q = 30000 # vehicle capacity
no_of_vehicles = 2 # free to modify
vehicle_speed = 33
df = pd.read_excel('VRP_reload_data.xlsx', engine='openpyxl')
Y, X = list(df["Y"]), list(df["X"])
coordinates = np.column_stack((X, Y))
et, lt, st = list(df["St"]), list(df["Et"]), list(df["Sert"])
n = len(coordinates)
depot, customers = coordinates[0:9, :], coordinates[9:, :]
M = 100**100 # big number
m = Model("MTVRPCTW")
x, f, t = {}, {}, {} #intialize the decision variables
'''distance matrix (34*34 array)'''
dist_matrix = np.empty([n, n])
for i in range(len(X)):
for j in range(len(Y)):
'''variable_1: X[i,j] =(0,1), i,j = Nodes'''
x[i, j] = m.addVar(vtype=GRB.BINARY, name="x%d,%d" % (i, j))
dist_matrix[i, j] = np.sqrt((X[i] - X[j]) ** 2 + (Y[i] - Y[j]) ** 2)
if i == j:
dist_matrix[i, j] = M ## big 'M'
continue
m.update()
'''variable_2: f[j] = cumulative demand'''
for j in range(n):
f[j] = m.addVar(lb=0, vtype=GRB.INTEGER, name="y%d" % (j)) # cumulative demand satisfied variable
t[j] = m.addVar(lb=0, vtype=GRB.INTEGER, name="z%d" % (j)) # cumulative time variable
m.update()
# vehicles leaving each node except start and end depot
for i in range(n):
if i == 0 or i == 8:
continue
m.addConstr(quicksum(x[(i, j)] for j in range(n)) == 1)
m.update()
# vehicles arriving to each node except start and end depot
for j in range(n):
if j == 0 or j == 8:
continue
m.addConstr(quicksum(x[(i, j)] for i in range(n)) == 1)
m.update()
# vehicles leaving start depot
m.addConstr(quicksum(x[(0, j)] for j in range(n)) == no_of_vehicles)
m.update()
# vehicles arriving to end depot
m.addConstr(quicksum(x[(i, 8)] for i in range(n)) == no_of_vehicles)
m.update()
'''constraint_5: capacity of vehicle, if the vehicle visits any node between 0 through 8 (station nodes), it reloads to Q,
else vehicle continues to unload'''
for i in range(n):
for j in range(0, n):
if j <= 8:
m.addConstr(f[j] == Q)
m.addConstr(f[j] <= Q - (dist_matrix[i, j] * (x[i, j])) + M*(1-x[(i, j)]))
else:
m.addConstr(f[j] <= Q)
m.addConstr(f[j] >= 0)
m.addConstr(f[j] <= f[i] - (dist_matrix[i, j] * (x[i, j])) + M*(1-x[(i, j)]))
m.update()
'''constraint_6: time-windows and also eliminating sub-tours'''
for i in range(n):
for j in range(n):
m.addConstr(t[j] >= t[i] + ((st[i] + (dist_matrix[i, j]/vehicle_speed)) * x[(i, j)]) - M*(1-x[(i, j)]))
m.addConstr(t[j] >= (et[j])) # service should start after the earliest service start time
m.addConstr(t[j] <= (lt[j])) # service can't be started after the latest service start time
m.update()
'''constraint 7: dropping visits constraints'''
for i in range(n):
for j in range(9):
if j == 0 or j == 8:
continue
m.addConstr((f[i] >= Q//2) >> x[(i, j)] == 0)
'''objective function'''
m.setObjective(quicksum(quicksum(quicksum(x[(i, j)]*dist_matrix[(i, j)] for k in range(no_of_vehicles))
for i in range(n)) for j in range(n)), GRB.MINIMIZE)
m.update()
'''optimize'''
m.optimize()
'''retrieve the solution'''
sol_y, sol_x, sol_z = m.getAttr('x', f), m.getAttr('x', x), m.getAttr('x', t)
X, Y, Z = np.empty([n, n]), np.empty([n]), np.empty([n])
for i in range(n):
Y[i] = sol_y[i]
Z[i] = sol_z[i]
for j in range(n):
X[i, j] = int(sol_x[i, j])
print('\nObjective is:', m.objVal)
print('\nDecision variable X (binary decision of travelling from one node to another):\n', X.astype('int32'))
# print('\nDecision variable z:(service start time of every customers in minutes)\n', Z.astype('int32')[1:])
print('\nDecision variable y (cumulative demand collected at every customer node):\n', Y.astype('int32')[1:])
-
Official comment
-
Gurobi Staff
This post is more than three years old. Some information may not be up to date. For current information, please check the Gurobi Documentation or Knowledge Base. If you need more help, please create a new post in the community forum. Or why not try our AI Gurobot?. -
-
-
Hi,
You could try to find out what causes the infeasibility as described in the Knowledge Base article How do I determine why my model is infeasible?
Best regards,
Jaromił0 -
-
Hi Mr. Najman,
I tried finding out the cause of infeasibility as described in that article through model.computeIIS() method and it seems that there is one constraint that is getting violated. But when I tried to create a feasibillity relaxation using model.feasrelax() method, I am not getting an output at all and I wonder what is happening as I am new to Gurobi. It would be great if you assist me with this issue of how to find the infeasibility in my model or how to create feasibility relaxation. I will attach the snippet of the feasibility relaxation that I tried. Looking forward to hearing from you. Thanks.
if m.status == GRB.INFEASIBLE:
var = m.getVars()
lbpen = [1.0]*m.numVars
ubpen = [10.0]*m.numVars
constr = []
for i in range(n):
for j in range(9):
if j == 0 or j == 8:
continue
# if f[i] >= Q//2:
# x1 = 1
# m.addGenConstrIndicator(x1, True, x[(i, j)] == 0)
constr.append(m.addConstr((f[i] >= Q//2) >> x[(i, j)] == 0))
m.feasRelax(0, False, var, lbpen, ubpen, constr, 1.0)
m.optimize()0 -
-
-
Hi,
Do you get any error message? The \(\texttt{feasRelax}\) line should read
rhspen = [1.0]*len(constr)
m.feasRelax(0, False, var, lbpen, ubpen, constr, rhspen)Best regards,
Jaromił0 -
-
Hi,
I tried this way that you mentioned in the previous message.
if m.status == GRB.INFEASIBLE:
var = m.getVars()
lbpen = [1.0]*m.numVars
ubpen = [10.0]*m.numVars
constr = []
for i in range(n):
for j in range(9):
if j == 0 or j == 8:
continue
# if f[i] >= Q//2:
# x1 = 1
# m.addGenConstrIndicator(x1, True, x[(i, j)] == 0)
constr.append(m.addConstr((f[i] >= Q//2) >> x[(i, j)] == 0))
rhspen = [1.0]*len(constr)
m.feasRelax(0, False, var, lbpen, ubpen, constr, rhspen)
m.optimize()And the output error message that I got is
Traceback (most recent call last):
line 121, in <module>
sol_y, sol_x, sol_z = m.getAttr('x', f), m.getAttr('x', x), m.getAttr('x', t)
File "src\gurobipy\model.pxi", line 1848, in gurobipy.Model.getAttr
File "src\gurobipy\attrutil.pxi", line 148, in gurobipy.__gettypedattrlist
gurobipy.GurobiError: Unable to retrieve attribute 'x'This is the same error I got when my model was infeasible, except that the gurobi output said my model is infeasible along with this error message. But right now, only this error message is appearing. Since this is the same error message I got before applying m.feasRelax() method, I am not sure of what this error message is conveying.
0 -
-
-
Hi,
The error message means that no solution point value can be retrieved. When do you execute the line
sol_y, sol_x, sol_z = m.getAttr('x', f), m.getAttr('x', x), m.getAttr('x', t)
Could you post the output before the error occurs? You could rewrite your code as
sol_y, sol_x, sol_z = {}, {}, {}
if m.status == GRB.OPTIMAL:
sol_y, sol_x, sol_z = m.getAttr('x', f), m.getAttr('x', x), m.getAttr('x', t)to avoid this error and only retrieve the solution point values when the optimization status is optimal. Please also note that it is not possible for the Community to reproduce this issue because you use excel input data.
Best regards,
Jaromił0 -
-
Hi,
Thank you for your reply. Below here is the excel input data that I use. Kindly use this to reproduce the issue. You can find the code example above in the first message of this issue. I am unable to find the reason for infeasibility. Although, your suggestion helped me to get rid of that error message without any output (as the status is not optimal in my case).
Node X Y St Et Sert 1 13200 13200 0 60 300 2 15190.56 13543.2 660 960 300 3 17181.12 13886.4 1560 1860 300 4 19166.4 14224.32 2460 2760 300 5 21172.8 14572.8 3360 3660 300 6 17101.92 13543.2 4260 4560 300 7 13041.6 12518.88 5160 5460 300 8 8976 11510.4 6060 6360 300 9 4915.68 10507.2 6960 7260 300 10 6864 5808 0 100000 300 11 3168 22176 7 100000 300 12 24288 21120 10 100000 300 13 6864 22176 14 100000 300 14 22704 2640 0 100000 300 15 25344 6864 0 100000 300 16 3168 1056 0 100000 300 17 26400 11088 0 100000 300 18 16368 6864 0 100000 300 19 10560 13200 0 100000 300 20 16896 20592 0 100000 300 21 21648 24816 0 100000 300 22 17424 3168 0 100000 300 23 17424 23760 0 100000 300 24 3696 14784 0 100000 300 25 21648 13728 0 100000 300 26 6336 5280 0 100000 300 27 3696 16896 0 100000 300 28 1584 1584 0 100000 300 29 22704 24816 0 100000 300 30 10560 0 0 100000 300 31 1056 8976 0 100000 300 32 2640 22176 0 100000 300 33 3696 2112 0 100000 300 34 21120 15312 0 100000 300 -1 -
-
-
Hi,
Thank you for posting the data and making the issue reproducible. The problem lies in the way you define your indicator constraints. Indicator constraints have to hold a binary variable which either equals 1 or 0 on the left side of the \(\texttt{>>}\) operator, cf. addGenConstrIndicator. This means that your constraint
m.addConstr((f[i] >= Q/2) >> x[(i, j)] == 0)
results in the constraint \(0 \leq -1\), which is trivially infeasible. You can see this by writing out the LP file and analyzing it via
m.write("myLP.lp")This wrong construction does not generate an error. However, correct usage of parentheses provides the expected error
m.addConstr((f[i] >= Q/2) >> (x[(i, j)] == 0))
gurobipy.GurobiError: Indicator constraints can only be triggered by a single binary variable at a given valueThe correct way to write your indicator constraint would be
m.addConstr((x[(i, j)] == 0) >> (f[i] >= Q/2))
Best regards,
Jaromił1 -
-
Hi,
Thank you for giving really a good insight. Now I am able to understand how indicator constraint should be formulated through Gurobi. Although, I have a question / minor issue in my model.
When I ran the output, the output shows that the nodes explored are 0 nodes (81 simplex iterations) in 0.02 seconds. isn't it odd that it shows that it explored 0 nodes? I have run a couple of other examples and it shows some reasonable values in 100s or 1000s in terms of exploring nodes. Can you let me know why this issue is happening in my code? I will attach the snapshot of my output below. Thank you for your help.
0 -
-
-
Hi,
Your problem suffers from severe numerical issues, see Does my model have numerical issues?
You are working with a Big-M value of \(100^{100}\). Please first scale your model properly as described in our Guidelines for Numerical Issues.
Best regards,
Jaromił0 -
-
Hi,
Thank you for providing me these links. It really gave me an insight. And as per the documentation, my model has large matrix and objective coefficient ranges. It exceeds more than 10^9 (10^14 approx. in my case) . I am attaching the screenshot of that below. And apart from this, I reduced my Big-M value to 100^7 and finally I am able to see the nodes getting explored. But I am not able to get the objective value though I get best bound value after the time_limit has reached. I kept the time limit to a maximum of 5 minutes and still I am getting the same output. Can you let me know:
1) How to reduce the coefficient ranges?
2) How to obtain the objective value? I know one way is to reformulate my model. Though I suspect some issue in my time window constraints, I couldn't see any fault in my formulation in it. Or should I set the Params' time limit to a much higher value (say 10 minutes)? Please let me know. Thank you so far and thanks in advance for your help.
0 -
-
-
Hi,
\(10^{14}\) and \(100^7\) are still way too large.
I kept the time limit to a maximum of 5 minutes and still I am getting the same output.
5 minutes with the coefficients this bad is very likely not enough.
1) How to reduce the coefficient ranges?
You can find some scaling strategies in our documentation on Advanced user scaling. You should question yourself why do you require the coefficients to be that large and simultaneously so small. If for example a coefficient depicts some revenue in $, do you really care about the variables with \(1\)$ coefficients if there are variables with \(10^{14}\)$ coefficients? You could just drop all coefficients below \(10^6\) and adjust the orders of magnitude, i.e., one unit would equal \(10^6\)$ instead of \(1\)$. If you think that you cannot improve the ranges any further, you could try running the optimization with the parameter NumericFocus set to 3. However, setting parameters often cannot fully mitigate the numerical issues caused by bad scaling.
2) How to obtain the objective value? I know one way is to reformulate my model. Though I suspect some issue in my time window constraints, I couldn't see any fault in my formulation in it. Or should I set the Params' time limit to a much higher value (say 10 minutes)? Please let me know. Thank you so far and thanks in advance for your help.
First, you should reformulate your model and fix the coefficient and rhs ranges. This should already improve the solution finding process. Increasing the time limit is very likely required in this case. Additionally, you might want to look at the list of most important parameters.
Best regards,
Jaromił0 -
Post is closed for comments.
Comments
12 comments