旅行販売員問題
巡回セールスマン問題（TSP）は、すべてのノードをちょうど1回ずつ訪問し、出発点に戻る最短の経路を求める問題です。 ノードは平面上に配置され、経路の長さはユークリッド距離で測定されると仮定します。
下図は、9つのノードと最適経路の例を示しています：
TSPのQUBO表現
経路はノードの順列で表現できる。 したがって、TSP解を符号化するために順列行列を用いる。
$X=(x_{i,j})$ ($0\leq i,j\leq n-1$) を $n\times n$ の二値行列とする。 行列 $X$ は置換行列であり、各行および各列には1つの要素のみが1となる（下図参照）。
$x_{k,i}$ を「巡回経路の $k$ 番目の位置がノード $i$ である」と解釈する。 したがって、$X$ の各行および各列はワンホットでなければならない。すなわち、以下の制約が成立する必要がある：
\[\begin{aligned} {\rm 行}:& \sum_{j=0}^{n-1}x_{i,j}=1 & (0\leq i\leq n-1)\\ {\rm 列}:& \sum_{i=0}^{n-1}x_{i,j}=1 & (0\leq j\leq n-1) \end{aligned}\]$d_{i,j}$ をノード $i$ と $j$ の間の距離と定義する。 すると、順列行列 $X$ の周長は次のように表せる：
\[\begin{aligned} {\rm 目的関数}: &\sum_{k=0}^{k-1} d_{i,j}x_{k,i}x_{(k+1)\bmod n,j} \end{aligned}\]この式は、位置 $k$ でノード $i$ を訪問し、次に位置 $(k+1)\bmod n$ でノード $j$ を訪問する際にのみ $d_{i,j}$ を加算するため、ツアーの総長に等しくなる。
TSPのためのQUBO++プログラム
上記の置換行列表現を用いて、TSPのQUBO++プログラムを以下のように記述できる：
#include "qbpp.hpp"
#include "qbpp_easy_solver.hpp"
class Nodes {
std::vector<std::pair<int, int>> nodes{{0, 0}, {0, 10}, {0, 20},
{10, 0}, {10, 10}, {10, 20},
{20, 0}, {20, 10}, {20, 20}};
public:
const std::pair<int, int>& operator[](std::size_t index) const {
return nodes[index];
}
std::size_t size() const { return nodes.size(); }
int dist(std::size_t i, std::size_t j) const {
auto [x1, y1] = nodes[i];
auto [x2, y2] = nodes[j];
const int dx = x1 - x2;
const int dy = y1 - y2;
return static_cast<int>(
std::llround(std::sqrt(static_cast<double>(dx * dx + dy * dy))));
}
};
int main() {
auto nodes = Nodes{};
auto x = qbpp::var("x", nodes.size(), nodes.size());
auto constraints = qbpp::sum(qbpp::vector_sum(x) == 1) +
qbpp::sum(qbpp::vector_sum(qbpp::transpose(x)) == 1);
auto objective = qbpp::expr();
for (size_t i = 0; i < nodes.size(); ++i) {
auto next_i = (i + 1) % nodes.size();
for (size_t j = 0; j < nodes.size(); ++j) {
for (size_t k = 0; k < nodes.size(); ++k) {
if (k != j) {
objective += nodes.dist(j, k) * x[i][j] * x[next_i][k];
}
}
}
}
auto f = objective + constraint * 1000;
f.simplify_as_binary();
auto solver = qbpp::easy_solver::EasySolver(f);
solver.time_limit(1.0);
auto sol = solver.search();
auto result = qbpp::onehot_to_int(sol(x));
std::cout << "Result: " << result << std::endl;
}
このプログラムでは、ノード0から8までの座標は Nodes オブジェクト nodesに格納される。
2次元配列 x を構築し、ワンホット制約と経路長目的関数を定義する。
これらを単一のQUBO式に統合する f として統合する。
実現可能性を優先させるため、ペナルティ重み（ここでは1000）を付加して制約を追加する。
その後、 f EasySolverを用いて1.0秒の時間制限で解きます。
結果の割り当て sol(x) は順列行列であり、 qbpp::onehot_to_int()によって整数リスト（順列）に変換され、出力される。
このプログラムは以下の出力を生成する:
Result: {7,8,5,2,1,0,3,4,6}
最初のノードを固定する
一般性を失うことなく、ノード0が巡路の開始ノードであると仮定できる。 TSP巡路は循環シフトに対して不変であるため、開始位置を固定しても最適巡路長は変わらない。
開始ノードを固定することで、QUBO式における二値変数の数を削減できる。 具体的には、ノード0が置換行列の0番目の位置に割り当てられることを強制する。 これを行うため、以下の二値変数を固定する：
\[\begin{aligned} x_{0,0} &= 1\\ x_{i,0} &= 0& (i\geq 1)\\ x_{0,j} &= 0& (j\geq 1) \end{aligned}\]これらの割り当てにより、ノード0は位置0にのみ出現し、他のノードが位置0に割り当てられることはありません。 結果として、実効的な問題規模が縮小され、一般的に局所探索ベースのソルバーによるQUBOの解法が容易になります。
最初のノードを固定するためのQUBO++プログラム
QUBO++では、固定変数割り当てを qbpp::replace() 関数を用いて固定変数の割り当てを適用できる：
qbpp::MapList ml;
ml.push_back({x[0][0], 1});
for (size_t i = 1; i < nodes.size(); ++i) {
ml.push_back({x[i][0], 0});
ml.push_back({x[0][i], 0});
}
auto g = qbpp::replace(f, ml);
g.simplify_as_binary();
auto solver = qbpp::easy_solver::EasySolver(g);
solver.time_limit(1.0);
auto sol = solver.search();
qbpp::Sol full_sol(f);
full_sol.set(sol);
full_sol.set(ml);
auto result = qbpp::onehot_to_int(full_sol(x));
std::cout << "Result: " << result << std::endl;
まず、 qbpp::MapList オブジェクト mlを作成します。
各代入は push_back() メンバ関数を使用して追加します。
次に、 qbpp::replace(f, ml)を呼び出します。これは、元のQUBO式に ml を元のQUBO式に代入して得られる新しい式を返します fに代入して得られる新しい式を返します。
結果の式は g に格納され、簡略化されます。
その後、 g に対するソルバーを作成し、解 sol。
ここで sol は簡約化された問題に対応するため、元の式に対する解を f を構築することで qbpp::Sol オブジェクト full_solを作成することで
元の式に対する解を再構築する。
ソルバー結果solと固定割り当てmlの両方が full_sol.
最後に、 full_sol(x) は qbpp::onehot_to_int() を用いて置換に変換され、出力される。
このプログラムはノード0から始まる以下のツアーを生成する:
Result: {0,3,6,7,8,5,2,1,4}