パッキング問題

ビンパッキング問題， カッティングストック問題， 2次元パッキング問題

準備

from IPython.display import Image, YouTubeVideo

YouTubeVideo("iLgDpPc-Sl0")

ビンパッキング問題

ビンパッキング問題(bin packing problem; 箱詰め問題)は，以下のように定義される問題である．

\(n\)個のアイテムから成る有限集合 \(N\) とサイズ（容量） \(B\) のビンが無限個準備されている． 個々のアイテム \(i \in N\) のサイズ \(0 \leq w_i \leq B\) は分かっているものとする． これら \(n\) 個のアイテムを，サイズ \(B\) のビンに詰めることを考えるとき， 必要なビンの数を最小にするような詰めかたを求めよ．

定式化

ビンの数の上限 \(U\) が与えられているものとする． アイテム \(i\) をビン \(j\) に詰めるとき \(1\) になる変数 \(x_{ij}\) と， ビン \(j\) の使用の可否を表す変数 \(y_j\) を用いることによって，ビンパッキング問題は，以下の整数最適化問題として記述できる．

\[ \begin{array}{l l l} minimize & \displaystyle\sum_{j=1}^U y_{j} & \\ s.t. & \displaystyle\sum_{j=1}^{U} x_{ij} =1 & \forall i =1,2,\ldots,n \\ & \displaystyle\sum_{i=1}^n w_i x_{ij} \leq B y_j & \forall j=1,2,\ldots,U \\ & x_{ij} \leq y_j & \forall i=1,2,\ldots,n, j=1,2,\ldots,U \\ & x_{ij} \in \{ 0,1 \} & \forall i=1,2,\ldots,n, j=1,2,\ldots,U \\ & y_j \in \{ 0,1 \} & \forall j=1,2,\ldots,U \end{array} \]

この定式化は，大規模な実際問題を解く際には使うべきではない．データが整数の場合には，枝フロー変数を用いたより強い定式化が提案されている． 詳細は， 以下のサイトを参照されたい．

https://vpsolver.fdabrandao.pt/

Pythonによる求解

離散的な値をもつアイテムのサイズの例題を生成し，上の定式化で求解してみる．

def bpp(s, B):
    """bpp: Martello and Toth's model to solve the bin packing problem.
    Parameters:
        - s: list with item widths
        - B: bin capacity
    Returns a model, ready to be solved.
    """
    n = len(s)
    U = n   # upper bound of the number of bins; 本来なら後述の近似解法を用いる len(FFD(s, B)) 
    model = Model("bpp")
    x, y = {}, {}
    for i in range(n):
        for j in range(U):
            x[i, j] = model.addVar(vtype="B", name="x(%s,%s)" % (i, j))
    for j in range(U):
        y[j] = model.addVar(vtype="B", name="y(%s)" % j)
    model.update()

    # assignment constraints
    for i in range(n):
        model.addConstr(quicksum(x[i, j] for j in range(U)) == 1, "Assign(%s)" % i)

    # bin capacity constraints
    for j in range(U):
        model.addConstr(
            quicksum(s[i] * x[i, j] for i in range(n)) <= B * y[j], "Capac(%s)" % j
        )

    # tighten assignment constraints
    for j in range(U):
        for i in range(n):
            model.addConstr(x[i, j] <= y[j], "Strong(%s,%s)" % (i, j))

    model.setObjective(quicksum(y[j] for j in range(U)), GRB.MINIMIZE)

    model.update()
    model.__data = x, y
    return model, U


def solvebinPacking(s, B):
    """solvebinPacking: use an IP model to solve the in Packing Problem.

    Parameters:
        - s: list with item widths
        - b: bin capacity

    Returns a solution: list of lists, each of which with the items in a roll.
    """
    n = len(s)
    model, U = bpp(s, B)
    x, y = model.__data

    model.optimize()

    bins = [[] for i in range(U)]
    for (i, j) in x:
        if x[i, j].X > 0.5:
            bins[j].append(s[i])
    for i in range(bins.count([])):
        bins.remove([])
    for b in bins:
        b.sort()
    bins.sort()
    return bins


def Discretebniform(n=10, lb=1, ub=99, b=100):
    """Discretebniform: create random, uniform instance for the bin packing problem."""
    b = 100
    s = [0] * n
    for i in range(n):
        s[i] = random.randint(lb, ub)
    return s, b

random.seed(256)
s, B = Discretebniform(10, 1, 50, 100)
print("items:", s)
print("bin size:", B)

bins = solvebinPacking(s, B)
print(len(bins), "bins:")
print(bins)

items: [32, 20, 28, 24, 25, 3, 30, 50, 28, 16]
bin size: 100
Restricted license - for non-production use only - expires 2025-11-24
Gurobi Optimizer version 11.0.1 build v11.0.1rc0 (mac64[arm] - Darwin 23.0.0 23A344)

CPU model: Apple M2 Max
Thread count: 12 physical cores, 12 logical processors, using up to 12 threads

Optimize a model with 120 rows, 110 columns and 410 nonzeros
Model fingerprint: 0xf2f60471
Variable types: 0 continuous, 110 integer (110 binary)
Coefficient statistics:
  Matrix range     [1e+00, 1e+02]
  Objective range  [1e+00, 1e+00]
  Bounds range     [1e+00, 1e+00]
  RHS range        [1e+00, 1e+00]
Found heuristic solution: objective 6.0000000
Presolve time: 0.00s
Presolved: 120 rows, 110 columns, 410 nonzeros
Variable types: 0 continuous, 110 integer (110 binary)

Root relaxation: objective 2.560000e+00, 127 iterations, 0.00 seconds (0.00 work units)

    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

     0     0    2.56000    0    7    6.00000    2.56000  57.3%     -    0s
