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.

References

1. 2d Bin Packing, https://yetanothermathprogrammingconsultant.blogspot.com/2021/02/2d-bin-packing.html
2. 2d bin packing using or-tools: AddNoOverlap2D and OnlyEnforceIf gives MODEL_INVALID, https://stackoverflow.com/questions/66286586/2d-bin-packing-using-or-tools-addnooverlap2d-and-onlyenforceif-gives-model-inva