投稿コンテンツ
デモアプリ
典型問題
マシン
コンテンツ一覧
ログインユーザー登録

Rectangle Packing Problem (矩形パッキング問題)

Kotaro Terada
#PyQUBO#矩形パッキング

Rectangle Packing Problem (矩形パッキング問題) のイジング実装

概要

矩形パッキング問題 (長方形詰め込み問題) をイジングモデルに定式化して解きます。本記事では、式の構築に PyQUBO、アニーリングに dwave-neal、Sequence-pairの面積計算・フロアプランの可視化に rectangle-packing-solver を利用しています。

ベースとなっている論文

  • [1] K. Terada, D. Oku, S. Kanamaru, S. Tanaka, M. Hayashi, M. Yamaoka, M. Yanagisawa, and N. Togawa, "An Ising model mapping to solve rectangle packing problem," in Proc. of 2018 International Symposium on VLSI Design, Automation and Test (VLSI-DAT), Hsinchu, Taiwan, Apr. 2018.

問題

論文 [1] の例題、矩形数 $N = 4$ の問題 (inst4) を取り上げる。面積の最適解は 56。

In [1]:
from IPython.display import Image, display
display(Image(url='https://user-images.githubusercontent.com/1787125/152990869-c841d42c-5acc-415c-84b1-dcb425d1d068.png', width=500))
In [2]:
rectangles = [
    (4, 6),  # (幅, 高さ)
    (4, 4),
    (2, 3),
    (5, 1),
]
N = len(rectangles)
print('N:', N)

# 各矩形の縦横の正規化
W_sum = sum([rectangles[i][0] for i in range(N)])
H_sum = sum([rectangles[i][1] for i in range(N)])
W, H = [], []
for i in range(N):
    W.append(rectangles[i][0] / (W_sum + H_sum))
    H.append(rectangles[i][1] / (W_sum + H_sum))
print('W:', W)
print('H:', H)
N: 4
W: [0.13793103448275862, 0.13793103448275862, 0.06896551724137931, 0.1724137931034483]
H: [0.20689655172413793, 0.13793103448275862, 0.10344827586206896, 0.034482758620689655]

準備

必要パッケージのインポート

In [3]:
import pyqubo
import neal
import rectangle_packing_solver as rps

import itertools
import sys

エネルギー関数の構築

矩形パッキング問題を解くイジングモデルのエネルギー関数は以下の式で定式化される。

$$H_{\textrm{RPP}} = \alpha H_{\textrm{A}} + \beta H_{\textrm{B}} + \gamma H_{\textrm{G}} + \delta H_{\textrm{area}} + \epsilon H_{\textrm{longest}} $$

それぞれの項の詳細な意味は論文[1]をご参照ください。

$\alpha$, $\beta$, $\gamma$, $\delta$, $\epsilon$ は制約重み係数 (ハイパーパラメータ) である。

In [4]:
x = pyqubo.Array.create('x', (N, N, N), 'BINARY')
y = pyqubo.Array.create('y', (N, N, N), 'BINARY')
z = pyqubo.Array.create('z', (N, N, N), 'BINARY')

$H_{\textrm{A}}$ の構築。

In [5]:
h_a_i = 0.0
for i in range(N):
    h_a_i += pyqubo.Constraint(
        (sum(x[i, j, k] for j, k in itertools.product(range(N), repeat=2)) - 1)**2, 
        label="h_a_i{}".format(i)
    )
h_a_j = 0.0
for j in range(N):
    h_a_j += pyqubo.Constraint(
        (sum(x[i, j, k] for i, k in itertools.product(range(N), repeat=2)) - 1)**2, 
        label="h_a_j{}".format(j)
    )
h_a_k = 0.0
for k in range(N):
    h_a_k += pyqubo.Constraint(
        (sum(x[i, j, k] for i, j in itertools.product(range(N), repeat=2)) - 1)**2, 
        label="h_a_k{}".format(k)
    )

$H_{\textrm{B}}$ の構築。

In [6]:
h_b = 0.0
for (i, j, k) in itertools.product(range(N), repeat=3):
    h_b += -2 * x[i, j, k] * y[i, j, k] + x[i, j, k] + y[i, j, k]
    h_b += -2 * x[i, j, k] * z[i, j, k] + x[i, j, k] + z[i, j, k]

$H_{\textrm{G}}$ の構築。

In [7]:
h_g_h = 0.0
for k1 in range(N):
    for (i1, j1) in itertools.product(range(N), repeat=2):
        for (i2, j2) in itertools.product(list(range(N)), repeat=2):
            if (i1 < i2) and (j1 > j2):
                for k2 in range(N):
                    h_g_h += y[i1, j1, k1] * y[i2, j2, k2]
        for (i2, j2) in itertools.product(list(range(N)), repeat=2):
            if (i1 > i2) and (j1 < j2):
                for k2 in range(N):
                    h_g_h += y[i1, j1, k1] * y[i2, j2, k2]

