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):

# used bin
for i in range(n):
for k in range(m):

# 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):

# 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.

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(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)]

At the end we define the objective:

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

for i in range(n):
for j in range(m):