modeling question
AnsweredI've solved a hierarchical facility location problem, but something is wrong.
this is a one of my answer.
I want the blue dots can be assigned to the green(Launch) as well. If the distance is shorter than passing through the recharge station(red).
I kept modify my model, it still doesn't work.
this my model
and the dis_func is a module I defined myself.
#dis_func
from random import randint
import math
def calculate_dis(location1,location2):
dis_bet={}
for i in location1:
for j in location2:
s1=math.sqrt(((location1[i])[0]-(location2[j])[0])**2+
((location1[i])[1]-(location2[j])[1])**2)
dis_bet.setdefault((i,j),s1)
return dis_bet
def calculate_dis_3(location1,location2,location3):
dis_bet_3={}
for i in location1:
for j in location2:
for k in location3:
s1=math.sqrt(((location1[i])[0]-(location2[j])[0])**2+
((location1[i])[1]-(location2[j])[1])**2)
s2=math.sqrt(((location2[j])[0]-(location3[k])[0])**2+
((location2[j])[1]-(location3[k])[1])**2)
dis_bet_3.setdefault((i,j,k),(s1+s2))
return dis_bet_3
def create_dots(length):
coor=[]
for i in range(length):
coor.append([])
for j in range(2):
s=randint(0, 100)
coor[i].append(s)
dots_dict={(i+1):coor[i] for i in range(length)}
return dots_dict
#gurobi model
from dis_func import calculate_dis,create_dots,calculate_dis_3
import gurobipy as gp
from gurobipy import *
from random import randint
import math
import matplotlib.pyplot as plt
import matplotlib.path as mpath
#create dots
demand=create_dots(10)
facility1=create_dots(20)
facility2=create_dots(20)
#calculate distance
#demand-facility1
dis_DF1=calculate_dis(facility1,demand)
#demand-facility2
dis_DF2=calculate_dis(facility2,demand)
#facility1-facility2
dis_FF=calculate_dis(facility2,facility1)
dis_FFD=calculate_dis_3(facility2,facility1,demand)
#cost of level 2 facility(Launch station)
CL=3000000
#cost of level 1 facility(recharge station)
CR=30000
#operation cost per unit
Co=100
#Maximum flying distance
E=40
m=Model("project_412")
#Uijk demadn k assigned to level 2 facility[i] and level 1 facility[j]
u=m.addVars(facility2,facility1,demand,vtype=GRB.BINARY,name='U')
'''
#the distance of demadn k assigned to level 2 facility[i] and level 1 facility[j]
d=m.addVars(dis_FFD,name="D")
'''
#Uik demand k directly assigned to level 1facility[i]
u1=m.addVars(facility2,demand,vtype=GRB.BINARY,name="U1")
'''
#the distance demand k directly assigned to level 1facility[i]
d1=m.addVars(dis_DF2,name="D1")
'''
#yi build a level 2 facility in [i]
y=m.addVars(facility2,vtype=GRB.BINARY,name="Y")
#zj build a level 1 facility in [j]
z=m.addVars(facility1,vtype=GRB.BINARY,name="Z")
#Objective
m.setObjective((quicksum(u1[i,k]*dis_DF2[i,k] for i in facility2 for k in demand)
+quicksum(u[i,j,k]*dis_FFD[i,j,k] for i in facility2 for j in facility1 for k in demand))*Co
+quicksum(y[i]*CL for i in facility2)
+quicksum(z[j]*CR for j in facility1)
,GRB.MINIMIZE)
#Constraint
#one demand only assigned once
for k in demand:
m.addConstr(quicksum(u[i,j,k]+u1[i,k]for i in facility2 for j in facility1)==1 ,name="only one" )
#demand is assigned to a pair of level 2 and level 1 facility or assigned to a level2 facility directly
for k in demand:
for i in facility2:
m.addConstr((quicksum(u[i,j,k] for j in facility1)+u1[i,k])==y[i] ,name="一定")
#only assigned to a facility which is build(level 1)
m.addConstrs((u[i,j,k]<=z[j] for i in facility2 for j in facility1 for k in demand),name="一級")
#only assigned to a facility which is build(level 2)
m.addConstrs((u[i,j,k]<=y[i] for i in facility2 for j in facility1 for k in demand),name="二級")
#only assigned to a facility which is build(level 2)
m.addConstrs((u1[i,k]<=y[i] for i in facility2 for k in demand),name="直接")
#non-negative constraint
m.addConstrs((u[i,j,k]>=0 for i in facility2 for j in facility1 for k in demand),name="非負")
m.addConstrs((u1[i,k]>=0 for i in facility2 for k in demand),name="非負2")
#flying range constraint(level 2- level 1)
for i in facility2:
for j in facility1:
m.addConstr((y[i]*z[j]*dis_FF[i,j])<=E,name="飛距限制1")
#flying range constraint(level 1-Demand)
for i in facility2:
for j in facility1:
for k in demand:
m.addConstr((z[j]*u[i,j,k]*dis_DF1[j,k])<=E ,name="飛距限制2")
#flying range constraint(level 2-Demand)
for i in facility2:
for k in demand:
m.addConstr((y[i]*u1[i,k]*dis_DF2[i,k])<=E ,name="飛距限制3")
#optimize
m.optimize()
#print answer
for v in m.getVars():
if v.X==1:
print('%s %g' % (v.varName, v.X))
#matplotlib
for i in facility2:
if y[i].X==1:
plt.plot((facility2[i])[0],(facility2[i])[1],"g^")
for k in demand:
plt.plot((demand[k])[0],(demand[k])[1],'bo')
for j in facility1:
if z[j].X==1:
plt.plot((facility1[j])[0],(facility1[j])[1],"ro")
for i in facility2:
for j in facility1:
for k in demand:
if u[i,j,k].X==1:
x=[(facility2[i])[0],(facility1[j])[0]]
y=[(facility2[i])[1],(facility1[j])[1]]
lines1=plt.plot(x,y,"k")
f=[(facility1[j])[0],(demand[k])[0]]
l=[(facility1[j])[1],(demand[k])[1]]
lines2=plt.plot(f,l,"k--")
#infeasible
#m.computeIIS()
#print(m.getAttr(GRB.Attr.IISConstr))
I don't know why my u1 variable are always none, my objective should automatically assigned demand to level 2 facility(green) if the distance (level 2-demand) is shorter than (level 2- level 1- demand), but every result I get, no matter how many demands and facility candidate points, are like the picture I post beginning of the post.
I've tried to add
#nested
for i in facility2:
for k in demand:
if dis_DF2[i,k]<E:
m.addConstr(u1[i,k]==y[i],name="強制")
and the model became infeasible
I've tried to restrict at least one u1 variable equals 1 and the model became infeasible again.
Is there a mistake in my model?
-
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?. -
-
-
It looks like the "only one" constraints are the problem:
for k in demand:
m.addConstr(quicksum(u[i,j,k] + u1[i,k] for i in facility2 for j in facility1) == 1, name="only one")The mathematical formulation of these constraints is
$$\begin{align*}\sum_{i \in F_2} \sum_{j \in F_1} (u_{ijk} + u^1_{ik}) = 1 \quad \forall k \in D.\end{align*}$$
The \( u^1 \) variables are enclosed in a summation over \( F_1 \). As a result, each \( u^1_{ik} \) variable shows up \( |F_1| = 20 \) times on the left-hand side of the constraint:
$$\begin{align*}\sum_{i \in F_2} \sum_{j \in F_1} u_{ijk} + |F_1| \sum_{i \in F_2} u^1_{ik} = 1 \quad \forall k \in D.\end{align*}$$
As a result, setting any of the binary \( u^1 \) variables to \( 1 \) results in an infeasible model. To fix this, move the \( u^1_{ik} \) terms outside of the summation over \( F_1 \):
$$\begin{align*}\sum_{i \in F_2} \bigg( u^1_{ik} + \sum_{j \in F_1} u_{ijk} \bigg) = \sum_{i \in F_2} \sum_{j \in F_1} u_{ijk} + \sum_{i \in F_2} u^1_{ik} = 1 \quad \forall k \in D.\end{align*}$$
You can write these constraints concisely using tupledict.sum():
for k in demand:
m.addConstr(u.sum('*', '*', k) + u1.sum('*', k) == 1, name='only_one')A few other comments:
- Instead of checking the value of a binary variable with \( \texttt{if y[i].X == 1} \), use something like \(\texttt{if y[i].X > 1 - 1e-5} \). Solvers use numerical tolerances to solve problems, so you may encounter values for integer variables that are slightly non-integral.
- The non-negativity constraints aren't necessary. By default, variables added to a model with Model.addVars() or Model.addVar() have a lower bound of \( 0 \).
0 -
-
Thanks for your help, I actually try the same thing before, using different form
for k in demand:
m.addConstr(quicksum(u[i,j,k] for i in facility2 for j in facility1)+quicksum(u1[i,k] for i in facility2)==1, name="only one")but somehow it didn't work , so I thought this kind of constraint is not correct, and after I try your answer, I decide to try my origin one again, because I think the concept is the same as your answer, and it works. I have no idea why it wasn't working then, but anyway thanks for your help!
0 -
Post is closed for comments.
Comments
3 comments