H    0     0                       3.0000000    2.56000  14.7%     -    0s
     0     0    2.56000    0    7    3.00000    2.56000  14.7%     -    0s

Explored 1 nodes (247 simplex iterations) in 0.02 seconds (0.00 work units)
Thread count was 12 (of 12 available processors)

Solution count 2: 3 6 

Optimal solution found (tolerance 1.00e-04)
Best objective 3.000000000000e+00, best bound 3.000000000000e+00, gap 0.0000%
3 bins:
[[3, 25, 28], [16, 24, 28, 32], [20, 30, 50]]

ヒューリスティクス

ビンの数の上界 \(U\) を求めるために，ヒューリスティクスを用いる．

first fit (FF) ヒューリスティクス

アイテムを \(1,2,\ldots,n\) の順にビンに詰めていく． このとき，アイテムは詰め込み可能な最小添字のビンに詰めるものとする． どのビンに入れてもサイズの上限 \(B\) を超えてしまうなら， 新たなビンを用意し，そこに詰める．

このヒューリスティクスは，最適なビン数の \(1.7\)倍以下のビン数しか使わないという保証をもつ．

first fit decreasing (FFD) ヒューリスティクス

アイテムのサイズを非減少順に並べ替えた後に， first fitヒューリスティクスを行う方法．

このヒューリスティクスは，最悪の場合の性能保証が \(11/9\) に漸近することが知られている．

first fitヒューリスティクス改良版として，best fitヒューリスティクスがある．

best fit (BF) ヒューリスティクス

アイテムを \(1,2,\ldots,n\) の順にビンに詰めていく． このとき，アイテムは詰め込み可能な最大の高さをもつビン（同点の場合には最小添字のビン）に詰めるものとする． どのビンに入れてもサイズの上限 \(1\) を超えてしまうなら，新たなビンを用意し，そこに詰める．

また，アイテムのサイズを非減少順に並べ替えた後にBFヒューリスティクスを適用したものを，best fit decreasing (BFD) ヒューリスティクスとよぶ．

上と同じランダムな問題例に対して適用してみる．

def FF(s, B):
    """First Fit heuristics for the bin Packing Problem.
    Parameters:
        - s: list with item widths
        - b: bin capacity
    Returns a list of lists with bin compositions.
    """
    remain = [B]  # keep list of empty space per bin
    sol = [[]]  # a list ot items (i.e., sizes) on each used bin
    for item in s:
        for (j, free) in enumerate(remain):
            if free >= item:
                remain[j] -= item
                sol[j].append(item)
                break
        else:  # does not fit in any bin
            sol.append([item])
            remain.append(B - item)
    return sol

def FFD(s, B):
    """First Fit Decreasing heuristics for the bin Packing Problem.
    Parameters:
        - s: list with item widths
        - b: bin capacity
    Returns a list of lists with bin compositions.
    """
    remain = [B]  # keep list of empty space per bin
    sol = [[]]  # a list ot items (i.e., sizes) on each used bin
    for item in sorted(s, reverse=True):
        for (j, free) in enumerate(remain):
            if free >= item:
                remain[j] -= item
                sol[j].append(item)
                break
        else:  # does not fit in any bin
            sol.append([item])
            remain.append(B - item)
    return sol

def BF(s, B):
    """Best Fit heuristics for the bin Packing Problem.
    Parameters:
        - s: list with item widths
        - b: bin capacity
    Returns a list of lists with bin compositions.
    """
    remain = [B]  # keep list of empty space per bin
    sol = [[]]  # a list ot items (i.e., sizes) on each used bin

    for item in s:
        best_remain = np.inf
        best_j = -1
        for (j, free) in enumerate(remain):
            rem = free - item
            if rem == 0:
                best_remain = rem
                best_j = j
                break
            elif rem > 0 and rem < best_remain:
                best_remain = rem
                best_j = j
        if best_j >= 0:
            remain[best_j] -= item
            sol[best_j].append(item)
        else:  # does not fit in any bin
            sol.append([item])
            remain.append(B - item)
    return sol

bins = FF(s, B)
print("FF", len(bins), "bins:")
print(bins)

bins = FFD(s, B)
print("FFD", len(bins), "bins:")
print(bins)

bins = BF(s, B)
print("BF", len(bins), "bins:")
print(bins)

カッティングストック問題

ビンパッキング（箱詰め）問題に類似の古典的な問題としてカッティングストック問題（cutting stock problem; 切断問題）がある．

\(m\)個の個別の幅をもった注文を，幅 \(b\) の原紙から切り出すことを考える． 注文 \(i=1,2,\ldots,m\) の幅 \(0 \leq w_i \leq b\) と注文数 \(q_i\) が与えられたとき， 必要な原紙の数を最小にするような切り出し方を求めよ．

ここでは，切断問題に対する Gilmore-Gomoryによる列生成法(column generation method)を構築する．

幅\(b\) の原紙からの \(m\)種類の注文の \(k\)番目の切断パターンを， \((t_1^k, t_2^k,\ldots,t_m^k)\) とする． ここで，\(t_i^k\) は注文 \(i\) が \(k\)番目の切断パターン \(k\) で切り出される数を表す． また，実行可能な切断パターン（箱詰め問題における詰め合わせ）とは，以下の式を満たすベクトル \((t_1^k, t_2^k,\ldots,t_m^k)\) を指す． \[ \sum_{i=1}^m t_i^k \leq b \]

実行可能な切断パターンの総数を \(K\) とする． 切断問題はすべての可能な切断パターンから， 注文 \(i\) を注文数 \(q_j\) 以上切り出し， かつ使用した原紙数を最小にするように 切断パターンを選択する問題になる． 切断パターン \(k\) を採用する回数を表す整数変数 \(x_k\) を 用いると，切断問題は以下の整数最適化問題として書くことができる．

