定式化
ビンの数の上限 $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} $$この定式化は，大規模な実際問題を解く際には使うべきではない．データが整数の場合には，枝フロー変数を用いたより強い定式化が提案されている． 詳細は， 以下のサイトを参照されたい．
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)
AMPLによる求解
AMPLモデル
param n; # number of items
param U; # maximum number of bins
param s {1..n};
param B;
var X {1..n, 1..U} binary;
var Y {1..U} binary;
minimize bins: sum {j in 1..U} Y[j];
subject to
Take {i in 1..n}: sum {j in 1..U} X[i,j] = 1;
Cap {j in 1..U}: sum {i in 1..n} s[i] * X[i,j] <= B * Y[j];
Activate {i in 1..n, j in 1..U}: X[i,j] <= Y[j];ampl = AMPL(Environment(AMPL_PATH))
ampl.setOption('solver', 'gurobi')
ampl.read("../ampl/bpp.mod")
ampl.param['n'] = n
ampl.param['U'] = U
ampl.param['s'] = s
ampl.param['B'] = B
# solve and report solution
ampl.solve()
bins = ampl.obj['bins']
print("Objective is:", bins.value())
X = ampl.var['X']
Y = ampl.var['Y']
for j in range(1,U+1):
v = Y[j].value()
if v > 0:
print("bin {}".format(j), end="\t")
for i in range(1,n+1):
v = X[i,j].value()
if v > 0:
print(s[i-1], end=" ")
print()
ヒューリスティクス
ビンの数の上界 $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:")
AMPLによる求解
AMPLモデル (csp_master.mod)
param K;
param m;
param q {1..m};
param t {1..K, 1..m};
var x {1..K} >= 0, integer;
minimize rolls: sum {k in 1..K} x[k];
subject to
Demand {i in 1..m}: sum {k in 1..K} t[k,i] * x[k] >= q[i];AMPLモデル (csp_knapsack.mod)
param m;
param B;
param w {1..m};
param lambda {1..m};
var y {1..m} >= 0, integer;
maximize z: sum {i in 1..m} lambda[i] * y[i];
subject to
Feasible: sum {i in 1..m} w[i] * y[i] <= B;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(pat)
# column generation loop
# initialize master problem
master = AMPL(Environment(AMPL_PATH))
master.option['solver'] = 'gurobi'
master.option['relax_integrality'] = 1
master.read("../ampl/csp_master.mod")
master.param['K'] = K
master.param['m'] = m
master.param['q'] = q
t_ = master.param['t']
for k in range(1, K + 1):
for i in range(1, m + 1):
t_[k,i] = t[k-1][i-1]
# initialize subproblem
kp = AMPL(Environment(AMPL_PATH))
kp.option['solver'] = 'gurobi'
kp.read("../ampl/csp_knapsack.mod")
kp.param['m'] = m
kp.param['B'] = B
kp.param['w'] = w
while True:
print("current patterns:")
for ti in t:
print(ti)
print()
master.solve()
rolls = master.obj['rolls']
print("master objective:", rolls.value())
demand_ = master.con['Demand']
lambda_ = {}
for i in range(1, m + 1):
lambda_[i] = demand_[i].dual()
print("lambda {} {}".format(i, lambda_[i]))
# setup knapsack subproblem
kp.param['lambda'] = lambda_
kp.solve()
z = kp.obj['z']
print("subproblem objective:", z.value())
if z.value() <= 1:
break
# update new pattern
pat = [0] * m # vector of number of orders to be packed into one roll (bin)
y = kp.var['y']
for i in range(1, m + 1):
v = int(round(y[i].value()))
print("y[{}] = {} : {}".format(i,v,y[i].value()))
pat[i-1] = int(v)
print("added pattern", pat)
t.append(pat)
K += 1
master.param['K'] = K
for i in range(1, m + 1):
t_[K, i] = pat[i-1]
master.option['relax_integrality'] = 0
master.solve()
rolls = master.obj['rolls']
print("master objective:", rolls.value())
x = master.var['x']
rolls = []
for k in range(1,K+1):
n = int(x[k].value() + .5) # number of times pattern is used
for j in range(n):
rolls.append(sorted([w[i] for i in range(m) if t[k-1][i] > 0 for j in range(t[k-1][i])]))
rolls.sort()
print(rolls)
定式化
ビンのサイズが一定でなく，次元ごとの上限が異なる問題は，変動サイズベクトルパッキング問題(variable size vector packing problem)とよばれ， データセンターにおけるプロセスのマシンへの割当や，トラックへの荷物の割当に応用をもつ． 以下では，ビン $j$ を使用したときの費用 $c_j$ を導入し，最適化問題として定式化する．
アイテム $i$ をビン $j$ に詰めるとき $1$ になる変数 $x_{ij}$ と， ビン $j$ の使用の可否を表す変数 $y_j$ を用いることによって，変動サイズベクトルパッキング問題は，以下の整数最適化問題として記述できる．
$$ \begin{array}{l l l} minimize & \displaystyle\sum_{j=1}^U c_j 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^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])
ヒューリスティクス
変動サイズベクトルパッキング問題に対しては，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ヒューリスティクス
未割当アイテムに対して以下の操作を繰り返す．
- （次元を統合した）サイズを計算
- サイズが最大のアイテムを，割り当て可能な最小の残りサイズをもつビンに割り当てる．（どのビンにも割り当て不能なら終了；実行不能）
実際問題を解く際には，単に実行不能を返すだけでなく，なるべく多くのアイテムを詰めたいので，終了判定条件を変える必要がある．
ビン中心型BFDヒューリスティクス
残りのビン集合に対して以下の操作を繰り返す． 終了後に未割当のアイテムが残っていたら実行不能と判定する．
- （次元を統合した）ビンのサイズを計算
- サイズが最小のビンを選択
- 未割当のアイテムに対してサイズを計算し，ビンに入る最大のものを入れる
以下に，アイテム中心型の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
BB = B.copy()
print("w=",w)
print("B=",B)
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])
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で図示してみる．
fn = folder + "rp200"
timelimit = 1
opt_val, xy, wh, fig = packingmeta(fn, timelimit)
print("obj = ", opt_val)
plotly.offline.plot(fig);
確率的ビンパッキング問題
ビンパッキング問題の応用では，アイテムに関する情報が事前に与えられておらず， 時間の経過とともにアイテムが到着し，そのサイズが判明するケースがある． このような仮定を課した問題をオンラインビンパッキング問題（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")