Solution not Feasible in Orienteering Problem in Python Gurobi
回答済みI am coding Orienteering problem using Gurobi in Python. Somehow I managed to eliminate syntax error and coded pure logic based on this famous equations for Orienteering problem as follows:
I am however getting solution is infeasible. Does that mean Gurobi cant solve this? Or somehow there is something else I am missing. I have written the code indicating which constraint matches which which equation so allow easy trouble shooting. Any assistance shall be appreciated.
import random
import gurobipy as grb
import math
n = 4
Distance = 50000000
def distance(points, i, j):
dx = points[i][0] - points[j][0]
dy = points[i][1] - points[j][1]
return math.sqrt(dx*dx + dy*dy)
random.seed(1)
points = []
for i in range(n):
points.append((random.randint(0,100),random.randint(0,100)))
opt_model = grb.Model(name="MILP Model")
# <= Variables
x_vars = {}
for i in range(n):
for j in range(n):
x_vars[i,j] = opt_model.addVar(vtype=grb.GRB.BINARY,
name='e'+str(i)+'_'+str(j))
u={}
for i in range(1,n):
u[i]=opt_model.addVar(vtype=grb.GRB.INTEGER,
name='e'+str(i))
# <= Constraint (Mandatory Nodes) Eq(1)
opt_model.addConstr((grb.quicksum(x_vars[1,j] for j in range(1,n))) == 1)
opt_model.addConstr((grb.quicksum(x_vars[i,n-1] for i in range(n-1))) == 1)
# <= Constraint (Distance) Eq(3)
for i in range(n-1):
opt_model.addConstr(grb.quicksum(x_vars[i,j]*distance(points, i, j) for j in range(1,n)) <= Distance)
x_vars[i,i].ub = 0
# <= Constraint (Equality & Coming in to node and going out should be 1 each) Eq(2)
for k in range(1, n-1):
opt_model.addConstr(grb.quicksum(x_vars[i,k] for i in range(n-1))
== grb.quicksum(x_vars[k,j] for j in range(1, n)) <=1)
# <= Constraint (Subtour elimination) Eq(4) Eq(5)
for i in range(1,n):
opt_model.addConstr(u[i] <= n)
for i in range(1,n):
opt_model.addConstr(u[i] >= 2)
for i in range(1,n):
for j in range(1,n):
opt_model.addConstr((u[i] - u[j] +1 == (n-1)*(1-x_vars[j,i])))
opt_model.update()
# <= objective is to maximize Eq(1)
objective = grb.quicksum(x_vars[i,j]
for i in range(1, n-1)
for j in range(1, n))
opt_model.ModelSense = grb.GRB.MAXIMIZE
opt_model.setObjective(objective)
opt_model.update()
opt_model.optimize()
solution = opt_model.getAttr('x', x_vars )
print solution-
正式なコメント
-
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 Abdul,
Whenever you get an infeasible model, it is useful to run IIS, to understand why your model is infeasible. Here is a link an article on Managing Infeasibility
0 -
-
Dear Jennifer,
Thanks for the feedback. I made the code little more compact and scannable and ran IIS but it says it cannot compute IIS on a feasible model
import random
import gurobipy as grb
import math
n = 4
Distance = 50000000
def distance(points, i, j):
dx = points[i][0] - points[j][0]
dy = points[i][1] - points[j][1]
return math.sqrt(dx*dx + dy*dy)
random.seed(1)
points = []
for i in range(n):
points.append((random.randint(0,100),random.randint(0,100)))
opt_model = grb.Model(name="MILP Model")
# <= Variables
x_vars = {}
for i in range(n):
for j in range(n):
x_vars[i,j] = opt_model.addVar(vtype=grb.GRB.BINARY,
name='e'+str(i)+'_'+str(j))
u={}
for i in range(1,n):
u[i]=opt_model.addVar(vtype=grb.GRB.INTEGER,
name='e'+str(i))
# <= Constraint (Mandatory Edges and excluding vertexes) Eq(1)
opt_model.addConstr((grb.quicksum(x_vars[1,j] for j in range(1,n))) == 1)
opt_model.addConstr((grb.quicksum(x_vars[i,n-1] for i in range(n-1))) == 1)
opt_model.addConstr((grb.quicksum(x_vars[i,i] for i in range(n-1))) == 0)
# <= Constraint (Distance) Eq(3)
for i in range(n-1):
opt_model.addConstr(grb.quicksum(x_vars[i,j]*distance(points, i, j) for j in range(1,n)) <= Distance)
# <= Constraint (Equality & Single edge in and out) Eq(2)
for k in range(1, n-1):
opt_model.addConstr(grb.quicksum(x_vars[i,k] for i in range(n-1))
== grb.quicksum(x_vars[k,j] for j in range(1, n)) <=1)
# <= Constraint (Subtour elimination) Eq(4) Eq(5)
for i in range(1,n):
opt_model.addConstr(2 <= u[i] <= n)
for i in range(1,n):
for j in range(1,n):
opt_model.addConstr((u[i] - u[j] +1 <= (n-1)*(1-x_vars[j,i])))
# <= objective (maximize) Eq(1)
objective = grb.quicksum(x_vars[i,j]
for i in range(1, n-1)
for j in range(1, n))
opt_model.ModelSense = grb.GRB.MAXIMIZE
opt_model.setObjective(objective)
opt_model.update()
opt_model.computeIIS()
opt_model.write("model.ilp")
solution = opt_model.getAttr('x', x_vars )
print solution0 -
-
-
Try running computeIIS() after optimize(). Check the status to see if the model is still infeasible, if it is, then run IIS.
m.optimize()
status = m.status
if status == GRB.Status.INFEASIBLE: m.computeIIS()
if m.IISMinimal:
print('IIS is minimal\n')
else:
print('IIS is not minimal\n')
print('\nThe following constraint(s) cannot be satisfied:')
for c in m.getConstrs():
if c.IISConstr:
print('%s' % c.constrName)1 -
-
-
The two version of your code produce different models. The first one is indeed infeasible since you are setting
x_vars[i,i].ub = 0
which fixes \(x\_vars[i,i]\) to 0, but on the other hand you also haveu[i] - u[j] +1 == (n-1)*(1-x_vars[j,i])
which for \(i=j\) becomes
(n-1) * x_vars[i,i] == n-2
and for \(n \neq 2\) this contradicts \(x\_vars[i,i] = 0\).
The second model is feasible and hence you cannot compute an IIS.
1 -
-
Well thanks for checking this out in detail. Yes now the solution is feasible. I have one last question though. How is it possible to select only the edges that are 1 i.e. part of the solution. I tried the following code but it doesnt select anything and prints all the edges.
select = grb.tuplelist((i,j) for i,j in x_vars.keys() if x_vars.values() > 0.5)
But it takes all the edges instead of the one with values greater then 0.5
0 -
-
-
If there is a solution (e.g. after optimize() has finished), you can access the variables' values in the current solution using the X attribute. E.g.:
select = grb.tuplelist((i,j) for i,j in x_vars.keys() if x_vars[i,j].X > 0.5)
1 -
-
Solved. Thanks a lot Dr. Silke!
0 -
投稿コメントは受け付けていません。
コメント
8件のコメント