Facility Layout Problem, non-overlapping rectangles, model appears to be infeasible
AnsweredDear,
I try to model a Facility layout problem, where:
- data is a matrix with the flow of materials
- data_in: flow from o1 to machines
- data_out: flow from machines to o2
- W: width of floor
- L: length
- w: widths machines
- l: lengths machines
- o1: obstacle in left bottom, coordinates of upper right corner
- o2: obstacle right on top, coordinates of lower left corner
I want to minimize the flow of materials multiplied by the distance that they travel. The most important constraints are the ones that make sure that there is no overlap, however, these are the ones that make my model infeasible. I know there are some solutions possible.
Any tips or help is appreciated!
Thank you in advance.
Kind regards, Carole
w = dict({'BF1':1,'BF2':1,'BF3':1,'BF4':1,'SM':1,'WT':1,'BL':1})
l = dict({'BF1':1,'BF2':1,'BF3':1,'BF4':1,'SM':2.5,'WT':2.5,'BL':3})
W = 9.15
L = 8
o1 = (1.5,3)
o2 = (6.15,6)
data = [[0, 0, 0, 0, 0, 486, 0, 486], [0, 0, 0, 0, 0, 882, 0, 882], [0, 0, 0, 0, 0, 495, 0, 495], [0, 0, 0, 0, 0, 495, 0, 495], [486, 882, 495, 495, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 1476, 0], [0, 0, 0, 0, 1476, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0]]
data_in = [486, 882, 495, 495, 0, 0, 0]
data_out = [486, 882, 495, 495, 0, 0, 0]
def solve_FLP(w,l,L,W,o1,o2,data, data_in, data_out):
m = gp.Model('FLP')
T = len(machines)
x = m.addVars(range(T),lb = 0, ub = W, vtype = GRB.CONTINUOUS, name='x') #x and y coordinates of centroids of machines
y = m.addVars(range(T),lb = 0, ub = L, vtype = GRB.CONTINUOUS, name='y')
da = m.addVars(range(T), vtype = GRB.BINARY, name='da')
db = m.addVars(range(T), vtype = GRB.BINARY, name='db')
delta = m.addVars(range(T), range(T), range(4), vtype = GRB.BINARY, name='delta')
deltao1 = m.addVars(range(T), range(2), vtype = GRB.BINARY, name='deltao1')
deltao2 = m.addVars(range(T), range(2), vtype = GRB.BINARY, name='deltao2')
x_dist0 = m.addVars(range(T),range(T),name = "x_dist0")
y_dist0 = m.addVars(range(T),range(T),name = "y_dist0")
dist0 = m.addVars(range(T),range(T),name = "dist0")
pow_x0 = m.addVars(range(T),range(T),name = "pow_x0")
pow_y0 = m.addVars(range(T),range(T),name = "pow_y0")
x_dist1 = m.addVars(range(T),name = "x_dist1")
y_dist1 = m.addVars(range(T),name = "y_dist1")
dist1 = m.addVars(range(T),name = "dist1")
pow_x1 = m.addVars(range(T),name = "pow_x1")
pow_y1 = m.addVars(range(T),name = "pow_y1")
x_dist2 = m.addVars(range(T),name = "x_dist2")
y_dist2 = m.addVars(range(T),name = "y_dist2")
dist2 = m.addVars(range(T),name = "dist2")
pow_x2 = m.addVars(range(T),name = "pow_x2")
pow_y2 = m.addVars(range(T),name = "pow_y2")
a = list(w.values())
b = list(l.values())
M = 100
# Orientation
for i in range(T):
m.addConstr(da[i] + db[i] == 1)
# no overlap with borders
for i in range(T):
m.addConstr(x[i] >= (da[i]*a[i] + db[i]*b[i])/2)
m.addConstr(x[i] <= W - (da[i]*a[i] + db[i]*b[i])/2)
m.addConstr(y[i] >= ((1-da[i])*a[i] + (1-db[i])*b[i])/2)
m.addConstr(y[i] <= L - ((1-da[i])*a[i] + (1-db[i])*b[i])/2)
#no overlap of machines
for i in range(T):
for j in range(T):
if i < j:
m.addConstr(x[i] + (da[i]*a[i] + db[i]*b[i])/2 <= x[j]-(da[j]*a[j] + db[j]*b[j])/2 + M*delta[(i,j,0)])
m.addConstr(x[j] + (da[j]*a[j] + db[j]*b[j])/2 <= x[i]-(da[i]*a[i] + db[i]*b[i])/2 + M*delta[(i,j,1)])
m.addConstr(y[i] + ((1-da[i])*a[i] + (1-db[i])*b[i])/2 <= y[j]-((1-da[j])*a[j] + (1-db[j])*b[j])/2 + M*delta[(i,j,2)])
m.addConstr(y[j] + ((1-da[j])*a[j] + (1-db[j])*b[j])/2 <= y[i]-((1-da[i])*a[i] + (1-db[i])*b[i])/2 + M*delta[(i,j,3)])
m.addConstr(gp.quicksum(delta[(i,j,k)] for k in range(4)) <= 3)
else:
m.addConstrs(delta[i,j,k] == 0 for k in range(4))
m.addConstr(gp.quicksum(delta[(i,j,k)] for k in range(4)) <= 3)
#no overlap machine <> obstacle
for i in range(T): #obstacle 1
m.addConstr(x[i] - (da[i]*a[i] + db[i]*b[i])/2 >= o1[0]+ M*(1-deltao1[(i,1)]))
m.addConstr(o1[1] <= y[i] - ((1-da[i])*a[i] + (1-db[i])*b[i])/2 + M*(1-deltao1[(i,1)]))
m.addConstr(gp.quicksum(deltao1[(i,k)] for k in range(2)) <= 1)
for i in range(T): #obstacle 2
m.addConstr(x[i] + (da[i]*a[i] + db[i]*b[i])/2 <= o2[0] + M*deltao2[(i,0)])
m.addConstr(y[i] + ((1-da[i])*a[i] + (1-db[i])*b[i])/2 <= o2[1] + M*deltao2[(i,1)])
m.addConstr(gp.quicksum(deltao2[(i,k)] for k in range(2)) <= 1)
# Rectilinear distance modeling
for i in range(T):
for j in range(T):
m.addConstr(pow_x0[(i,j)] == x[j] - x[i])
m.addGenConstrAbs(x_dist0[(i,j)],pow_x0[(i,j)] ,"absconstr")
m.addConstr(pow_y0[(i,j)] == y[j] - y[i])
m.addGenConstrAbs(y_dist0[(i,j)],pow_y0[(i,j)],"absconstr")
m.addConstr(dist0[(i,j)] == x_dist0[(i,j)] + y_dist0[(i,j)])
for i in range(T):
m.addConstr(pow_x1[i] == o1[0] - x[i])
m.addGenConstrAbs(x_dist1[i],pow_x1[i] ,"absconstr")
m.addConstr(pow_y1[i] == o1[1] - y[i])
m.addGenConstrAbs(y_dist1[i],pow_y1[i],"absconstr")
m.addConstr(dist1[i] == x_dist1[i] + y_dist1[i])
for i in range(T):
m.addConstr(pow_x2[i] == o2[0] - x[i])
m.addGenConstrAbs(x_dist2[i],pow_x2[i] ,"absconstr")
m.addConstr(pow_y2[i] == o2[1] - y[i])
m.addGenConstrAbs(y_dist2[i],pow_y2[i],"absconstr")
m.addConstr(dist2[i] == x_dist2[i] + y_dist2[i])
m.update()
obj = gp.quicksum(gp.quicksum(dist0[(i,j)]*data[i][j] for i in range(T)) for j in range(T)) + gp.quicksum(dist1[i]*data_in[i] for i in range(T)) + gp.quicksum(dist2[i]*data_out[i] for i in range(T))
m.setObjective(obj, GRB.MINIMIZE)
m.optimize()
print(m.status)
if m.status == GRB.OPTIMAL:
print('Solution:')
print('cost = ', m.objVal)
xvals = { f'x{k}' : v.X for k,v in x.items() }
print(xvals)
yvals = { f'y{k}' : v.X for k,v in y.items() }
print(yvals)
distvals = { f'dist21{k}' : v.X for k,v in dist2.items() }
print(distvals)
print(a)
print(b)
da_values = { f'da{k}' : v.X for k,v in da.items() }
print(da_values)
db_values = { f'db{k}' : v.X for k,v in db.items() }
print(db_values)
delta_o1 = { f'delta01{k}' : v.X for k,v in deltao1.items() }
print(delta_o1)
print(o1)
print(o2)
else:
print('No solution found.')
solve_FLP(w,l,L,W,o1,o2,data, data_in, data_out)-
-
Gurobi Staff
Hi Carole,
Please note that the code you provided is not a minimal reproducible example, because all input data is missing.
Regarding the infeasibility, you might want to have a look at the Knowledge Base article How do I determine why my model is infeasible?
Best regards,
Jaromił0 -
-
Dear,
Thank you for the link, however I did not manage to find the problem yet.
I updated my previous post and added the data.
Thank you in advance.
Kind regards,
Carole
0 -
-
-
HI Carole,
Thank you updating the code. Unfortunately, there is still the definition of the \(\texttt{machines}\) object missing.
I tried defining it as
machines = [0,1,2,3,4,5,6]
Would this reproduce the issue?
With the above machine setting, you can compute an IIS via
[...]
else:
m.computeIIS()
m.write("iis.ilp")
print('No solution found.')and analyze the \(\texttt{iis.ilp}\) file.
The IIS compute by Gurobi in this case is given as
\ Model FLP_copy
\ LP format - for model browsing. Use MPS format to capture full model detail.
Minimize
Subject To
R292: x[4] - 0.5 da[4] - 1.25 db[4] + 100 deltao1[4,1] >= 101.5
R293: y[4] + 0.5 da[4] + 1.25 db[4] - 100 deltao1[4,1] >= -95.25
R481: x[4] + pow_x1[4] = 1.5
R482: y[4] + pow_y1[4] = 3
Bounds
x[4] free
y[4] free
Binaries
da[4] db[4] deltao1[4,1]
EndSolving equality constraints \(\texttt{R481,R482}\) for \(\texttt{x[4],y[4]}\) and substituting in \(\texttt{R292,R293}\) gives
- pow_x1[4] - 0.5 da[4] - 1.25 db[4] + 100 deltao1[4,1] >= 100
- pow_y1[4] + 0.5 da[4] + 1.25 db[4] - 100 deltao1[4,1] >= -98.25In the above, \(\texttt{deltao1[4,1]}\) is forced to \(1\) and \(\texttt{pow_x1[4]=da[4]=db[4]=}\)\(0\), because \(\texttt{pow_x1[4]}\) can take only non-negative values.
This makes inequality \(\texttt{R293}\) read
- pow_y1[4] >= 1.75
which is infeasible because \(\texttt{pow_y1[4]}\) is non-negative.
Best regards,
Jaromił0 -
-
Dear Jaromil,
Thank you for investigating my problem, I understand what goes wrong but I do not know how to fix it.
My intention with all the pow_ variables was to take the difference of the coordinates and I wanted to use the addGenConstrAbs afterwards in order to take the absolute value of the pow_ variables. This was to calculate the rectilinear distances. However, if they can only take positive values, this will indeed make it infeasible.
Is there any way to solve this?
Kind regards,
Carole
0 -
-
-
Hi Carole,
My intention with all the pow_ variables was to take the difference of the coordinates
In this case the \(\texttt{pow_}\) variables can be negative. Assume you have 2 coordinates \(x_1=2,x_2=5\). If you compute \(pow_{x_1,x_2} = x_1 - x_2\) then \(pow_{x_1,x_2}\) is negative. You can then take the abs value of the \(pow_{x_1,x_2}\) variable.
Best regards,
Jaromił0 -
Please sign in to leave a comment.
Comments
5 comments