ナップサック問題

重さと価値を持つアイテムの集合と,重量制限のあるナップサックが与えられたとき,ナップサック問題は,総重量が容量以内に収まるようにしつつ,総価値を最大化するアイテムの部分集合を選択することを目的とします.

\(w_i\) と \(v_i\) (\(0\leq i\leq n-1\)) をそれぞれアイテム \(i\) の重さと価値とします. \(S\in \lbrace 0, 1, \ldots n-1\rbrace\) を選択されたアイテムの集合とします.

\[\begin{split} \begin{aligned} \text{Maximize:} & \sum_{i\in S} v_i \\ \text{Subject to:} & \sum_{i\in S} w_i \leq W \end{aligned} \end{split}\]

ここで \(W\) はナップサックの重量容量です.

QUBO定式化

この問題をQUBOとして定式化するために,\(n\) 個のバイナリ変数 \(x_i\in\lbrace 0,1\rbrace\) (\(0\leq i\leq n-1\)) の集合 \(X\) を導入します.ここで,アイテム \(i\) が選択されるのは \(x_i=1\) のときかつそのときに限ります.

上記の定式化は次のように書き換えられます：

\[\begin{split} \begin{aligned} \text{Maximize:} & \sum_{i=0}^{n-1} v_ix_i \\ \text{Subject to:} & \sum_{i=0}^{n-1} w_ix_i \leq W \end{aligned} \end{split}\]

Hi-QUBOプログラム

制約はHi-QUBOが提供する範囲演算子を用いて表現できます. 結果として得られるQUBO目的関数は次のように定義されます：

\[ \begin{aligned} f(X) &= -\sum_{i=0}^{n-1} v_ix_i + P\times (0\leq \sum_{i=0}^{n-1} w_ix_i \leq W) \end{aligned} \]

QUBOソルバーは目的関数を最小化するため,元の最大化目的は符号を反転しています. 定数 \(P\) は制約を強制するための十分大きなペナルティパラメータです.

以下のHi-QUBOプログラムは,Exhaustive Solverを用いて10個のアイテムのナップサック問題を解きます：

knapsack-program1.cpp

#include <qbpp/qbpp.hpp>
#include <qbpp/exhaustive_solver.hpp>

int main() {
  auto w = qbpp::array({10, 20, 30, 5, 8, 15, 12, 7, 17, 18});
  auto v = qbpp::array({60, 100, 120, 60, 80, 150, 110, 70, 150, 160});
  int capacity = 50;

  auto x = qbpp::var("x", w.size());

  auto constraint = 0 <= qbpp::sum(w * x) <= capacity;
  auto objective = qbpp::sum(v * x);

  auto f = -objective + 1000 * constraint;
  f.simplify_as_binary();

  auto solver = qbpp::ExhaustiveSolver(f);
  auto sols = solver.search({{"best_energy_sols", 0}});
  for (size_t i = 0; i < sols.size(); ++i) {
    const auto& sol = sols.sols[i];
    std::cout << "[Solution " << i << "]" << std::endl;
    std::cout << "Energy = " << sol.energy() << std::endl;
    std::cout << "Constraint  = " << constraint.body(sol) << std::endl;
    std::cout << "Objective  = " << sol(objective) << std::endl;
    for (size_t j = 0; j < w.size(); ++j) {
      if (sol(x[j]) == 1) {
        std::cout << "Item " << j << ": weight = " << w[j]
                  << ", value =  " << v[j] << std::endl;
      }
    }
  }
}

このプログラムでは,式 constraint と objective を別々に構築し,ペナルティ係数 1000 を用いて最終的なQUBO式 f に結合しています. 次に,Exhaustive Solver を f に適用し,すべての最適解を列挙します.

以下の出力は,エネルギー,制約値,目的関数値を含む最適解を示しています：

[Solution 0]
Energy = -480
Constraint  = 50
Objective  = 480
Item 3: weight = 5, value =  60
Item 5: weight = 15, value =  150
Item 6: weight = 12, value =  110
Item 9: weight = 18, value =  160
[Solution 1]
Energy = -480
Constraint  = 50
Objective  = 480
Item 3: weight = 5, value =  60
Item 4: weight = 8, value =  80
Item 6: weight = 12, value =  110
Item 7: weight = 7, value =  70
Item 9: weight = 18, value =  160

このインスタンスには2つの最適解があり,いずれも総価値 480 を達成しつつ,容量制約をちょうど満たしていることがわかります.

qbpp::cons() で容量制約を表す

