整数ナップサック問題に対する動的最適化
最初に考えるのは整数ナップサック問題(integer knapsack problem)とよばれる問題である．
$n$ 個のアイテムから成る有限集合 $N$， 各々の $i \in N$ の重量 $s_i \in \mathbf{Z}_+$ と価値 $v_i \in \mathbf{Z}_+$， およびナップサックの重量の上限 $b \in \mathbf{Z}_+$ が与えられたとき, アイテムの重量の合計が $b$ を超えないようなアイテムの詰め合わせの中で， 価値の合計が最大のものを求めよ．ただし，各アイテムは何個でもナップサックに詰めて良いものと仮定する．
この問題は，非負の整数の値をとる変数 $x_i \in \{0,1,2,\cdots\}=\mathbf{Z}_+$ を用いて，以下のように定式化できる． $$ \begin{array}{l l l} maximize & \sum_{i \in N} v_i x_i & \\ s.t. & \sum_{i \in N} s_i x_i \leq b & \\ & x_i \in \mathbf{Z}_+ & i \in N \end{array} $$
動的最適化の特徴は，自明な部分問題からはじめて順次大きな部分問題を解いていく点にある． 整数ナップサック問題の場合には，ナップサックの重量上限を制限した部分問題を考えれば良い．
ナップサックに入っているアイテムの重量が $\theta ~(=0,1,2,\cdots,b)$ 以下のときの最大価値を $f(\theta)$ と書くことにしよう． もとの問題の最適値は，重量が $b$ 以下の中で最大の価値であるから $f(b)$ である．
$f(0)=0$ は自明であり，さらに $f(b)=-\infty ~(b<0)$ と定義しておく． $f(\theta)$ は，$\theta$ から $s_i$ だけ軽い重量が詰まっているナップサックにおける最大価値 $f(\theta-s_i)$ に $v_i$ の価値を加えることによって得られる最大の価値であるので， $$ f(\theta) = \max_{i \in N} \left\{ f(\theta - s_i)+v_i, 0 \right\} \ \ \ \theta =1,2,\ldots,b $$ の関係式が得られる．これが再帰方程式になる．
再帰方程式を辞書によるメモ化によって高速化したものが，動的最適化アルゴリズムである． 辞書に，ナップサックに入れたアイテムの情報を保管することによって，最適解も得ることができる．
このアルゴリズムは，各 $\theta (=0,1,\ldots,b)$ に対してアイテムの数 $n$ だけの計算時間がかかるので， 全体として $O(n b)$ 時間かかる．
この計算時間は，一見すると入力サイズの多項式関数に見えるが， 実はナップサックの重量の上限 $b$ の正確な入力サイズは，$b$ でなく $\lceil \log_2 b \rceil$ なのである． ここで $\lceil \cdot \rceil$ は天井関数(ceiling function)であり，$\cdot$ 以上の最小の整数（切り上げ）を意味する． $\lceil \log_2 b \rceil$ を $\beta$ とおくと，動的計画アルゴリズムの計算時間は $O(n 2^{\beta})$ となる． これは，正確な入力サイズ $\beta$ 多項式関数ではない． ちなみに， このようなアルゴリズムを擬多項式時間(pseudo-polynomial time)アルゴリズムとよぶ．
再帰方程式をコードにすると以下のように書ける．
def ikpdp(s, v, b):
def f(b):
if b == 0:
return 0, -1
if b < 0:
return -999999, -1
if b in memo:
return memo[b]
else:
max_value = 0
prev = -1
for i, size in enumerate(s):
if f(b - size)[0] + v[i] > max_value:
max_value = f(b - size)[0] + v[i]
prev = i
memo[b] = max_value, prev
return memo[b]
memo = {}
opt_val, prev = f(b)
x = [0 for i in range(len(s))]
while True:
val, prev = memo[b]
x[prev] += 1
b -= s[prev]
if b <= 0:
break
return opt_val, x
s = [2, 3, 4, 5]
v = [16, 19, 23, 28]
b = 7
opt_val, x = ikpdp(s, v, b)
print("Opt. value=", opt_val)
print("Sol.=", x)
0-1ナップサック問題に対する動的最適化
この問題は上で考えた整数ナップサック問題とほぼ同じであるが，各アイテムが $1$個ずつしかない点が異なる．，
$n$個のアイテムから成る有限集合 $N$， 各々の $i \in N$ の重量 $s_i \in \mathbf{Z}_+$ と価値 $v_i \in \mathbf{Z}_+$， およびナップサックの重量の上限 $b \in \mathbf{Z}_+$ が与えられたとき, アイテムの重量の合計が $b$ を超えないようなアイテムの詰め合わせの中で， 価値の合計が最大のものを求めよ．ただし，各アイテムはナップサックに高々$1$個しか詰めることができないと仮定する．
$0$-$1$ナップサック問題は，アイテム $i (\in N)$ をナップサックに詰めるとき $1$， それ以外のとき $0$ になる $0$-$1$ 変数$x_i$を使うと，以下の$0$-$1$整数最適化問題として定式化できる． $$ \begin{array}{l l l} maximize & \sum_{i \in N} v_i x_i & \\ s.t. & \sum_{i \in N} s_i x_i \leq b & \\ & x_i \in \{0,1\} & i \in N \end{array} $$
再帰方程式を得るために， ナップサックに入れる対象を $k$番目までのアイテムに限定したとき， 重量の上限が $\theta$ 以下で価値の合計の最大のものを求める部分問題を考え，その最適値を $f(k,\theta)$ とする． もとの問題の最適値は $f(n,b)$ で計算できる．
部分問題間の関係より，$f(k,\theta)$ は，以下の再帰方程式によって計算可能である． $$ f(k,\theta) = \left\{ \begin{array}{ll} f(k-1,\theta) & 0 \leq \theta < s_k \\ \max \left\{ f(k-1,\theta), v_k + f(k-1, \theta - s_k) \right\} & s_k \leq \theta \end{array} \right\} k=1,2, \ldots,n $$
上の再帰方程式を条件 $$ f(0, \theta)= \left\{ \begin{array}{ll} 0 & 0 \leq \theta < s_0 \\ v_0 & s_0 \leq \theta \leq b \end{array} \right. $$ の下で解くことによって $f(k,\theta)$ を $O(k \theta)$時間で得ることができる． したがって，$0$-$1$ナップサック問題の最適値を $O(n b)$ 時間で求めることができる．
def kpdp(s, v, b):
def f(k, b):
if b < 0:
return -9999, 0
if k == 0:
if b >= s[0]:
memo[0, b] = v[0], 1
return memo[0, b]
else:
return 0, 0
if (k, b) in memo:
return memo[k, b]
else:
if f(k - 1, b)[0] < f(k - 1, b - s[k])[0] + v[k]:
max_value = f(k - 1, b - s[k])[0] + v[k]
sol = 1
else:
max_value = f(k - 1, b)[0]
sol = 0
memo[k, b] = max_value, sol
return memo[k, b]
memo = {}
opt_val, sol = f(len(s) - 1, b)
x = [0 for i in range(len(s))]
for k in range(len(s) - 1, -1, -1):
val, sol = memo[k, b]
if sol == 1:
x[k] += 1
b -= s[k] * sol
if b <= 0:
break
return opt_val, x
s = [2, 3, 4, 5]
v = [16, 19, 23, 28]
b = 7
opt_val, x = kpdp(s, v, b)
print("Opt. value=", opt_val)
print("Sol.=", x)
多制約ナップサック問題
多制約（$0$-$1$）ナップサック問題は，以下のように定義される．
$n$ 個のアイテムからなる有限集合 $N$，$m$ 本の制約の添え字集合 $M$， 各々のアイテム $j \in N$ の価値 $v_j (\geq 0)$， アイテム $j \in N$ の制約 $i \in M$ に対する重み $a_{ij}(\geq 0)$， および制約 $i \in M$ に対する制約の上限値 $b_i (\geq 0)$ が与えられたとき， 選択したアイテムの重みの合計が各制約 $i \in M$ の上限値 $b_i$ を超えないという条件の下で， 価値の合計を最大にするように $N$ からアイテムを選択する問題．
ナップサック問題は，アイテム $j (\in N)$ をナップサックに詰めるとき $1$， それ以外のとき $0$ になる $0$-$1$ 変数 $x_j$ を使うと，以下のように整数最適化問題として定式化できる． $$ \begin{array}{l l l} maximize & \sum_{j \in N} v_j x_j & \\ s.t. & \sum_{j \in N} a_{ij} x_j \leq b_i & \forall i \in M \\ & x_j \in \{0,1\} & \forall j \in N \end{array} $$
この問題は，制約が $1$本の問題（ナップサック問題）でも$NP$-困難である． ナップサック問題は，分枝限定法や動的最適化で容易に解くことができるが，制約の数が増えた場合には解くことが困難になる． 多くの市販の汎用の（混合）整数最適化ソルバーは，線形最適化問題に対する単体法や内点法と， 最も標準的に用いられる厳密解法である分枝限定法から構成されているが， 中規模の問題例でも求解が困難になることがある． 動的最適化は，与えられた数値データが小さな整数の場合には高速であるが，制約の数が増えたり，数値データがある程度大きな整数や実数であったりすると，求解が困難になる．
試しに，区間 $(0,1]$ の一様乱数 $U(0,1]$ を用いて以下のような（比較的難しいと言われている）問題例を作成してみた． 制約式の係数 $a_{ij}$ を $1- 1000 \log_2 U(0,1]$ とし，右辺定数 $b_i$ を $0.25 \sum_{j} a_{ij}$， 目的関数の係数 $v_j$ を $10 \sum_{i} a_{ij}/m+10 U(0,1]$ とする．
2007年に出版した本のために，そのとき使用していた市販の数理最適化ソルバーで試したことがある． 変数の数 $n=100$，制約の数 $m=5$ の問題例を，求解したところ，$2$ 時間かけても最適解を得ることができず，メモリ上限を超過してしまった． 本書を書くために，最新の混合整数最適化ソルバーGurobiで試したところ，この程度の問題例なら難なく最適解を得ることができた． もちろん，問題例の規模が増大すると厳密解を得ることが難しくなるが，制限時間内に良好な近似解を得るには十分な性能である．
多制約ナップサック問題のベンチマーク問題例は，以下のサイト(OR Library)から入手することができる．
ただし，このベンチマーク問題例は古く，最新の数理最適化ソルバーには簡単すぎる．
a = {}
b, v = [], []
m = 5
n = 100
for i in range(m):
sumA = 0.0
for j in range(n):
a[i, j] = 1.0 - math.log(random.random(), 2)
sumA += a[i, j]
b.append(0.25 * sumA)
for j in range(n):
sumA = 0.0
for i in range(m):
sumA += a[i, j]
v.append(10.0 * sumA / m + 10.0 * random.random())
import coptpy as cp
from coptpy import COPT, quicksum, tuplelist, multidict, tupledict
GRB = COPT
class Model():
def __init__(self, name=""):
env = cp.Envr()
self.ObjVal = None
self.model = env.createModel(name)
def addVar(self, *args, **kwargs):
return self.model.addVar(*args, **kwargs)
def addConstr(self, *args, **kwargs):
return self.model.addConstr(*args, **kwargs)
def setObjective(self, *args, **kwargs):
return self.model.setObjective(*args, **kwargs)
def getVars(self, *args, **kwargs):
return self.model.getVars(*args, **kwargs)
def getConstrs(self, *args, **kwargs):
return self.model.getConstrs(*args, **kwargs)
def getSOS(self, *args, **kwargs):
return self.model.SOS(*args, **kwargs)
def setParam(self, *args, **kwargs):
return self.model.setParam(*args, **kwargs)
def optimize(self):
ret = self.model.solve()
self.ObjVal = self.model.objval
self.status = self.model.status
return ret
J = list(range(n))
I = list(range(m))
model = Model("mkp")
x = {}
for j in J:
x[j] = model.addVar(vtype="B", name=f"x({j})")
#model.update()
constr = {}
for i in I:
for j in J:
constr[i,j] = model.addConstr(quicksum(a[i, j] * x[j] for j in J) <= b[i], f"Capacity({i},{j})")
model.setParam(GRB.Param.TimeLimit, 10.0)
model.setObjective(quicksum(v[j] * x[j] for j in J), GRB.MAXIMIZE)
model.optimize()
if model.status == GRB.OPTIMAL:
print("Objective value =", model.ObjVal)
print("Solution:")
for j in x:
print(f"x[{j}]={x[j].X}")