\[ \begin{array}{l l l} minimize & \displaystyle\sum_{k=1}^K x_k & \\ s.t. & \displaystyle\sum_{k=1}^K t_i^k x_k \geq q_i & \forall i=1,2,\ldots,m \\ & x_k \in \mathbf{Z}_+ & \forall k=1,2,\ldots,K \end{array} \]

これを主問題(master problem)と呼ぶ． 主問題の線形最適化緩和問題を考え，その最適双対変数ベクトルを \(\lambda\) とする． このとき，被約費用が負の列（実行可能な切断パターン）を求める問題は，以下の整数ナップサック問題になる． \[ \begin{array}{l l l} maximize & \displaystyle\sum_{i=1}^m \lambda_i y_i & \\ s.t. & \displaystyle\sum_{i=1}^m w_i y_i \leq b & \\ & y_i \in \mathbf{Z}_+ & \forall i=1,2,\ldots,m \end{array} \]

EPS = 1.0e-6
def solveCuttingStock(w, q, b):
    """solveCuttingStock: use column generation (Gilmore-Gomory approach).
    Parameters:
        - w: list of item's widths
        - q: number of items of a width
        - b: bin/roll capacity
    Returns a solution: list of lists, each of which with the cuts of a roll.
    """
    t = []  # patterns
    m = len(w)
    # generate initial patterns with one size for each item width
    for (i, width) in enumerate(w):
        pat = [0] * m  # vector of number of orders to be packed into one roll (bin)
        pat[i] = int(b / width)
        t.append(pat)

    K = len(t)
    master = Model("master LP")  # master LP problem
    x = {}
    for k in range(K):
        x[k] = master.addVar(vtype="I", name="x(%s)" % k)
    master.update()

    orders = {}
    for i in range(m):
        orders[i] = master.addConstr(
            quicksum(t[k][i] * x[k] for k in range(K) if t[k][i] > 0) >= q[i],
            "Order(%s)" % i,
        )

    master.setObjective(quicksum(x[k] for k in range(K)), GRB.MINIMIZE)

    master.update()  # must update before calling relax()
    master.Params.OutputFlag = 0  # silent mode

    while True:
        relax = master.relax()
        relax.optimize()
        pi = [c.Pi for c in relax.getConstrs()]  # keep dual variables

        knapsack = Model("KP")  # knapsack sub-problem
        knapsack.ModelSense = -1  # maximize
        y = {}
        for i in range(m):
            y[i] = knapsack.addVar(lb=0, ub=q[i], vtype="I", name="y(%s)" % i)
        knapsack.update()

        knapsack.addConstr(quicksum(w[i] * y[i] for i in range(m)) <= b, "Width")
        knapsack.setObjective(quicksum(pi[i] * y[i] for i in range(m)), GRB.MAXIMIZE)

        knapsack.Params.OutputFlag = 0  # silent mode
        knapsack.optimize()

        if knapsack.ObjVal < 1 + EPS:  # break if no more columns
            break

        pat = [int(y[i].X + 0.5) for i in y]  # new pattern
        t.append(pat)

        # add new column to the master problem
        col = Column()
        for i in range(m):
            if t[K][i] > 0:
                col.addTerms(t[K][i], orders[i])
        x[K] = master.addVar(obj=1, vtype="I", name="x(%s)" % K, column=col)
        master.update()  # must update before calling relax()
        K += 1

    master.optimize()

    rolls = []
    for k in x:
        for j in range(int(x[k].X + 0.5)):
            rolls.append(
                sorted([w[i] for i in range(m) if t[k][i] > 0 for j in range(t[k][i])])
            )
    rolls.sort()
    return rolls


def CuttingStockExample():
    """CuttingStockExample: create toy instance for the cutting stock problem."""
    b = 110  # roll width (bin size)
    w = [20, 45, 50, 55, 75]  # width (size) of orders (items)
    q = [48, 35, 24, 10, 8]  # quantitiy of orders
    return w, q, b

def mkCuttingStock(s):
    """mkCuttingStock: convert a bin packing instance into cutting stock format"""
    w, q = [], []  # list of different widths (sizes) of items, their quantities
    for item in sorted(s):
        if w == [] or item != w[-1]:
            w.append(item)
            q.append(1)
        else:
            q[-1] += 1
    return w, q

def mkbinPacking(w, q):
    """mkbinPacking: convert a cutting stock instance into bin packing format"""
    s = []
    for j in range(len(w)):
        for i in range(q[j]):
            s.append(w[j])
    return s

w, q, b = CuttingStockExample()
rolls = solveCuttingStock(w, q, b)
print(len(rolls), "rolls:")
print(rolls)

同じ問題例をビンパッキング問題に変換し，first fit decreasingヒューリスティクスで解いてみる．

s = mkbinPacking(w, q)
bins = FFD(s,b)
print("FFD", len(bins), "bins:")

変動サイズベクトルパッキング問題

ビンのサイズが一定でなく，次元ごとの上限が異なる問題は，変動サイズベクトルパッキング問題(variable size vector packing problem)とよばれ， データセンターにおけるプロセスのマシンへの割当や，トラックへの荷物の割当に応用をもつ．

定式化

ビン \(j\) を使用したときの費用 \(c_j\) を導入し，最適化問題として定式化する．

アイテム \(i\) をビン \(j\) に詰めるとき \(1\) になる変数 \(x_{ij}\) と， ビン \(j\) の使用の可否を表す変数 \(y_j\) を用いることによって，変動サイズベクトルパッキング問題は，以下の整数最適化問題として記述できる．

