巡回セールスマン問題

巡回セールスマン問題（TSP）は、すべての頂点をちょうど1回ずつ訪問して出発点に戻る最短巡回路を求める問題です。 頂点は平面上に配置され、巡回路の長さはユークリッド距離で測られるものとします。

TSPのQUBO定式化

巡回路は頂点の順列として表現できます。 そこで、TSPの解を符号化するために置換行列を使用します。 \(X = (x_{ij})~(0 \leq i,j \leq n-1)\) を \(n \times n\) のバイナリ値の行列とします。行列 \(X\) は置換行列であり、各行と各列にちょうど1つの1が含まれます。 \(x_{ki}\) を「巡回路の \(k\) 番目の位置が頂点 \(i\) である」と解釈します。 したがって、\(X\) のすべての行とすべての列はone-hotでなければなりません。すなわち以下の制約が成り立つ必要があります:

\[\begin{split} \mathrm{row} : \sum_{j=0}^{n-1}x_{ij} &= 1 ~(0 \leq i \leq n-1),\\ \mathrm{column} : \sum_{i=0}^{n-1}x_{ij} &= 1 ~(0 \leq j \leq n-1). \end{split}\]

\(d_{ij}\) を頂点 \(i\) と \(j\) の間の距離とします。 置換行列 \(X\) に対する巡回路の長さは次のように書けます:

\[ \mathrm{objective} : \sum_{k=0}^{n-1} d_{ij}x_{ki}x_{(k+1)~{\rm mod}~n,~j}. \]

この式は、頂点 \(i\) が位置 \(k\) で訪問され、頂点 \(j\) が次の位置（\((k+1)~\mathrm{mod}~n\)）で訪問されるときにちょうど \(d_{ij}\) を加算するので、巡回路の総距離に等しくなります。

TSPのQUBO++プログラム 上記の置換行列による定式化を用いて、TSPのQUBO++プログラムを以下のように記述できます:

#define MAXDEG 2
#include <qbpp/qbpp.hpp>
#include <qbpp/easy_solver.hpp>
#include <qbpp/graph.hpp>

class Nodes {
  std::vector<std::pair<int, int>> nodes{{10, 12},  {33, 125},  {12, 226},
                                         {121, 11}, {108, 142}, {111, 243},
                                         {220, 4},  {210, 113}, {211, 233}};

 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 constraint = qbpp::sum(qbpp::vector_sum(x, 1) == 1) +
                    qbpp::sum(qbpp::vector_sum(x, 0) == 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);
  qbpp::Params params;
  params.set("time_limit", "1.0");
  auto sol = solver.search(params);
  auto tour = qbpp::onehot_to_int(sol(x));

  std::cout << "Tour: " << tour << "\n";
  auto graph = qbpp::graph::GraphDrawer();
  for (size_t i = 0; i < nodes.size(); ++i) {
    graph.add(qbpp::graph::Node(i).position(nodes[i].first, nodes[i].second));
  }
  for (size_t i = 0; i < nodes.size(); ++i) {
    auto from = tour[i];
    auto to = tour[(i + 1) % nodes.size()];
    graph.add(qbpp::graph::Edge(from, to).color("red").penwidth(2).directed());
  }
  graph.write("tsp_solution.svg");
}

このプログラムでは、頂点 0 から 8 の座標が Nodes オブジェクト nodes に格納されています。 バイナリ変数の2次元配列 x を作成し、one-hot制約と巡回路長の目的関数を構成します。 これらの項は、制約にペナルティ重み（ここでは 1000）を付けて加算することで、1つのQUBO式 f にまとめられます。実行可能性が優先されます。

次に、1.0秒の制限時間で EasySolver を使って f を解きます。 得られた割り当て sol(x) は置換行列を形成します。 この行列は qbpp::onehot_to_int() を使って整数のリスト（順列）tour に変換され、出力されます。 最後に、計算された tour が有向グラフとして描画され、ファイル tsp_solution.svg に保存されます。

このプログラムは以下の出力を生成します:

Tour: {7,8,5,2,4,1,0,3,6}

最初の頂点の固定

一般性を失うことなく、頂点0を巡回路の出発点と仮定できます。 TSPの巡回路は巡回シフトに対して不変であるため、出発位置を固定しても最適巡回路長は変わりません。

出発頂点を固定することで、QUBO式中のバイナリ変数の数を削減できます。 具体的には、置換行列において頂点0を位置0に割り当てることを強制します。 そのために、以下のバイナリ変数を固定します:

これらの割り当てにより、頂点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);
  qbpp::Params params;
  params.set("time_limit", "1.0");
  auto sol = solver.search(params);

  auto full_sol = qbpp::Sol(f).set(sol).set(ml);
  auto tour = qbpp::onehot_to_int(full_sol(x));
  std::cout << "Tour: " << tour << "\n";

まず、変数の固定割り当てを格納する qbpp::MapList オブジェクト ml を作成します。 各割り当ては push_back() メンバ関数で追加されます。

