Packing rectangles in a fixed size rectangle

In this question, we consider how to pack a collection of polyominos into a rectangle of a fixed size.

In my case, I'm concerned with the special case where all polyominos are rectangles and while the problem is clearly still hard (it has subset sum as a special case), I am wondering if this additional structure is enough to come up with more efficient solutions than the excellent SAT/MILP-based ones that are given in the original questions.

A lot of questions already exist on rectangle packing, but most seem concerned with either letting the size of the target vary, or ensuring that all rectangles from the available set are used. Here, we just want to take exactly the same approach as in the beforementioned question. That is, we

have a target rectangle of fixed size,
we have a collection of rectangles we can place, placement is without overlap, and each rectangle can be placed multiple times (although the case where each rectangle can only be placed once is also of interest, and can be used to solve the other case by having each rectangle repeated sufficiently many times), but we don't have to use all of them (unlike in e.g. what Korf calls the "containment problem"),
we want to minimize the area of the gap after placement.

For example, if the target rectangle is 25x25 and the placeable rectangles are 5x4, 6x7, and 3x10, the area of the gap in an optimal solution is 13. This can be verified by using e.g. the integer programming solution with inputs of the form

dims = [(5, 4), (6, 7), (3, 10)]
tiles = [[(x, y) for x in range(dimx) for y in range(dimy)] for dimx, dimy in dims]

A solution is found relatively quickly but this approach struggles with larger inputs (with, say, the same placeable rectangles but a 89x89 target). Can we do better?

One thing I considered was using Z3 and introducing variables corresponding to the location of each rectangle (as opposed to having variables for each location and each rectangle type, thus reducing the total number of variables at the expense of making their ranges larger and constraints more complicated). I.e. do something like the below. This works but is much slower than the answers in the linked question.

from z3 import *

# Example input
dims = [(5, 4), (5, 4), (6, 7), (3, 10)]
n = len(dims)
areas = [x * y for x, y in dims]
rows = 20
cols = 20

s = Optimize()

# Introduce a few variables: active[i] is True iff rectangle i is in use,
# and (xs[i], ys[i]) are its coordinates
actives = [Bool(f'active{i}') for i in range(n)]
xs = [Int(f'x{i}') for i in range(n)]
ys = [Int(f'y{i}') for i in range(n)]

# Require that all rectangles are within the bounds
for i, (x, y) in enumerate(dims):
    s.add(xs[i] >= 0, xs[i] <= rows - x)
    s.add(ys[i] >= 0, ys[i] <= cols - y)

# Help the solver a bit by breaking the symmetry: If two consecutive
# rectangles are identical, only use the second if the first is in use
for i, (d1, d2) in enumerate(zip(dims, dims[1:])):
    if d1 == d2:
        s.add(Implies(actives[i + 1], actives[i]))

# Disallow overlap
for i in range(n):
    for j in range(i + 1, n):
        (xi, yi), (xj, yj) = dims[i], dims[j]
        # Introduce a variable that's true iff rectangles i and j do not overlap
        freex = Or(xs[i] >= xs[j] + xj, xs[j] >= xs[i] + xi)
        freey = Or(ys[i] >= ys[j] + yj, ys[j] >= ys[i] + yi)
        free = Or(freex, freey)
        # If rectangles i and j are in use, require them to not overlap
        s.add(Implies(And(actives[i], actives[j]), free))

# The objective is the sum of the areas of all included rectangles
obj = Sum([actives[i] * areas[i] for i in range(n)])
s.maximize(obj)
s.check()

Here, an alternative to using Optimizer is to use Solver and keep adding constraints based on the current best solution until none exist; i.e. replace the last three lines with:

while s.check() != unsat:
    model = s.model()
    covered = sum(is_true(model[actives[i]]) * areas[i] for i in range(n))
    print('Objective:', rows*cols - covered)
    s.add(obj > covered)

Using the equivalent-ish formulation with MIP solver Xpress also doesn't help much:

actives = [xp.var(vartype=xp.binary) for i in range(n)]
xs = [xp.var(lb=0, ub=rows - dims[i][0], vartype=xp.integer) for i in range(n)]
ys = [xp.var(lb=0, ub=cols - dims[i][1], vartype=xp.integer) for i in range(n)]
p.addVariable(actives, xs, ys)

for i, (d1, d2) in enumerate(zip(dims, dims[1:])):
    if d1 == d2:
        p.addConstraint(actives[i] >= actives[i+1])

for i in range(n):
    for j in range(i + 1, n):
        (xi, yi), (xj, yj) = dims[i], dims[j]
        c1 = xs[i] - xs[j] - xj
        c2 = xs[j] - xs[i] - xi
        c3 = ys[i] - ys[j] - yj
        c4 = ys[j] - ys[i] - yi
        # Use max to represent disjunction
        p.addConstraint(xp.max(-actives[i], -actives[j], c1, c2, c3, c4) >= 0)

obj = -xp.Sum([areas[i] * actives[i] for i in range(n)])
p.setObjective(obj)
p.solve()

Recommended topics

Hot tags