\[ \begin{array}{l l l} minimize & \sum_{j=1}^U c_j y_{j} & \\ s.t. & \sum_{j=1}^{U} x_{ij} =1 & \forall i =1,2,\ldots,n \\ & \sum_{i=1}^n w_i^k x_{ij} \leq B_j^k y_j & \forall j=1,2,\ldots,U, k=1,2,\ldots,d \\ & x_{ij} \leq y_j & \forall i=1,2,\ldots,n, j=1,2,\ldots,U \\ & x_{ij} \in \{ 0,1 \} & \forall i=1,2,\ldots,n, j=1,2,\ldots,U \\ & y_j \in \{ 0,1 \} & \forall j=1,2,\ldots,U \end{array} \]

以下のプログラムでは，実行不能性を避けるために，割り当てできないアイテムの数を最小化し，その中で使用したビンの費用を最小化するように変更している．

def vbpp(s, B, c, penalty=1000.0):
    n = len(s)
    U = len(c)
    model = Model("bpp")
    # setParam("MIPFocus",1)
    x, y, z = {}, {}, {}
    for i in range(n):
        z[i] = model.addVar(vtype="B", name=f"z({i})")
        for j in range(U):
            x[i, j] = model.addVar(vtype="B", name=f"x({i},{j})")
    for j in range(U):
        y[j] = model.addVar(vtype="B", name=f"y({j})")
    model.update()
    # assignment constraints
    for i in range(n):
        model.addConstr(quicksum(x[i, j] for j in range(U)) + z[i] == 1, f"Assign({i})")
    # tighten assignment constraints
    for j in range(U):
        for i in range(n):
            model.addConstr(x[i, j] <= y[j], f"Strong({i},{j})")
    # bin capacity constraints
    for j in range(U):
        for k in range(d):
            model.addConstr(
                quicksum(s[i, k] * x[i, j] for i in range(n)) <= B[j, k] * y[j],
                f"Capac({j},{k})",
            )
    model.setObjective(
        quicksum(penalty * s[i, k] * z[i] for i in range(n) for k in range(d))
        + quicksum(c[j] * y[j] for j in range(U)),
        GRB.MINIMIZE,
    )

    model.update()
    model.__data = x, y, z
    return model

np.random.seed(123)
n = 5   #アイテム数
U = 2    #ビンの数
d = 2    #次元数
lb = 3   #アイテムサイズの下限
ub = 19   #アイテムサイズの上限
blb = 15 #ビン容量の下限
bub = 20 #ビン容量の上限
penalty = 1000
s = np.random.randint(lb, ub, (n, d))
B = np.random.randint(blb, bub, (U, d))
c = np.random.randint(500, 1000, U)

model = vbpp(s, B, c, penalty)
model.optimize()

x, y, z = model.__data
bins = [[] for j in range(U)]
for (i, j) in x:
    if x[i, j].X > 0.5:
        bins[j].append(s[i])
unassigned = []
for i in z:
    if z[i].X > 0.5:
        unassigned.append(s[i])

print("unassigned items=", unassigned)
for j in range(U):
    weight = np.zeros(d, dtype=int)
    for s in bins[j]:
        weight += s
    print(weight, "<=", B[j])

Gurobi Optimizer version 11.0.1 build v11.0.1rc0 (mac64[arm] - Darwin 23.0.0 23A344)

CPU model: Apple M2 Max
Thread count: 12 physical cores, 12 logical processors, using up to 12 threads

Optimize a model with 19 rows, 17 columns and 59 nonzeros
Model fingerprint: 0x30b9b293
Variable types: 0 continuous, 17 integer (17 binary)
Coefficient statistics:
  Matrix range     [1e+00, 2e+01]
  Objective range  [6e+02, 3e+04]
  Bounds range     [1e+00, 1e+00]
  RHS range        [1e+00, 1e+00]
Found heuristic solution: objective 107000.00000
Presolve removed 19 rows and 17 columns
Presolve time: 0.00s
Presolve: All rows and columns removed

Explored 0 nodes (0 simplex iterations) in 0.01 seconds (0.00 work units)
Thread count was 1 (of 12 available processors)

Solution count 2: 55239 107000 

Optimal solution found (tolerance 1.00e-04)
Best objective 5.523900000000e+04, best bound 5.523900000000e+04, gap 0.0000%
unassigned items= [array([17,  5]), array([9, 4]), array([ 6, 13])]
[17 16] <= [18 16]
[15  5] <= [16 15]

ヒューリスティクス

変動サイズベクトルパッキング問題に対しては，BFD (Best Fit Decreasing) の変形が提案されている．

与えられた \(U\) 本のビンに容量を超えない割当が可能かどうかを判定する問題を考える．

残りのアイテムの集合を \(NR\)，残りのビンの集合を \(BR\) とする．

次元\(k\)のアイテムのサイズの和を \(W_k\) とする． \[ W_k = \sum_{i \in NR} w_i^k \]

ビン \(j\) の \(k\) 次元の残り容量を \(r_j^k\) とする． 次元\(k\)の残り容量の和を \(R_k\) とする．

\[ R_k = \sum_{j \in BR} r_j^k \]

\(W_k\) が大きい次元ほど， \(R_k\) が小さい次元ほど重要度が高いと考え， 次元 \(k\) の重要度（次元を統合したサイズ）を表す \(s_k\) を導入する． \(s_k\) は以下の3通りが考えられる．

\[ \begin{array}{lll} s_k & = & W_k \\ s_k & = & 1/R_k \\ s_k &= & W_k/R_k \end{array} \]

より一般的に，パラメータ \(0 \leq \alpha \leq 2\) を用いて以下のように定義する． \[ s_k = (W_k)^{\alpha} / (R_k)^{2-\alpha} \]

この尺度を用いてビン \(j\) の（次元を統合した）サイズ \[ \sum_{k=1}^d s_k r_j^k \ \ \ \forall j \in BR \] とアイテム \(i\) のサイズ \[ \sum_{k=1}^d s_k w_i^k \ \ \ \forall i \in NR \] を計算する．

この新たなサイズを用いてBFDヒューリスティクスは，アイテム中心型とビン中心型の2通りが定義できる．

