Skip to main content

Facility Layout Problem, non-overlapping rectangles, model appears to be infeasible

Answered
Carole De Sloovere
Edited
  • Gurobi-versary
  • First Comment
  • First Question

Dear,

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)

Related to

0

Comments

5 comments

Date Votes
  • Jaromił Najman
    • 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
  • Carole De Sloovere
    Edited
    • Gurobi-versary
    • First Comment
    • First Question

    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
  • Jaromił Najman
    • Gurobi Staff
    Edited

    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]
    End

    Solving 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.25

    In 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
  • Carole De Sloovere
    • Gurobi-versary
    • First Comment
    • First Question

    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
  • Jaromił Najman
    • Gurobi Staff

    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.

Didn't find what you were looking for?

New post