Saturday, February 20, 2021

2d bin packing with Google OR-Tools CP-SAT, how to fix?

In [1] I looked at some MIP formulations for a 2d bin packing problem. OR-Tools has a nice model.AddNoOverlap2D constraint function. So, I was excited to try this out. (In these times I get quickly enthusiastic about these things). 

The idea is as follows. ortools has this construct Constraint.OnlyEnforce If(boolean). This is essentially an indicator constraint: \((\color{darkred}x\ge 5)\text{.OnlyEnforceIf($\color{darkred}\delta$)}\) means \[\color{darkred}\delta=1 \Rightarrow \color{darkred}x \ge 5\] This inspires me to develop a model along the lines of: \[\begin{align} \min& \sum_k \color{darkred}u_k\\  & \color{darkred}\beta_{i,j} = 0 \Rightarrow \color{darkred}b_i\ne\color{darkred}b_j \\ & \color{darkred}\beta_{i,j} = 1 \Rightarrow \mathbf{NoOverlap2D}(i,j) \\ & \color{darkred}u_k = 0 \Rightarrow \color{darkred}b_i \lt k \end{align}\] Here \(\color{darkred}b_i\) is the bin number where item \(i\) is placed and \(\color{darkred}\beta_{i,j}\) is true if  \(\color{darkred}b_i=\color{darkred}b_j\).


OK, so I started coding:


from ortools.sat.python import cp_model

#---------------------------------------------------
# data 
#---------------------------------------------------

data = {'bin':{'h':60,'w':40},
        'cat1':{'w': 7,'h':12,'items':10},
        'cat2':{'w': 9,'h': 3,'items':10},
        'cat3':{'w': 5,'h':14,'items':10},
        'cat4':{'w':13,'h': 9,'items':10},
        'cat5':{'w': 6,'h': 8,'items': 5},
        'cat6':{'w':20,'h': 5,'items': 5}}        

#
# extract data for easier access
#

# bin width and height
H = data['bin']['h']
W = data['bin']['w']

# h,w,cat for each item
h = [data[cat]['h'] for cat in data if cat!='bin' for i in range(data[cat]['items'])]
w = [data[cat]['w'] for cat in data if cat!='bin' for i in range(data[cat]['items'])]
cat = [cat for cat in data if cat!='bin' for i in range(data[cat]['items'])]
n = len(h)  # number of items
m = 10      # number of bins

#---------------------------------------------------
# or-tools model 
#---------------------------------------------------


model = cp_model.CpModel()

#
# variables
#

# x1,x2 and y1,y2 are start and end
x1 = [model.NewIntVar(0,W-w[i],f'x1.{i}') for i in range(n)]
x2 = [model.NewIntVar(w[i],W,f'x2.{i}') for i in range(n)]

y1 = [model.NewIntVar(0,H-h[i],f'y1.{i}') for i in range(n)]
y2 = [model.NewIntVar(h[i],H,f'y2.{i}') for i in range(n)]

# interval variables
xival = [model.NewIntervalVar(x1[i],w[i],x2[i],f'xival{i}') for i in range(n)]
yival = [model.NewIntervalVar(y1[i],h[i],y2[i],f'yival{i}') for i in range(n)]

# bin numbers
b = [model.NewIntVar(0,m,f'b{i}') for i in range(n)]

# b2[(i,j)] = true if b[i]=b[j] for i<j
b2 = {(i,j):model.NewBoolVar(f'b2.{i}.{j}') for j in range(n) for i in range(j)}

# used bins
u = [model.NewBoolVar(f'u{k}') for k in range(m)]


#
# constraints
#

# no overlap for items in same bin
for j in range(n):
  for i in range(j):
    model.Add(b[i] != b[j]).OnlyEnforceIf(b2[(i,j)].Not())
    model.AddNoOverlap2D([xival[i],xival[j]],[yival[i],yival[j]]).OnlyEnforceIf(b2[(i,j)])

# used bin
for i in range(n):
  for k in range(m):
    model.Add(b[i]<k).OnlyEnforceIf(u[k].Not())

# objective
model.Minimize(sum([u[k] for k in range(m)]))


#
# solve model
#
solver = cp_model.CpSolver()
rc = solver.Solve(model)
print(rc)
print(solver.StatusName())