h_g_v = 0.0
for k1 in range(N):
    for (i1, j1) in itertools.product(range(N), repeat=2):
        for (i2, j2) in itertools.product(list(range(N)), repeat=2):
            if (i1 > i2) and (j1 > j2):
                for k2 in range(N):
                    h_g_v += z[i1, j1, k1] * z[i2, j2, k2]
        for (i2, j2) in itertools.product(list(range(N)), repeat=2):
            if (i1 < i2) and (j1 < j2):
                for k2 in range(N):
                    h_g_v += z[i1, j1, k1] * z[i2, j2, k2]

$H_{\textrm{area}}$ の構築。

In [8]:
h_area = 0.0
for (yi, yj, yk) in itertools.product(range(N), repeat=3):
    for (zi, zj, zk) in itertools.product(range(N), repeat=3):
        h_area += W[yk] * H[zk] * y[yi, yj, yk] * z[zi, zj, zk]

$H_{\textrm{longest}}$ の構築。

In [9]:
h_longest = 0.0
for (i, j, k) in itertools.product(range(N), repeat=3):
    h_longest -= W[k] * y[i, j, k]
    h_longest -= H[k] * z[i, j, k]

各項の重み係数値の設定。

論文に従い $\alpha = 2, \beta = 1, \gamma = 1, \delta = 1, \epsilon = 1$ とする。

In [10]:
alpha = pyqubo.Placeholder('alpha')
beta = pyqubo.Placeholder('beta')
gamma = pyqubo.Placeholder('gamma')
delta = pyqubo.Placeholder('delta')
epsilon = pyqubo.Placeholder('epsilon')

hamiltonian = (
    alpha * (h_a_i + h_a_j + h_a_k) +
    beta * h_b +
    gamma * (h_g_h + h_g_v) +
    delta * h_area + 
    epsilon * h_longest
)
model = hamiltonian.compile()
feed_dict = {'alpha': 2.0, 'beta': 1.0, 'gamma': 1.0, 'delta': 1.0, 'epsilon': 1.0}
bqm = model.to_bqm(feed_dict=feed_dict)

アニーリング実行

dwave-neal を利用したローカルでのソフトウェアアニーリングの実行。 アニーリングステップ数とシード値を固定しているので毎回同じサンプリング結果になるはずです。

In [11]:
sa = neal.SimulatedAnnealingSampler()
sampleset = sa.sample(bqm, num_reads=100, num_sweeps=2000, seed=1234).aggregate()

アニーリング結果の確認

In [12]:
# Seq-pair を求めるヘルパー関数
def get_seq_pair(state):
    oblique_grid = [[-1 for _ in range(N)] for _ in range(N)]
    for k in range(N):
        for (i, j) in itertools.product(range(N), repeat=2):
            if state['x[{}][{}][{}]'.format(i, j, k)] == 1:
                oblique_grid[i][j] = k
                break

    gp, gn= [], []
    for j in range(N):
        for i in range(N):
            if oblique_grid[i][j] != -1:
                gp.append(oblique_grid[i][j])
                break
    for i in range(N):
        for j in range(N):
            if oblique_grid[i][j] != -1:
                gn.append(oblique_grid[i][j])
                break

    return (gp, gn)
In [13]:
print('We have {} samples\n'.format(len(sampleset)))

# 各サンプルの結果
best_area = sys.float_info.max
best_seq_pair = None
for idx, sample in enumerate(sampleset.record):
    state = dict(zip(sampleset.variables, sample.sample))
    seq_pair = get_seq_pair(state)
    if (len(seq_pair[0]) != 4) or (len(seq_pair[1]) != 4):
        continue

    # `rectangle-packing-solver` パッケージを利用して面積計算
    rps_seqpair = rps.SequencePair(pair=seq_pair)
    floorplan = rps_seqpair.decode(problem=rps.Problem(rectangles=rectangles))

    if floorplan.area < best_area:
        best_area = floorplan.area
        best_seq_pair = seq_pair
        print(f'#{idx}  seq_pair: {seq_pair}  area: {floorplan.area}')
We have 93 samples

#0  seq_pair: ([2, 0, 3, 1], [1, 3, 0, 2])  area: 70
#3  seq_pair: ([1, 3, 2, 0], [2, 0, 3, 1])  area: 66
#30  seq_pair: ([0, 1, 3, 2], [3, 0, 2, 1])  area: 56

面積最小の解をピックアップしてフロアプランを表示

In [14]:
sequence_pair = rps.SequencePair(pair=best_seq_pair)
floorplan = sequence_pair.decode(problem=rps.Problem(rectangles=rectangles))
solution = rps.Solution(sequence_pair=sequence_pair, floorplan=floorplan)
rps.Visualizer().visualize(solution=solution, path='./floorplan.png')

display(Image('./floorplan.png', width=500))

最適解が得られました。

Kotaro Terada @kotaro

このコンテンツにコメントはありません。

empty
ユーザーのみコメントを投稿できます。
ユーザー登録