网格上的二维装箱

2D bin packing on a grid

我有一个 n × m 网格和一个 polyominos 的集合。我想知道是否可以将它们打包到网格中:不允许重叠或旋转。

我希望像大多数包装问题一样,这个版本是 NP-hard 并且难以近似,所以我不期待任何疯狂的事情,而是一种算法可以在 25 × 25 左右的网格上找到合理的包装并且相当综合10×10左右就好了。 (我的图块大部分是四联纸牌——四个方块——但它们可能有 5-9 个以上的方块。)

我会接受任何人提供的任何东西:算法、论文、可以改编的现有程序。

我不知道它是否对你有用,但我在 Python 中编写了一个小的粗略框架。它还没有放置 polyminos,但功能在那里 - 检查死空空间是原始的,但需要更好的方法。话又说回来,也许都是垃圾...

import functools
import itertools

M = 4 # x
N = 5 # y

field = [[9999]*(N+1)]+[[9999]+[0]*N+[9999] for _ in range(M)]+[[9999]*(N+1)]

def field_rd(p2d):
    return field[p2d[0]+1][p2d[1]+1]

def field_add(p2d,val):
    field[p2d[0]+1][p2d[1]+1] += val

def add2d(p,k):
    return p[0]+k[0],p[1]+k[1]

def norm(polymino_2d):
    x0,y0 = min(x for x,y in polymino_2d),min(y for x,y in polymino_2d)
    return tuple(sorted(map(lambda p: add2d(p,(-x0,-y0)), polymino_2d)))

def create_cutoff(occupied):
    """Receive a polymino and create the outer area of squares which could be cut off by a placement of this polymino"""
    cutoff = set(itertools.chain.from_iterable(map(lambda p: add2d(p,(x,y)),occupied) for (x,y) in [(-1,0),(1,0),(0,-1),(0,1)])) #(-1,-1),(-1,0),(-1,1),(0,1),(1,1),(1,0),(1,-1)]))
    return tuple(cutoff.difference(occupied))

def is_occupied(p2d):
    return field_rd(p2d) == 0

def is_cutoff(p2d):
    return not is_occupied(p2d) and all(map(is_occupied,map(lambda p: add2d(p,p2d),[(-1,0),(1,0),(0,-1),(0,1)])))

def polym_colliding(p2d,occupied):
    return any(map(is_occupied,map(lambda p: add2d(p,p2d),occupied)))

def polym_cutoff(p2d,cutoff):
    return any(map(is_cutoff,map(lambda p: add2d(p,p2d),cutoff)))

def put(p2d,occupied,polym_nr):
    for p in occupied:
        field_add(add2d(p2d,p),polym_nr)

def remove(p2d,occupied,polym_nr):
    for p in polym:
        field_add(add2d(p2d,p),-polym_nr)

def place(p2d,polym_nr):
    """Try to place a polymino at point p2d. If it fits without cutting off unreachable single cells return True else False"""
    occupied = polym[polym_nr][0]
    if polym_colliding(p2d,occupied):
        return False
    put(p2d,occupied,polym_nr)
    cutoff = polym[polym_nr][1]
    if polym_cutoff(p2d,cutoff):
        remove(p2d,occupied,polym_nr)
        return False
    return True

def NxM_array(N,M):
    return [[0]*N for _ in range(M)]


def generate_all_polyminos(n):
    """Create all polyminos with size n"""
    def gen_recur(polymino,i,result):
        if i > 1:
            new_pts = set(itertools.starmap(add2d,itertools.product(polymino,[(-1,0),(1,0),(0,-1),(0,1)])))
            new_pts = new_pts.difference(polymino)
            for p in new_pts:
                gen_recur(polymino.union({p}),i-1,result)
        else:
            result.add(norm(polymino))
    #---------------------------------------
    all_polyminos = set()
    gen_recur({(0,0)},n,all_polyminos)
    return all_polyminos

print("All possible Tetris blocks (all orientations): ",generate_all_polyminos(4))

一种方法是使用整数规划。我将使用 python pulp 包来实现它,尽管几乎所有编程语言都可以使用这些包。