Work in progress, so I expect to see problems. When I run this I see:

OK. That is not helpful. After a bit of trial and error, it looks like OnlyEnforceIf is not implemented for the AddNoOverlap2D constraints. Indeed, we can verify this with solver.parameters.log_search_progress = True [2].  At first, I did not see anything (I was running this in a Jupyter notebook on colab.google.com; Jupyter notebooks are notorious for not displaying logs and error messages). Running it locally from the command line, I saw:

Enforcement literal not supported in constraint: enforcement_literal: 256 no_overlap_2d { x_intervals: 0 x_intervals: 1 y_intervals: 50 y_intervals: 51 }

This confirms this is not supported.

I am not sure how to fix this with a simple work-around. If desperate, I could try to use some of the MIP formulations we did earlier.

Well, a bit too early to give up. Here is an alternative approach based on the proposal by Aliaksei Vaidzelevich in the comments below.


Attempt number 2


Instead of using the NoOverlap2D constraints within each bin, we just create one large area. It's width is \(\color{darkblue}m\cdot \color{darkblue}W\) and the height is just \(\color{darkblue}H\). Here \(\color{darkblue}W\) and \(\color{darkblue}H\) are the width and height of each bin, and \(\color{darkblue}m\) is the number if bins we consider. The high-level model can look like:\[\begin{align} \min\>& \color{darkred} z  \\ & \color{darkred}z = \max_i{\color{darkred}b_i} \\  & \color{darkred}x'_i = \color{darkred}x_i + \color{darkred}b_i \cdot \color{darkblue}W \\ & \mathbf{NoOverlap2D}(\color{darkred}x',\color{darkred}y)  \\ &  0 \le \color{darkred}x_i \le \color{darkblue}W-\color{darkblue}w_i \\ & 0 \le \color{darkred}y_i \le \color{darkblue}H-\color{darkblue}h_i \end{align}\] 


This solves almost immediately. The code looks like:


%%time
from ortools.sat.python import cp_model
import pandas

#---------------------------------------------------
# data 
#---------------------------------------------------

data = {'bin':{'h':60,'w':40},
        'cat1':{'w': 7,'h':12,'items':10},
        'cat2':{'w': 9,'h': 3,'items':10},
        'cat3':{'w': 5,'h':14,'items':10},
        'cat4':{'w':13,'h': 9,'items':10},
        'cat5':{'w': 6,'h': 8,'items': 5},
        'cat6':{'w':20,'h': 5,'items': 5}}      

#
# extract data for easier access
#

# bin width and height
H = data['bin']['h']
W = data['bin']['w']

# h,w,cat for each item
h = [data[cat]['h'] for cat in data if cat!='bin' for i in range(data[cat]['items'])]
w = [data[cat]['w'] for cat in data if cat!='bin' for i in range(data[cat]['items'])]
cat = [cat for cat in data if cat!='bin' for i in range(data[cat]['items'])]
n = len(h)  # number of items
m = 10  # max number of bins

#---------------------------------------------------
# or-tools model 
#---------------------------------------------------

model = cp_model.CpModel()

#
# variables
#

# x and y
x = [model.NewIntVar(0,W-w[i],f'x{i}') for i in range(n)]

xb1 = [model.NewIntVar(0,m*W-w[i],f'xb1.{i}') for i in range(n)]
xb2 = [model.NewIntVar(w[i],m*W,f'xb2.{i}') for i in range(n)]

y1 = [model.NewIntVar(0,H-h[i],f'y1.{i}') for i in range(n)]
y2 = [model.NewIntVar(h[i],H,f'y2.{i}') for i in range(n)]

# interval variables
xival = [model.NewIntervalVar(xb1[i],w[i],xb2[i],f'xival{i}') for i in range(n)]
yival = [model.NewIntervalVar(y1[i],h[i],y2[i],f'yival{i}') for i in range(n)]

# bin numbers
b = [model.NewIntVar(0,m-1,f'b{i}') for i in range(n)]

# objective
z = model.NewIntVar(0,m-1,'z')

#
# constraints
#
for i in range(n):
  model.Add(xb1[i] == x[i] + b[i]*W)
  model.Add(xb2[i] == xb1[i] + w[i])

model.AddNoOverlap2D(xival,yival)

model.AddMaxEquality(z,[b[i] for i in range(n)])

# objective
model.Minimize(z)    

#
# solve model
#
solver = cp_model.CpSolver()
# log does not work inside a Jupyter notebook
solver.parameters.log_search_progress = True
rc = solver.Solve(model)
print(f"return code:{rc}")
print(f"status:{solver.StatusName()}")

if rc == 4:
    df = pandas.DataFrame({ 
        'bin' : [solver.Value(b[i]) for i in range(n)],
        'x'   : [solver.Value(x[i]) for i in range(n)],
        'y'   : [solver.Value(y1[i]) for i in range(n)],
        'w'   : w,
        'h'   : h})
    display(df)


The results look like:



and the total turnaround time is:




Here is a picture of the results. The picture is wide as I started out with a maximum of 10 bins. In the end, we only need two tightly packed bins.




Notes:
  • The interesting thing in this model is that I don't even really use \(\color{darkblue}m=10\) (the maximum number of bins). It is just an upper bound on \(\color{darkred}z\) and \(\color{darkred}b_i\). We can make it large, say \(\color{darkblue}m=1000\)  without any problems.
  • The performance is really bad when I replaced the one big NoOverlap2D constraint with a bunch of pairwise ones. Mathematically this is the same, but for the solver it is not. That may indicate the original formulation I tried would not be very good.
  • Ortools only works with integer variables and integer data. You may need to scale things so that you can work with integers.
  • For some data sets, things may take a longer time. It is a good idea to add a time limit, so the solver can return a good solution. It can be optimal but just not proven optimal.  

Here is the result for a larger data set:


This model was declared "feasible" after hitting a time limit of 10 minutes. I don't think there is any way we can improve upon this by a whole bin. Still, this is very good performance compared to MIP formulations solved with a MIP solver.

2 comments:

  1. Hello Erwin.

    Thank you for your blog. I always read your posts with interest. It may be natural to use optional intervals, but it seems AddNoOverlap2D doesn’t support them. So we must use regular intervals.

    x1 = [model.NewIntVarFromDomain(cp_model.Domain.FromIntervals(
    [[j * W, (j + 1) * W - w[i]] for j in range(m)]), f'x1[{i}]') for i in range(n)]
    x2 = [model.NewIntVarFromDomain(cp_model.Domain.FromIntervals(
    [[j * W + w[i], (j + 1) * W] for j in range(m)]), f'x2[{i}]') for i in range(n)]
    y1 = [model.NewIntVar(0, H - h[i], f'y1[{i}]') for i in range(n)]
    y2 = [model.NewIntVar(h[i], H, f'y2[{i}]') for i in range(n)]

    Now let us define literals such that lit[i][j] = True iff the ith interval belongs to the jth bin:

    lit = [[model.NewBoolVar(f'lit[{i}][{j}]') for j in range(m)] for i in range(n)]

    Now we assert that every rectangle belongs to exactly one bin. Also we must refine left and right boundaries of every rectangle:

    for i in range(n):
    model.Add(sum(lit[i]) == 1)
    model.Add(sum(j * problem.w_bin * lit[i][j] for j in range(m)) <= x1[i])
    model.Add(sum((j + 1) * problem.w_bin * lit[i][j] for j in range(m)) >= x2[i])

    Then let us define regular intervals and use AddNoOverlap2D function:

    x_interval = [model.NewIntervalVar(x1[i], w[i], x2[i], f'x_interval[{i}]')
    for i in range(n)]
    y_interval = [model.NewIntervalVar(y1[i], h[i], y2[i], f'y_interval[{i}]')
    for i in range(n)]
    model.AddNoOverlap2D(x_interval, y_interval)

    At the end we define the objective:

    obj = model.NewIntVar(0, m, 'obj')

    for i in range(n):
    for j in range(m):
    model.Add(obj >= j + 1).OnlyEnforceIf(lit[i][j])

    model.Minimize(obj)

    Now we are ready to solve the problem:

    solver = cp_model.CpSolver()
    solver.parameters.num_search_workers = 8
    status = solver.Solve(model)

    It takes about 0.5 sec to solve the problem on my laptop.

    ReplyDelete