次に、qbpp::replace(f, ml) を呼び出します。これは ml で指定された固定値を元のQUBO式 f に代入して得られる新しい式を返します。 結果の式は g に格納され、簡約化されます。

次に、g に対するソルバーを作成して解 sol を得ます。 sol は縮小された問題に対応するため、f に対する qbpp::Sol オブジェクトを作成し、ソルバーの出力 sol と固定割り当て ml の両方を設定します。 結果の full_sol は x のすべての変数に対する完全な割り当てを格納します。

最後に、full_sol(x) で表される置換行列が qbpp::onehot_to_int() を使って順列に変換され、出力されます。

このプログラムは頂点0から始まる以下の巡回路を生成します:

Tour: {0,3,6,7,8,5,2,1,4}

TSPのPyQBPPプログラム

import math
import pyqbpp as qbpp

nodes = [(10, 12),  (33, 125),  (12, 226),
         (121, 11), (108, 142), (111, 243),
         (220, 4),  (210, 113), (211, 233)]

def dist(i, j):
    dx = nodes[i][0] - nodes[j][0]
    dy = nodes[i][1] - nodes[j][1]
    return round(math.sqrt(dx * dx + dy * dy))

n = len(nodes)
x = qbpp.var("x", n, n)

constraint = qbpp.sum(qbpp.vector_sum(x, 1) == 1) + qbpp.sum(qbpp.vector_sum(x, 0) == 1)

objective = qbpp.expr()
for i in range(n):
    next_i = (i + 1) % n
    for j in range(n):
        for k in range(n):
            if k != j:
                objective += dist(j, k) * x[i][j] * x[next_i][k]

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

solver = qbpp.EasySolver(f)
solver.time_limit(1.0)
sol = solver.search()

# Extract tour from permutation matrix
tour = []
for i in range(n):
    for j in range(n):
        if sol(x[i][j]) == 1:
            tour.append(j)
            break
print(f"Tour: {tour}")

このプログラムでは、ノードの座標をリストに格納しています。 バイナリ変数の2次元配列 x を作成し、one-hot制約と巡回路長の目的関数を構築します。

このプログラムの出力は以下のとおりです：

Tour: [7, 8, 5, 2, 4, 1, 0, 3, 6]

最初のノードの固定

一般性を失うことなく、ノード0を巡回路の開始ノードとすることができます。 開始ノードを固定することで、QUBO式のバイナリ変数の数を削減できます。

import pyqbpp as qbpp

ml = [(x[0][0], 1)]
ml += [(x[i][0], 0) for i in range(1, n)]
ml += [(x[0][i], 0) for i in range(1, n)]
ml = qbpp.MapList(ml)

g = qbpp.replace(f, ml)
g.simplify_as_binary()

solver = qbpp.EasySolver(g)
solver.time_limit(1.0)
sol = solver.search()

full_sol = qbpp.Sol(f).set(sol)
full_sol = full_sol.set(ml)

# Extract tour
tour = []
for i in range(n):
    for j in range(n):
        if full_sol(x[i][j]) == 1:
            tour.append(j)
            break
print(f"Tour: {tour}")

まず、変数の固定値を格納するペアのリスト ml を作成します。 次に replace(f, ml) を呼び出し、固定値を代入した新しい式を取得します。 sol は縮小された問題に対応するため、f に対する Sol オブジェクトを作成し、ソルバーの出力 sol と固定値 ml の両方を設定します。

このプログラムはノード0から始まる以下の巡回路を出力します：

Tour: [0, 3, 6, 7, 8, 5, 2, 1, 4]

C++ QUBO++との比較

C++ QUBO++	PyQBPP
qbpp::onehot_to_int(sol(x))	sol(x[i][j]) による手動ループ
qbpp::MapList ml	ml = [(x[0][0], 1)]
ml.push_back({x[0][0], 1})	ml.append((x[0][0], 1))
qbpp::replace(f, ml)	replace(f, ml)
qbpp::Sol(f).set(sol).set(ml)	Sol(f).set([sol, ml])

matplotlibによる可視化

以下のコードはTSPの解を可視化します：

import matplotlib.pyplot as plt

plt.figure(figsize=(6, 6))
for i, (nx_, ny) in enumerate(nodes):
    plt.plot(nx_, ny, "ko", markersize=8)
    plt.annotate(str(i), (nx_, ny), textcoords="offset points", xytext=(5, 5))

for i in range(n):
    fr = tour[i]
    to = tour[(i + 1) % n]
    plt.annotate("", xy=(nodes[to][0], nodes[to][1]),
                 xytext=(nodes[fr][0], nodes[fr][1]),
                 arrowprops=dict(arrowstyle="->", color="#e74c3c", lw=2))
plt.title("TSP Tour")
plt.savefig("tsp.png", dpi=150, bbox_inches="tight")
plt.show()

巡回路はノードを結ぶ赤い有向矢印で表示されます。