容量制約は,同じ範囲式を qbpp::cons() で囲むことで制約として印を付けられます. 上のプログラムからの変更はこれだけです.constraint を qbpp::cons(0 <= qbpp::sum(w * x) <= capacity) にするだけで,constraint.body(sol) を含む残りの部分はそのままです.バンドルされたソルバーはこれを制約として扱い, Exhaustive Solver では容量を満たす選択だけが列挙されます:

knapsack-program2.cpp

#include <qbpp/qbpp.hpp>
#include <qbpp/exhaustive_solver.hpp>

int main() {
  auto w = qbpp::array({10, 20, 30, 5, 8, 15, 12, 7, 17, 18});
  auto v = qbpp::array({60, 100, 120, 60, 80, 150, 110, 70, 150, 160});
  int capacity = 50;

  auto x = qbpp::var("x", w.size());

  auto constraint = qbpp::cons(0 <= qbpp::sum(w * x) <= capacity);
  auto objective = qbpp::sum(v * x);

  auto f = -objective + 1000 * constraint;
  f.simplify_as_binary();

  auto solver = qbpp::ExhaustiveSolver(f);
  auto sols = solver.search({{"best_energy_sols", 0}});
  for (size_t i = 0; i < sols.size(); ++i) {
    const auto& sol = sols.sols[i];
    std::cout << "[Solution " << i << "]" << std::endl;
    std::cout << "Energy = " << sol.energy() << std::endl;
    std::cout << "Constraint  = " << constraint.body(sol) << std::endl;
    std::cout << "Objective  = " << sol(objective) << std::endl;
    for (size_t j = 0; j < w.size(); ++j) {
      if (sol(x[j]) == 1) {
        std::cout << "Item " << j << ": weight = " << w[j]
                  << ", value =  " << v[j] << std::endl;
      }
    }
  }
}

このプログラムは,範囲演算子版と同じ2つの最適解を出力します:

[Solution 0]
Energy = -480
Constraint  = 50
Objective  = 480
Item 3: weight = 5, value =  60
Item 5: weight = 15, value =  150
Item 6: weight = 12, value =  110
Item 9: weight = 18, value =  160
[Solution 1]
Energy = -480
Constraint  = 50
Objective  = 480
Item 3: weight = 5, value =  60
Item 4: weight = 8, value =  80
Item 6: weight = 12, value =  110
Item 7: weight = 7, value =  70
Item 9: weight = 18, value =  160

2つの定式化は等価です.制約を qbpp::cons() で囲むことで,バンドルされた ソルバーがそれを制約として扱うため,より大きなナップサック問題も扱いやすく なります.

重さと価値を持つアイテムの集合と,重量制限のあるナップサックが与えられたとき,ナップサック問題は,総重量が容量以内に収まるようにしつつ,総価値を最大化するアイテムの部分集合を選択することを目的とします.

\(w_i\) と \(v_i\)（\(0\leq i\leq n-1\)）をそれぞれアイテム \(i\) の重さと価値とします. \(S\in \lbrace 0, 1, \ldots n-1\rbrace\) を選択されたアイテムの集合とします.

\[\begin{split} \begin{aligned} \text{Maximize:} & \sum_{i\in S} v_i \\ \text{Subject to:} & \sum_{i\in S} w_i \leq W \end{aligned} \end{split}\]

ここで \(W\) はナップサックの重量容量です.

QUBO定式化

この問題をQUBOとして定式化するために,\(n\) 個のバイナリ変数 \(x_i\in\lbrace 0,1\rbrace\)（\(0\leq i\leq n-1\)）の集合 \(X\) を導入します.ここで,アイテム \(i\) が選択されるのは \(x_i=1\) のときかつそのときに限ります.

上記の定式化は次のように書き換えられます：

\[\begin{split} \begin{aligned} \text{Maximize:} & \sum_{i=0}^{n-1} v_ix_i \\ \text{Subject to:} & \sum_{i=0}^{n-1} w_ix_i \leq W \end{aligned} \end{split}\]

PyQBPPプログラム

制約は PyQBPP が提供する範囲演算子 (lo <= expr) & (qbpp.same <= hi) を用いて表現できます. 結果として得られるQUBO目的関数は次のように定義されます：

\[ \begin{aligned} f(X) &= -\sum_{i=0}^{n-1} v_ix_i + P\times (0\leq \sum_{i=0}^{n-1} w_ix_i \leq W) \end{aligned} \]