アイテム中心型BFDヒューリスティクス

未割当アイテムに対して以下の操作を繰り返す． 1. （次元を統合した）サイズを計算 2. サイズが最大のアイテムを，割り当て可能な最小の残りサイズをもつビンに割り当てる．（どのビンにも割り当て不能なら終了；実行不能）

実際問題を解く際には，単に実行不能を返すだけでなく，なるべく多くのアイテムを詰めたいので，終了判定条件を変える必要がある．

ビン中心型BFDヒューリスティクス

残りのビン集合に対して以下の操作を繰り返す． 終了後に未割当のアイテムが残っていたら実行不能と判定する．

1. （次元を統合した）ビンのサイズを計算
2. サイズが最小のビンを選択
3. 未割当のアイテムに対してサイズを計算し，ビンに入る最大のものを入れる

以下に，アイテム中心型のBFDヒューリスティクスのコードと，上と同じ例題を解いたときの結果を示す．

def BFD(w, B, alpha=1.0):
    """
    Item Based Best Fit Decreasing (BFD) heuristcs for variable size vector packing problem
    """
    n, d = w.shape
    U = len(B)
    unassigned = []
    bins = [[] for j in range(U)]
    while len(w) > 0:
        W = w.sum(axis=0)  # アイテムに対する残りサイズの和（次元ごと）
        R = B.sum(axis=0)  # ビンに対する残り容量の和（次元ごと）
        s = (W**alpha) / (R + 0.000001) ** (2.0 - alpha)  # 次元の重要度
        BS = B @ s  # ビンのサイズ（次元統合後）
        WS = w @ s  # アイテムのサイズ（次元統合後）
        max_item_idx = np.argmax(WS)
        max_item = w[max_item_idx]
        for j in np.argsort(BS):
            remain = B[j] - max_item
            if np.all(remain >= 0):  # 詰め込み可能
                B[j] = remain
                bins[j].append(max_item)
                break
        else:
            unassigned.append(max_item)
        w = np.delete(w, max_item_idx, axis=0)
    return bins, unassigned

np.random.seed(123)
n = 10
U = 4
d = 2
lb = 3
ub = 10
blb = 15
bub = 20
w = np.random.randint(lb, ub, (n, d))
B = np.random.randint(blb, bub, (U, d))
BB = B.copy()
bins, unassigned = BFD(w, B, alpha=0.0)
print("Unassigned=", unassigned)
for j in range(U):
    weight = np.zeros(d, dtype=int)
    for s in bins[j]:
        weight += s
    print(weight, "<=", BB[j])

Unassigned= [array([3, 9]), array([3, 4])]
[15  9] <= [15 16]
[16  9] <= [18 19]
[15 13] <= [15 15]
[18  8] <= [19 16]

2次元パッキング問題

2次元（長方形；矩形）パッキング問題 (2-dimensioanl rectangular packing problem) は, 与えられた長方形を平面上に互いに重ならないように配置する問題である．

メタヒューリスティクスによる最適化

長方形集合 \(I=\{1,2,\ldots,n\}\) の各要素 \(i\) は \(m_i\) 種類のモードが与えられ, 各モード \(k(=1,2,\ldots,m_i)\) に対して， 幅 \(w_i^k\) と高さ \(h_i^k\) が与えられる. 配置は, 各長方形について 1つのモードを選択し, さらに左下の位置の座標値 \((x,y)\) を与えることで定まる. 配置はなるべくコンパクトなことが望ましく， \(x,y\) 軸に対して区分的線形な費用関数を定義することによって，目的関数を定める．

以下の論文に掲載されているメタヒューリスティクスをPythonから呼び出して使用する．

S. Imahori, M. Yagiura, T. Ibaraki:
Improved Local Search Algorithms for the Rectangle Packing Problem
with General Spatial Costs,
European Journal of Operational Research 167 (2005) 48-67

このコード自体はオープンソースではないが，研究用なら筆者に連絡すれば利用可能である．詳細は，付録1の OptPack を参照されたい．

以下，ソルバーのパラメータである．