基本思想是为每个图块的每个可能放置位置定义一个决策变量。如果决策变量取值 1,则将其关联的图块放在那里。如果它的值为 0,则它不会放在那里。因此,objective 是最大化决策变量的总和乘以变量块中的方块数——这对应于在棋盘上放置最大可能的方块数。

我的代码实现了两个约束:

  • 每个方块只能放置一次(下面我们会放宽这个限制)
  • 每个方格最多只能有一个方块

这是 4x5 网格上一组五个固定四联骨牌的输出:

import itertools
import pulp
import string

def covered(tile, base):
    return {(base[0] + t[0], base[1] + t[1]): True for t in tile}

tiles = [[(0,0), (1,0), (0,1), (0,2)],
         [(0,0), (1,0), (2,0), (3,0)],
         [(1,0), (0,1), (1,1), (2,0)],
         [(0,0), (1,0), (0,1), (1,1)],
         [(1,0), (0,1), (1,1), (2,1)]]
rows = 25
cols = 25
squares = {x: True for x in itertools.product(range(rows), range(cols))}
vars = list(itertools.product(range(rows), range(cols), range(len(tiles))))
vars = [x for x in vars if all([y in squares for y in covered(tiles[x[2]], (x[0], x[1])).keys()])]
x = pulp.LpVariable.dicts('tiles', vars, lowBound=0, upBound=1, cat=pulp.LpInteger)
mod = pulp.LpProblem('polyominoes', pulp.LpMaximize)
# Objective value is number of squares in tile
mod += sum([len(tiles[p[2]]) * x[p] for p in vars])
# Don't use any shape more than once
for tnum in range(len(tiles)):
    mod += sum([x[p] for p in vars if p[2] == tnum]) <= 1
# Each square can be covered by at most one shape
for s in squares:
    mod += sum([x[p] for p in vars if s in covered(tiles[p[2]], (p[0], p[1]))]) <= 1
# Solve and output
mod.solve()
out = [['-'] * cols for rep in range(rows)]
chars = string.ascii_uppercase + string.ascii_lowercase
numset = 0
for p in vars:
    if x[p].value() == 1.0:
        for off in tiles[p[2]]:
            out[p[0] + off[0]][p[1] + off[1]] = chars[numset]
        numset += 1
for row in out:
    print(''.join(row))

得到如下最优解:

AAAB-
A-BBC
DDBCC
DD--C

如果我们允许重复(注释掉限制为每个形状一个副本的约束),那么我们可以完全平铺网格:

ABCDD
ABCDD
ABCEE
ABCEE

对于 10x10 网格,它几乎可以立即工作:

ABCCDDEEFF
ABCCDDEEFF
ABGHHIJJKK
ABGHHIJJKK
LLGMMINOPP
LLGMMINOPP
QQRRSTNOUV
QQRRSTNOUV
WWXXSTYYUV
WWXXSTYYUV

该代码在 100 秒的运行时间中获得了 25x25 网格的最佳解决方案,但不幸的是,我的输出代码没有足够的字母和数字来打印解决方案。

