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
%%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:
- 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.
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
- 2d Bin Packing, https://yetanothermathprogrammingconsultant.blogspot.com/2021/02/2d-bin-packing.html
- 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
Hello Erwin.
ReplyDeleteThank 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.
This is excellent! Thanks.
Delete