• size: 長方形数
• pft_x: x軸コストがあるかどうか (0: あり, 1: なし (w,xgr を利用))
• pft_y: y軸コストがあるかどうか (0: あり, 1: なし (h,ygr を利用))
• w: 入れ物の幅 (x 軸コストが無い場合のみ意味あり)
• xgr: x方向, 入れ物からの超過に対するペナルティの傾き
• h: 入れ物の高さ (y 軸コストが無い場合のみ意味あり)
• ygr: y方向, 入れ物からの超過に対するペナルティの傾き
• time: 最大計算時間 (終了条件)
• move: 最大移動回数 (終了条件)
• lopt: 最大局所最適解数 (終了条件, MLS, ILSのみ意味あり)
• h_type: メタ戦略の種類 (1: 局所探索, 2: MLS(多スタート局所探索）, 3: ILS（反復局所探索）)

ベンチマーク問題例を解き，結果をPlotlyで図示してみる．

---

source

packingmeta

 packingmeta (fn, timelimit)

#fn = folder + "ami33"
fn = folder + "rp200"
timelimit = 1
opt_val, xy, wh, fig = packingmeta(fn, timelimit)
print("obj = ", opt_val)
plotly.offline.plot(fig);

出力（一部）

BEST obj = 1313.0, ( pmax = 675.0, qmax = 638.0 ) time = 1.00 eva1 = 1313.0, eva2 = 89.0, eva3 = 135834.0

ID: (x, y) [w, h] 1: (175.0, 164.0) [68.0, 45.0] 2: (267.0, 328.0) [45.0, 70.0] 3: (538.0, 0.0) [66.0, 42.0] 4: (563.0, 120.0) [68.0, 42.0] 5: (210.0, 269.0) [65.0, 42.0] … 198: (102.0, 164.0) [58.0, 75.0] 199: (222.0, 552.0) [55.0, 80.0] 200: (277.0, 503.0) [85.0, 55.0]

Image("../figure/2dbpp.png")

cp-sat による最適化

２D Packing モデルクラス

---

source

plot_solution

 plot_solution (ax, instance:__main__.Instance,
                solution:__main__.Solution)

---

source

RectangleKnapsackWithRotationsModel

 RectangleKnapsackWithRotationsModel (instance:__main__.Instance)

Initialize self. See help(type(self)) for accurate signature.

---

source

Solution

 Solution (placements:list[typing.Optional[__main__.Placement]])

*Usage docs: https://docs.pydantic.dev/2.8/concepts/models/

A base class for creating Pydantic models.

Attributes: class_vars: The names of classvars defined on the model. private_attributes: Metadata about the private attributes of the model. signature: The signature for instantiating the model.

__pydantic_complete__: Whether model building is completed, or if there are still undefined fields.
__pydantic_core_schema__: The pydantic-core schema used to build the SchemaValidator and SchemaSerializer.
__pydantic_custom_init__: Whether the model has a custom `__init__` function.
__pydantic_decorators__: Metadata containing the decorators defined on the model.
    This replaces `Model.__validators__` and `Model.__root_validators__` from Pydantic V1.
__pydantic_generic_metadata__: Metadata for generic models; contains data used for a similar purpose to
    __args__, __origin__, __parameters__ in typing-module generics. May eventually be replaced by these.
__pydantic_parent_namespace__: Parent namespace of the model, used for automatic rebuilding of models.
__pydantic_post_init__: The name of the post-init method for the model, if defined.
__pydantic_root_model__: Whether the model is a `RootModel`.
__pydantic_serializer__: The pydantic-core SchemaSerializer used to dump instances of the model.
__pydantic_validator__: The pydantic-core SchemaValidator used to validate instances of the model.

__pydantic_extra__: An instance attribute with the values of extra fields from validation when
    `model_config['extra'] == 'allow'`.
__pydantic_fields_set__: An instance attribute with the names of fields explicitly set.
__pydantic_private__: Instance attribute with the values of private attributes set on the model instance.*

---

source

Placement

 Placement (x:int, y:int, rotated:bool=False)

*Usage docs: https://docs.pydantic.dev/2.8/concepts/models/

A base class for creating Pydantic models.

Attributes: class_vars: The names of classvars defined on the model. private_attributes: Metadata about the private attributes of the model. signature: The signature for instantiating the model.

__pydantic_complete__: Whether model building is completed, or if there are still undefined fields.
__pydantic_core_schema__: The pydantic-core schema used to build the SchemaValidator and SchemaSerializer.
__pydantic_custom_init__: Whether the model has a custom `__init__` function.
__pydantic_decorators__: Metadata containing the decorators defined on the model.
    This replaces `Model.__validators__` and `Model.__root_validators__` from Pydantic V1.
__pydantic_generic_metadata__: Metadata for generic models; contains data used for a similar purpose to
    __args__, __origin__, __parameters__ in typing-module generics. May eventually be replaced by these.
__pydantic_parent_namespace__: Parent namespace of the model, used for automatic rebuilding of models.
__pydantic_post_init__: The name of the post-init method for the model, if defined.
__pydantic_root_model__: Whether the model is a `RootModel`.
__pydantic_serializer__: The pydantic-core SchemaSerializer used to dump instances of the model.
__pydantic_validator__: The pydantic-core SchemaValidator used to validate instances of the model.

__pydantic_extra__: An instance attribute with the values of extra fields from validation when
    `model_config['extra'] == 'allow'`.
__pydantic_fields_set__: An instance attribute with the names of fields explicitly set.
__pydantic_private__: Instance attribute with the values of private attributes set on the model instance.*

---

source

Instance

 Instance (container:__main__.Container,
           rectangles:list[__main__.Rectangle])

*Usage docs: https://docs.pydantic.dev/2.8/concepts/models/

A base class for creating Pydantic models.

Attributes: class_vars: The names of classvars defined on the model. private_attributes: Metadata about the private attributes of the model. signature: The signature for instantiating the model.

__pydantic_complete__: Whether model building is completed, or if there are still undefined fields.
__pydantic_core_schema__: The pydantic-core schema used to build the SchemaValidator and SchemaSerializer.
__pydantic_custom_init__: Whether the model has a custom `__init__` function.
__pydantic_decorators__: Metadata containing the decorators defined on the model.
    This replaces `Model.__validators__` and `Model.__root_validators__` from Pydantic V1.
__pydantic_generic_metadata__: Metadata for generic models; contains data used for a similar purpose to
    __args__, __origin__, __parameters__ in typing-module generics. May eventually be replaced by these.
__pydantic_parent_namespace__: Parent namespace of the model, used for automatic rebuilding of models.
__pydantic_post_init__: The name of the post-init method for the model, if defined.
__pydantic_root_model__: Whether the model is a `RootModel`.
__pydantic_serializer__: The pydantic-core SchemaSerializer used to dump instances of the model.
__pydantic_validator__: The pydantic-core SchemaValidator used to validate instances of the model.

__pydantic_extra__: An instance attribute with the values of extra fields from validation when
    `model_config['extra'] == 'allow'`.
__pydantic_fields_set__: An instance attribute with the names of fields explicitly set.
__pydantic_private__: Instance attribute with the values of private attributes set on the model instance.*

---

source

Container

 Container (width:int, height:int)

*Usage docs: https://docs.pydantic.dev/2.8/concepts/models/

A base class for creating Pydantic models.

Attributes: class_vars: The names of classvars defined on the model. private_attributes: Metadata about the private attributes of the model. signature: The signature for instantiating the model.

__pydantic_complete__: Whether model building is completed, or if there are still undefined fields.
__pydantic_core_schema__: The pydantic-core schema used to build the SchemaValidator and SchemaSerializer.
__pydantic_custom_init__: Whether the model has a custom `__init__` function.
__pydantic_decorators__: Metadata containing the decorators defined on the model.
    This replaces `Model.__validators__` and `Model.__root_validators__` from Pydantic V1.
__pydantic_generic_metadata__: Metadata for generic models; contains data used for a similar purpose to
    __args__, __origin__, __parameters__ in typing-module generics. May eventually be replaced by these.
__pydantic_parent_namespace__: Parent namespace of the model, used for automatic rebuilding of models.
__pydantic_post_init__: The name of the post-init method for the model, if defined.
__pydantic_root_model__: Whether the model is a `RootModel`.
__pydantic_serializer__: The pydantic-core SchemaSerializer used to dump instances of the model.
__pydantic_validator__: The pydantic-core SchemaValidator used to validate instances of the model.

__pydantic_extra__: An instance attribute with the values of extra fields from validation when
    `model_config['extra'] == 'allow'`.
__pydantic_fields_set__: An instance attribute with the names of fields explicitly set.
__pydantic_private__: Instance attribute with the values of private attributes set on the model instance.*

---

source

Rectangle

 Rectangle (width:int, height:int, value:int=1)

*Usage docs: https://docs.pydantic.dev/2.8/concepts/models/

A base class for creating Pydantic models.

Attributes: class_vars: The names of classvars defined on the model. private_attributes: Metadata about the private attributes of the model. signature: The signature for instantiating the model.

__pydantic_complete__: Whether model building is completed, or if there are still undefined fields.
__pydantic_core_schema__: The pydantic-core schema used to build the SchemaValidator and SchemaSerializer.
__pydantic_custom_init__: Whether the model has a custom `__init__` function.
__pydantic_decorators__: Metadata containing the decorators defined on the model.
    This replaces `Model.__validators__` and `Model.__root_validators__` from Pydantic V1.
__pydantic_generic_metadata__: Metadata for generic models; contains data used for a similar purpose to
    __args__, __origin__, __parameters__ in typing-module generics. May eventually be replaced by these.
__pydantic_parent_namespace__: Parent namespace of the model, used for automatic rebuilding of models.
__pydantic_post_init__: The name of the post-init method for the model, if defined.
__pydantic_root_model__: Whether the model is a `RootModel`.
__pydantic_serializer__: The pydantic-core SchemaSerializer used to dump instances of the model.
__pydantic_validator__: The pydantic-core SchemaValidator used to validate instances of the model.

__pydantic_extra__: An instance attribute with the values of extra fields from validation when
    `model_config['extra'] == 'allow'`.
__pydantic_fields_set__: An instance attribute with the names of fields explicitly set.
__pydantic_private__: Instance attribute with the values of private attributes set on the model instance.*

fig, ax = plt.subplots(1)
plot_solution(ax, instance, model.solution)
plt.show()

確率的ビンパッキング問題

ビンパッキング問題の応用では，アイテムに関する情報が事前に与えられておらず， 時間の経過とともにアイテムが到着し，そのサイズが判明するケースがある． このような仮定を課した問題をオンラインビンパッキング問題（online bin packing problem）とよぶ． より正確に言うと，先に到着するアイテムについての情報をもたずに入れるビンを決め，さらに1度入れたら，それを移動することができないという制限を課した問題である．

この条件を少し緩和した問題として，到着するアイテムのサイズの確率分布が既知であると仮定した問題が考えられる． このような問題を確率的ビンパッキング問題(stochastic bin packing problem)とよぶ．

実務における問題は，オンラインと確率的の両者の性質をもつ．以下の分類にしたがって，解法を選ぶことが望ましい．

• 分布が未知で未来の情報がまったくない： オンラインビンパッキングのヒューリスティクスを用いる．
• 分布が定常：確率分布を推定し，確率的ビンパッキング問題として解く．
• 分布が非定常：オンラインで到着するアイテムのサイズから適応的に分布を推定し，ヒューリスティクスを適用する．
• 近い未来の情報は既知で，ある時点以降の情報がまったくない： 既知の部分までを確定的な問題として求解して，ローリングホライズン方式を適用する．
• 近い未来の情報は既知で，先に行くにしたがって情報が不確実になる： 既知の部分までを確定的な問題として求解して，ローリングホライズン方式を適用するが，確率的な情報を利用して，最終状態を適正なものに近づける．

アイテムのサイズを離散値とした確率的ビンパッキング問題を考える．

ビンの容量を正の整数 \(B\) とする． アイテムのサイズの種類は \(J\) 種類であり，小さい方から順に \(s_1 <s_2 <\cdots <s_J\) とする． ここで \(s_j (j=1,\ldots,J)\) は \([1,B]\) の整数とする．サイズが \(s_j\) のアイテムが発生する確率を \(p_j\) とする． ここでは，\(p_j\) は非負の有理数であり，その和は \(1\) であるとする． \(i\)番目のアイテムは，独立に確率 \(p_j\) でサイズ \(s_j\) になるように決められる．

サイズ \(s_j\) のアイテムを高さ \(h\) のビンに割り当てる確率を表す実数変数 \(x_{jh}\) を用いた以下の線形最適化問題を導入する．

\[ \begin{array}{l l l } minimize & \sum_{h=1}^{B-1} (B-h) \cdot \left(\sum_{j=1}^J x_{j,h-s_j}-\sum_{j=1}^J x_{jh} \right) & \\ s.t. & \sum_{h=0}^{B-1} x_{jh} =p_j & j=1,\ldots,J \\ & \sum_{j=1}^J x_{jh} \leq \sum_{j=1}^J x_{j,h-s_j} & h=1,\ldots,B-1 \\ & x_{jh} \geq 0 & j=1,\ldots,J,h=0,\ldots,B-1 \\ & x_{jh} =0 & j=1,\ldots,J,h=B-s_j,\ldots,B-1 \end{array} \]

目的関数は，ビンの余り \((B-h)\) が超過して生成される確率に，余り \((B-h)\) を乗じたものを最小化することを表す． 最初の制約は，サイズ \(s_j\) のアイテムがいずれかのビンに割り当てられる確率が，発生する確率に等しいことを表す． 2番目の制約は，高さ \(h\) が生成される確率（式の右辺）が，消滅する確率（式の左辺）以上であることを表す．（この式の余裕変数(slack)が，\(h\) が余分になる確率を表す．） 3番目の制約は，変数の非負条件であり，最後の制約はビンのサイズの上限を超えてしまう割当確率が \(0\) であることを表す．

容量 \(10\) のビンに対して，\(1,2,5,6\) のサイズをもつアイテムが，それぞれ確率 \(0.3,0.3,0.3,0.1\) で発生する場合を考える． アイテムのサイズはSciPyのrv_discrete関数を用いて生成する．

B = 10
s = [1, 2, 5, 6]
p = [0.3, 0.3, 0.3, 0.1]
J = len(s)
size = rv_discrete(values=(s, p))

model = Model("discrete bpp")
x, slack = {}, {}
for j in range(J):
    for h in range(B):
        if s[j] + h <= B:
            x[j, h] = model.addVar(vtype="C", name=f"x[{j},{h}]")
for h in range(1, B):
    slack[h] = model.addVar(vtype="C", name=f"slack[{h}]")
model.update()
for j in range(J):
    model.addConstr(quicksum(x[j, h] for h in range(B) if (j, h) in x) == p[j])
for h in range(1, B):
    model.addConstr(
        quicksum(x[j, h] for j in range(J) if (j, h) in x) + slack[h]
        == quicksum(x[j, h - s[j]] for j in range(J) if h >= s[j])
    )
model.setObjective(quicksum((B - h) * slack[h] for h in range(1, B)), GRB.MINIMIZE)

model.optimize()

目的関数 \(0\) が残り容量の期待値であり，漸近的に完全パッキングが可能であることを表す． また，どのアイテムをどの高さのビンに（確率的に）割り振るかの指針を得ることができる．

for (j, h) in x:
    if x[j, h].X > 0.001:
        print("size", s[j], "item should assign the bin with height", h, " with prob. ", x[j, h].X)

オンラインビンパッキング問題

アイテムがすべて既知として問題を解く状況をオフラインとよび，将来発生するアイテムに対する情報をもたずに最適化する状況をオンラインとよぶ． 前に示した first fit decreasingヒューリスティクスはオフラインでしか使えず，first fitヒューリスティクスはオンラインで使えるが性能が悪い．

以下では， オンラインの環境下で，アイテムサイズが離散分布の場合に良い性能を示す2乗和法（sum-of-square algorithm）を紹介し， best fitヒューリスティクスと比較する．

途中までアイテムを詰めた状態を考える． このとき，アイテムのサイズは整数であるので，空でなくかつ一杯になっていないビンの高さ \(h\) は \([1,B-1]\) の整数である． 高さが \(h\) のビンの数を \(M(h)\) と書く． オンラインの環境下においては，将来どのようなサイズが到着するか分からないと仮定するが， このような状態における策として，\(M(h)\) がなるべく均等になるように準備しておく方法が考えられる． もし，すべての高さのビンがあれば，どのようなアイテムが到着しても，ビンの高さがちょうど \(B\) になるように詰めることができるからである．

均等にするための常套手段として，2乗和を最小化する方法がある． 和が一定の数の分割において，2乗和が最小になるのは，すべての数が同じ場合である． （たとえば，\(9\) の3分割では，\(3^2+3^2+3^2=27\) が最小で，\(9^2+0^2+0^2=81\) が最大になる．） 2乗和法では，このアイディアに基づいて，以下の指標を最小にするビンにアイテムを入れる． \[ SS= \sum_{h=1}^{B-1} M(h)^2 \]

サイズ \(w\) のアイテムを，高さ \(h\) のビンに入れたときの上の指標の増加量 \(\Delta\) は， \[ \Delta= \left\{ \begin{array}{l l } 2M(w)+1 & 新しいビンを生成するとき （h=0） \\ -M(B-w)+1 & ビンが一杯になったとき （h=B-w） \\ 2(M(h+w)-M(h))+2 & その他 \end{array} \right. \] と計算できる．したがって，指標 \(SS\) を最小にする高さを求めるには，2乗和を再計算する 必要はなく，\(\Delta\) を最小にする \(h\) を求めれば良い．これは \(O(B)\) 時間でできる．

容量 \(10\) のビンに対して，\(1,2,5,6\) のサイズをもつアイテムが，それぞれ確率 \(0.3,0.3,0.3,0.1\) で発生する場合を考える．

B = 10
s = [1, 2, 5, 6]
p = [0.3, 0.3, 0.3, 0.1]
J = len(s)
size = rv_discrete(values=(s, p))
n = 10000

M = np.zeros(B + 1, dtype=int)
w = size.rvs(size=n)
for i in range(n):
    min_delta = 2 * w[i] + 1
    best_h = 0
    for h in range(1, B + 1 - w[i]):
        if M[h] == 0:
            continue
        if h + w[i] == B:
            delta = -M[h] + 1
        else:
            delta = 2 * (M[h + w[i]] - M[h]) + 2
        if delta < min_delta:
            min_delta = delta
            best_h = h
    if best_h > 0:
        M[best_h] -= 1
    M[best_h + w[i]] += 1
print("M=", M)
print(sum(M), "bins")

best fitヒューリスティクスと比較してみよう．

sol = BF(w, B)
print(len(sol),"bins")