这是一个类似于原型的 SAT-solver 方法,它解决了:

  • a-priori 固定 polyomino 模式(参见代码中的 Constants / Input
    • 如果允许旋转,则必须将旋转的棋子添加到集合中
  • 每个 polyomino 都可以放置 0-inf 次
  • 除此之外没有计分机制:
    • 未覆盖的图块数量已最小化!

考虑到用于组合优化的经典现成方法 (SAT, CP, MIP),这个方法可能会扩展得最好(有根据的猜测)。在设计自定义启发式时也很难被击败!

如果需要,这些幻灯片提供了一些基于 practical introduction to SAT-solvers in practice. Here we are using CDCL 的求解器,它们是 完整的(如果有的话,总能在有限的时间内找到解决方案;永远是如果有none就可以证明有限时间内无解;当然记忆力也有一定的作用!)。

通常很难合并更复杂的(线性)每图块评分函数。这是 (M)IP 方法可以做得更好的地方。但就纯搜索而言,SAT 求解通常要快得多。

我的 polyomino 集的 N=25 问题需要 ~ 1 秒(并且可以很容易地将其并行化到多个粒度级别 -> SAT 求解器(线程-param) vs. outer-loop; 后者将在后面解释).

当然以下内容成立:

  • 因为这是一个 NP-hard 问题,所以会有容易和不容易的实例
  • 我没有对许多不同的多联骨牌进行科学基准测试
    • 有些集合比其他集合更容易求解是可以预料的
  • 这是 一个可能的 SAT 公式(不是最简单的!)无穷多个
    • 每个公式都有优点和缺点

想法

一般方法是创建一个 决策问题 并将其转换为 CNF, which is then solved by highly efficient SAT-solvers (here: cryptominisat; CNF will be in DIMCAS-CNF format),这将用作黑盒求解器(无参数调整! ).

由于我们的目标是优化填充图块的数量并且我们正在使用决策问题,因此我们需要一个外循环,添加最小图块使用约束并尝试解决它。如果不成功,请减少此数字。所以通常我们会多次调用 SAT 求解器(从头开始!)。

有许多不同的 CNF 可能的公式/转换。在这里,我们使用(二进制)决策变量 X 表示 放置 placement 是一个类似于 polyomino, x_index, y_index 的元组(该索引标记了某些模式的左上角字段)。变量个数和所有多联骨牌的可能放置个数之间存在一对一的映射。

核心思想是:在 space 所有可能的放置组合中搜索一个解决方案,这不会使某些约束失效。

此外,我们还有决策变量 Y,它指示正在填充的图块。有M*N个这样的变数。

访问所有可能的位置时,很容易计算每个图块索引 (M*N) 的碰撞集。给定一些固定的图块,我们可以检查哪些位置可以填充这个,并将问题限制为只有 select <=1 个。这在 X 上有效。在 (M)IP 世界中,这可能被称为碰撞的凸包。

n<=k-约束在 SAT 求解中无处不在,许多不同的公式都是可能的。朴素编码一般需要指数数量的子句,这很容易变得不可行。使用新变量,可以进行许多变量子句权衡(参见 Tseitin-encoding)。我正在重用一个(旧代码;我的代码只有 python2 的唯一原因),它在过去对我很有用。它基于将基于硬件的反逻辑描述到 CNF 中,并提供良好的经验和理论性能(参见论文)。当然还有很多选择。

此外,我们需要强制 SAT 求解器不要使所有变量都为负。我们必须添加描述以下内容的约束(这是一种方法):

  • 如果使用某个字段:必须至少有一个放置处于活动状态 (poly + x + y),这会导致覆盖该字段!
    • 这是一个基本的逻辑蕴涵,很容易表述为一个潜在的大逻辑或

然后只缺少核心循环,尝试填充N个字段,然后N-1直到成功。这再次使用前面提到的 n<=k 公式。

代码

这是 python2 代码,需要 SAT 求解器 cryptominisat 5 在脚本 运行 来自的目录中。

我也在使用 python 优秀的科学堆栈中的工具。

# PYTHON 2!
import math
import copy
import subprocess
import numpy as np
import matplotlib.pyplot as plt      # plotting-only
import seaborn as sns                # plotting-only
np.set_printoptions(linewidth=120)   # more nice console-output

""" Constants / Input
        Example: 5 tetrominoes; no rotation """
M, N = 25, 25
polyominos = [np.array([[1,1,1,1]]),
              np.array([[1,1],[1,1]]),
              np.array([[1,0],[1,0], [1,1]]),
              np.array([[1,0],[1,1],[0,1]]),
              np.array([[1,1,1],[0,1,0]])]

""" Preprocessing
        Calculate:
        A: possible placements
        B: covered positions
        C: collisions between placements
"""
placements = []
covered = []
for p_ind, p in enumerate(polyominos):
    mP, nP = p.shape
    for x in range(M):
        for y in range(N):
            if x + mP <= M:          # assumption: no zero rows / cols in each p
                if y + nP <= N:      # could be more efficient
                    placements.append((p_ind, x, y))
                    cover = np.zeros((M,N), dtype=bool)
                    cover[x:x+mP, y:y+nP] = p
                    covered.append(cover)                           
covered = np.array(covered)

collisions = []
for m in range(M):
    for n in range(N):
        collision_set = np.flatnonzero(covered[:, m, n])
        collisions.append(collision_set)

""" Helper-function: Cardinality constraints """
# K-ARY CONSTRAINT GENERATION
# ###########################
# SINZ, Carsten. Towards an optimal CNF encoding of boolean cardinality constraints.
# CP, 2005, 3709. Jg., S. 827-831.

def next_var_index(start):
    next_var = start
    while(True):
        yield next_var
        next_var += 1

class s_index():
    def __init__(self, start_index):
        self.firstEnvVar = start_index

    def next(self,i,j,k):
        return self.firstEnvVar + i*k +j

def gen_seq_circuit(k, input_indices, next_var_index_gen):
    cnf_string = ''
    s_index_gen = s_index(next_var_index_gen.next())

    # write clauses of first partial sum (i.e. i=0)
    cnf_string += (str(-input_indices[0]) + ' ' + str(s_index_gen.next(0,0,k)) + ' 0\n')
    for i in range(1, k):
        cnf_string += (str(-s_index_gen.next(0, i, k)) + ' 0\n')

    # write clauses for general case (i.e. 0 < i < n-1)
    for i in range(1, len(input_indices)-1):
        cnf_string += (str(-input_indices[i]) + ' ' + str(s_index_gen.next(i, 0, k)) + ' 0\n')
        cnf_string += (str(-s_index_gen.next(i-1, 0, k)) + ' ' + str(s_index_gen.next(i, 0, k)) + ' 0\n')
        for u in range(1, k):
            cnf_string += (str(-input_indices[i]) + ' ' + str(-s_index_gen.next(i-1, u-1, k)) + ' ' + str(s_index_gen.next(i, u, k)) + ' 0\n')
            cnf_string += (str(-s_index_gen.next(i-1, u, k)) + ' ' + str(s_index_gen.next(i, u, k)) + ' 0\n')
        cnf_string += (str(-input_indices[i]) + ' ' + str(-s_index_gen.next(i-1, k-1, k)) + ' 0\n')

    # last clause for last variable
    cnf_string += (str(-input_indices[-1]) + ' ' + str(-s_index_gen.next(len(input_indices)-2, k-1, k)) + ' 0\n')

    return (cnf_string, (len(input_indices)-1)*k, 2*len(input_indices)*k + len(input_indices) - 3*k - 1)

def gen_at_most_n_constraints(vars, start_var, n):
    constraint_string = ''
    used_clauses = 0
    used_vars = 0
    index_gen = next_var_index(start_var)
    circuit = gen_seq_circuit(n, vars, index_gen)
    constraint_string += circuit[0]
    used_clauses += circuit[2]
    used_vars += circuit[1]
    start_var += circuit[1]

    return [constraint_string, used_clauses, used_vars, start_var]

def parse_solution(output):
    # assumes there is one
    vars = []
    for line in output.split("\n"):
        if line:
            if line[0] == 'v':
                line_vars = list(map(lambda x: int(x), line.split()[1:]))
                vars.extend(line_vars)
    return vars

def solve(CNF):
    p = subprocess.Popen(["cryptominisat5.exe"], stdin=subprocess.PIPE, stdout=subprocess.PIPE)
    result = p.communicate(input=CNF)[0]
    sat_line = result.find('s SATISFIABLE')
    if sat_line != -1:
        # solution found!
        vars = parse_solution(result)
        return True, vars
    else:
        return False, None

""" SAT-CNF: BASE """
X = np.arange(1, len(placements)+1)                                     # decision-vars
                                                                        # 1-index for CNF
Y = np.arange(len(placements)+1, len(placements)+1 + M*N).reshape(M,N)
next_var = len(placements)+1 + M*N                                      # aux-var gen
n_clauses = 0

cnf = ''                                                                # slow string appends
                                                                        # int-based would be better
# <= 1 for each collision-set
for cset in collisions:
    constraint_string, used_clauses, used_vars, next_var = \
        gen_at_most_n_constraints(X[cset].tolist(), next_var, 1)
    n_clauses += used_clauses
    cnf += constraint_string

# if field marked: one of covering placements active
for x in range(M):
    for y in range(N):
        covering_placements = X[np.flatnonzero(covered[:, x, y])]  # could reuse collisions
        clause = str(-Y[x,y])
        for i in covering_placements:
            clause += ' ' + str(i)
        clause += ' 0\n'
        cnf += clause
        n_clauses += 1

print('BASE CNF size')
print('clauses: ', n_clauses)
print('vars: ', next_var - 1)

""" SOLVE in loop -> decrease number of placed-fields until SAT """
print('CORE LOOP')
N_FIELD_HIT = M*N
while True:
    print(' N_FIELDS >= ', N_FIELD_HIT)
    # sum(y) >= N_FIELD_HIT
    # == sum(not y) <= M*N - N_FIELD_HIT
    cnf_final = copy.copy(cnf)
    n_clauses_final = n_clauses

    if N_FIELD_HIT == M*N:  # awkward special case
        constraint_string = ''.join([str(y) + ' 0\n' for y in Y.ravel()])
        n_clauses_final += N_FIELD_HIT
    else:
        constraint_string, used_clauses, used_vars, next_var = \
            gen_at_most_n_constraints((-Y).ravel().tolist(), next_var, M*N - N_FIELD_HIT)
        n_clauses_final += used_clauses

    n_vars_final = next_var - 1
    cnf_final += constraint_string
    cnf_final = 'p cnf ' + str(n_vars_final) + ' ' + str(n_clauses) + \
        ' \n' + cnf_final  # header

    status, sol = solve(cnf_final)
    if status:
        print(' SOL found: ', N_FIELD_HIT)

        """ Print sol """
        res = np.zeros((M, N), dtype=int)
        counter = 1
        for v in sol[:X.shape[0]]:
            if v>0:
                p, x, y = placements[v-1]
                pM, pN = polyominos[p].shape
                poly_nnz = np.where(polyominos[p] != 0)
                x_inds, y_inds = x+poly_nnz[0], y+poly_nnz[1]
                res[x_inds, y_inds] = p+1
                counter += 1
        print(res)

        """ Plot """
        # very very ugly code; too lazy
        ax1 = plt.subplot2grid((5, 12), (0, 0), colspan=11, rowspan=5)
        ax_p0 = plt.subplot2grid((5, 12), (0, 11))
        ax_p1 = plt.subplot2grid((5, 12), (1, 11))
        ax_p2 = plt.subplot2grid((5, 12), (2, 11))
        ax_p3 = plt.subplot2grid((5, 12), (3, 11))
        ax_p4 = plt.subplot2grid((5, 12), (4, 11))
        ax_p0.imshow(polyominos[0] * 1, vmin=0, vmax=5)
        ax_p1.imshow(polyominos[1] * 2, vmin=0, vmax=5)
        ax_p2.imshow(polyominos[2] * 3, vmin=0, vmax=5)
        ax_p3.imshow(polyominos[3] * 4, vmin=0, vmax=5)
        ax_p4.imshow(polyominos[4] * 5, vmin=0, vmax=5)
        ax_p0.xaxis.set_major_formatter(plt.NullFormatter())
        ax_p1.xaxis.set_major_formatter(plt.NullFormatter())
        ax_p2.xaxis.set_major_formatter(plt.NullFormatter())
        ax_p3.xaxis.set_major_formatter(plt.NullFormatter())
        ax_p4.xaxis.set_major_formatter(plt.NullFormatter())
        ax_p0.yaxis.set_major_formatter(plt.NullFormatter())
        ax_p1.yaxis.set_major_formatter(plt.NullFormatter())
        ax_p2.yaxis.set_major_formatter(plt.NullFormatter())
        ax_p3.yaxis.set_major_formatter(plt.NullFormatter())
        ax_p4.yaxis.set_major_formatter(plt.NullFormatter())

        mask = (res==0)
        sns.heatmap(res, cmap='viridis', mask=mask, cbar=False, square=True, linewidths=.1, ax=ax1)
        plt.tight_layout()
        plt.show()
        break

    N_FIELD_HIT -= 1  # binary-search could be viable in some cases
                      # but beware the empirical asymmetry in SAT-solvers:
                      #    finding solution vs. proving there is none!

输出控制台

BASE CNF size
('clauses: ', 31509)
('vars: ', 13910)
CORE LOOP
(' N_FIELDS >= ', 625)
(' N_FIELDS >= ', 624)
(' SOL found: ', 624)
[[3 2 2 2 2 1 1 1 1 1 1 1 1 2 2 1 1 1 1 1 1 1 1 2 2]
 [3 2 2 2 2 1 1 1 1 1 1 1 1 2 2 2 2 2 2 1 1 1 1 2 2]
 [3 3 3 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 1 1 1 1 2 2]
 [2 2 3 1 1 1 1 1 1 1 1 2 2 2 2 1 1 1 1 2 2 2 2 2 2]
 [2 2 3 3 3 2 2 2 2 2 2 2 2 2 2 1 1 1 1 2 2 2 2 2 2]
 [1 1 1 1 3 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 1 1 2 2]
 [1 1 1 1 3 3 3 2 2 1 1 1 1 2 2 2 2 2 2 2 2 1 1 1 1]
 [2 2 1 1 1 1 3 2 2 2 2 2 2 2 2 1 1 1 1 2 2 2 2 2 2]
 [2 2 2 2 2 2 3 3 3 2 2 2 2 1 1 1 1 2 2 2 2 2 2 2 2]
 [2 2 2 2 2 2 2 2 3 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2]
 [2 2 1 1 1 1 2 2 3 3 3 2 2 2 2 2 2 1 1 1 1 2 2 2 2]
 [1 1 1 1 1 1 1 1 2 2 3 2 2 1 1 1 1 1 1 1 1 1 1 1 1]
 [2 2 3 1 1 1 1 3 2 2 3 3 4 1 1 1 1 2 2 1 1 1 1 2 2]
 [2 2 3 1 1 1 1 3 1 1 1 1 4 4 3 2 2 2 2 1 1 1 1 2 2]
 [2 2 3 3 5 5 5 3 3 1 1 1 1 4 3 2 2 1 1 1 1 1 1 1 1]
 [2 2 2 2 4 5 1 1 1 1 1 1 1 1 3 3 3 2 2 1 1 1 1 2 2]
 [2 2 2 2 4 4 2 2 1 1 1 1 1 1 1 1 3 2 2 1 1 1 1 2 2]
 [2 2 2 2 3 4 2 2 2 2 2 2 1 1 1 1 3 3 3 2 2 2 2 2 2]
 [3 4 2 2 3 5 5 5 2 2 2 2 1 1 1 1 2 2 3 2 2 2 2 2 2]
 [3 4 4 3 3 3 5 5 5 5 1 1 1 1 2 2 2 2 3 3 3 2 2 2 2]
 [3 3 4 3 1 1 1 1 5 1 1 1 1 4 2 2 2 2 2 2 3 2 2 2 2]
 [2 2 3 3 3 1 1 1 1 1 1 1 1 4 4 4 2 2 2 2 3 3 0 2 2]
 [2 2 3 1 1 1 1 1 1 1 1 5 5 5 4 4 4 1 1 1 1 2 2 2 2]
 [2 2 3 3 1 1 1 1 1 1 1 1 5 5 5 5 4 1 1 1 1 2 2 2 2]
 [2 2 1 1 1 1 1 1 1 1 1 1 1 1 5 1 1 1 1 1 1 1 1 2 2]]

输出图

此参数化不能覆盖一个字段!

具有更多模式集的其他示例

Square M=N=61(素数 -> 直觉:更难),其中 base-CNF 有 450.723 个子句和 185.462 个变量。有最佳包装!

Non-square M,N =83,131(双素数),其中 base-CNF 有 1.346.511 个子句和 553.748 个变量。有最佳包装!