QUBOソルバーは目的関数を最小化するため,元の最大化目的は符号を反転しています. 定数 \(P\) は制約を強制するための十分大きなペナルティパラメータです.

以下のPyQBPPプログラムは,Exhaustive Solverを用いて10個のアイテムのナップサック問題を解きます：

knapsack-program1.py

import pyqbpp as qbpp

w = qbpp.array([10, 20, 30, 5, 8, 15, 12, 7, 17, 18])
v = qbpp.array([60, 100, 120, 60, 80, 150, 110, 70, 150, 160])
capacity = 50

x = qbpp.var("x", shape=len(w))

constraint = (0 <= qbpp.sum(w * x)) & (qbpp.same <= capacity)
objective = qbpp.sum(v * x)

f = -objective + 1000 * constraint
f.simplify_as_binary()

solver = qbpp.ExhaustiveSolver(f)
result = solver.search(best_energy_sols=0)
for idx, sol in enumerate(result.sols):
    print(f"[Solution {idx}]")
    print(f"Energy = {sol.energy}")
    print(f"Constraint = {sol(constraint.body)}")
    print(f"Objective = {sol(objective)}")
    for j in range(len(w)):
        if sol(x[j]) == 1:
            print(f"Item {j}: weight = {w[j]}, value = {v[j]}")

このプログラムでは,式 constraint と objective を別々に構築し,ペナルティ係数 1000 を用いて最終的なQUBO式 f に結合しています. 次に,Exhaustive Solver を f に適用し,すべての最適解を列挙します.

以下の出力は,エネルギー,制約値,目的関数値を含む最適解を示しています：

[Solution 0]
Energy = -480
Constraint = 50
Objective = 480
Item 3: weight = 5, value = 60
Item 5: weight = 15, value = 150
Item 6: weight = 12, value = 110
Item 9: weight = 18, value = 160
[Solution 1]
Energy = -480
Constraint = 50
Objective = 480
Item 3: weight = 5, value = 60
Item 4: weight = 8, value = 80
Item 6: weight = 12, value = 110
Item 7: weight = 7, value = 70
Item 9: weight = 18, value = 160

このインスタンスには2つの最適解があり,いずれも総価値 480 を達成しつつ,容量制約をちょうど満たしていることがわかります.

qbpp.cons() で容量制約を表す

容量制約は,同じ範囲式を qbpp.cons() で囲むことで制約として印を付けられます. 上のプログラムからの変更はこれだけです.constraint を qbpp.cons((0 <= qbpp.sum(w * x)) & (qbpp.same <= capacity)) にするだけで, sol(constraint.body) を含む残りの部分はそのままです.バンドルされたソルバーは これを制約として扱い,Exhaustive Solver では容量を満たす選択だけが列挙されます:

knapsack-program2.py

import pyqbpp as qbpp

w = qbpp.array([10, 20, 30, 5, 8, 15, 12, 7, 17, 18])
v = qbpp.array([60, 100, 120, 60, 80, 150, 110, 70, 150, 160])
capacity = 50

x = qbpp.var("x", shape=len(w))

constraint = qbpp.cons((0 <= qbpp.sum(w * x)) & (qbpp.same <= capacity))
objective = qbpp.sum(v * x)

f = -objective + 1000 * constraint
f.simplify_as_binary()

solver = qbpp.ExhaustiveSolver(f)
result = solver.search(best_energy_sols=0)
for idx, sol in enumerate(result.sols):
    print(f"[Solution {idx}]")
    print(f"Energy = {sol.energy}")
    print(f"Constraint = {sol(constraint.body)}")
    print(f"Objective = {sol(objective)}")
    for j in range(len(w)):
        if sol(x[j]) == 1:
            print(f"Item {j}: weight = {w[j]}, value = {v[j]}")

このプログラムは,範囲演算子版と同じ2つの最適解を出力します:

[Solution 0]
Energy = -480
Constraint = 50
Objective = 480
Item 3: weight = 5, value = 60
Item 5: weight = 15, value = 150
Item 6: weight = 12, value = 110
Item 9: weight = 18, value = 160
[Solution 1]
Energy = -480
Constraint = 50
Objective = 480
Item 3: weight = 5, value = 60
Item 4: weight = 8, value = 80
Item 6: weight = 12, value = 110
Item 7: weight = 7, value = 70
Item 9: weight = 18, value = 160

2つの定式化は等価です.制約を qbpp.cons() で囲むことで,バンドルされた ソルバーがそれを制約として扱うため,より大きなナップサック問題も扱いやすく なります.