QUBO＋＋入門：巡回セールスマン問題

巡回セールスマン問題のQUBO化

巡回セールスマン問題（Traveling Salesman Problem, TSP）とは，複数のノードと全ノード間の距離が与えられたとき，すべてのノードを一度ずつ訪問し，最初のノードに戻るツアー（巡回路）のうち，総移動距離が最小となるものを求める問題です．本稿では，簡単のため，2次元平面上の整数座標を持つ
$n$個のノードが与えられている場合を考えます．ノード間の距離としては，ユークリッド距離を整数に丸めた値を使用します．

この巡回セールスマン問題を，QUBO++を用いてQUBO式で表現し，その解を求めるプログラムを設計します．例えば，以下は $n=8$ の場合における入力ノードの例です．この8つのノードを巡り，元のノードに戻るツアーを巡回セールスマン問題の解とし，そのうち総移動距離が最小となるものを探索します．

巡回セールスマン問題の解は，全ノードの順列で表すことができます．そこで，

• 前の記事：QUBO＋＋入門：行列操作（順列生成）
  で使用した順列生成のためのQUBO式を利用します．二値変数で構成される大きさ $n\times n$の行列
  $x$を用いて，$n$個のノードの順列を表現します．この$x$の各行および各列がone-hotベクトルとなり，順列を表します．順列に基づく総移動距離は，次の式で計算できます

$$\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}\sum_{k=0}^{n-1} d_{j,k}x_{i,j}x_{(i+1)\bmod n,k}$$

ここで，$d_{i,j}$はノード$i$とノード$j$の距離を表します．$x{i,j}=1$かつ$x_{(i+1)\bmod n,k}=1$の場合，すなわち，順列で$i$番目のノードが$j$で，次の$(i+1)\bmod n$がノード$k$であるとき，移動距離$d_{j,k}$が加算されます．この総移動距離が最小となる$x$の値を求めれば，それが巡回セールスマンの最適解を表します．

QUBO＋＋プログラム

先の9ノードの場合に巡回セールスマン問題を解くQUBO＋＋プログラムは以下のようになります．
二値変数が64個となり，Exhaustive Solverで解くには多すぎるため，Easy Solverを使用します．

#include "qbpp.hpp"
#include "qbpp_easy_solver.hpp"

int distance(const std::pair<int, int> &p1, const std::pair<int, int> &p2) {
  return static_cast<int>(
      std::round(std::sqrt((p1.first - p2.first) * (p1.first - p2.first) +
                           (p1.second - p2.second) * (p1.second - p2.second))));
}

int main() {
  std::vector<std::pair<int, int>> nodes = {
      {0, 0}, {10, 0}, {20, 0}, {0, 10}, {10, 10}, {20, 10}, {0, 20}, {10, 20}};

  auto x = qbpp::var("x", nodes.size(), nodes.size());
  auto permutation = qbpp::sum(qbpp::sum(x) == 1) +
                     qbpp::sum(qbpp::sum(qbpp::transpose(x)) == 1);

  auto dist = qbpp::expr();
  for (size_t i = 0; i < nodes.size(); i++) {
    for (size_t j = 0; j < nodes.size(); j++) {
      for (size_t k = 0; k < nodes.size(); k++) {
        if (j == k) continue;
        dist += distance(nodes[j], nodes[k]) * x[i][j] *
                x[(i + 1) % nodes.size()][k];
      }
    }
  }

  auto f = simplify_as_binary(permutation * 100 + dist);

  auto quad_model = qbpp::QuadModel(f);
  std::cout << "var_count = " << quad_model.var_count() << std::endl;
  std::cout << "linear_term_count = " << quad_model.term_count(1) << std::endl;
  std::cout << "quadratic_term_count = " << quad_model.term_count(2)
            << std::endl;

  auto solver = qbpp::easy_solver::EasySolver(quad_model);
  solver.set_time_limit(1);
  auto sol = solver.search();

  std::cout << sol << std::endl;
  auto tour = qbpp::onehot_to_int(sol.get(x));
  std::cout << str(tour, "tour") << std::endl;
}

• 関数distanceは，2つのノードの座標を引数として受け取り，その距離を計算して返します．
• 変数nodesには8個のノードの座標が格納されます．
• xを二値変数の$9\times 9$行列として生成します．
• permutaitonはxの各行と各列がone-hotベクトル，つまり順列を表すときに最小値0をとるQUBO式を代入しています．
• 三重のforループを用いて，順列により表されるツアーの総移動距離を計算するQUBO式distを求めています．
• QUBO式fは，permutationを優先するために100を掛けた後，distを加算しています．
• fをqbpp::QuadModelオブジェクトquad_modelに変換し，変数，一次項，二次項の個数を表示しています．
• Easy Solverを用いて，解を探索しています．ここで制限時間を1秒に設定しています．
• 解xが保持する順列をtourに代入し，表示しています．

このプログラムを実行すると次の結果が得られます．

var_count = 64
linear_term_count = 64
quadratic_term_count = 1792
80:{{x[0][0],0},{x[0][1],0},{x[0][2],0},{x[0][3],0},{x[0][4],1},{x[0][5],0},{x[0][6],0},{x[0][7],0},{x[1][0],0},{x[1][1],0},{x[1][2],0},{x[1][3],0},{x[1][4],0},{x[1][5],0},{x[1][6],0},{x[1][7],1},{x[2][0],0},{x[2][1],0},{x[2][2],0},{x[2][3],0},{x[2][4],0},{x[2][5],0},{x[2][6],1},{x[2][7],0},{x[3][0],0},{x[3][1],0},{x[3][2],0},{x[3][3],1},{x[3][4],0},{x[3][5],0},{x[3][6],0},{x[3][7],0},{x[4][0],1},{x[4][1],0},{x[4][2],0},{x[4][3],0},{x[4][4],0},{x[4][5],0},{x[4][6],0},{x[4][7],0},{x[5][0],0},{x[5][1],1},{x[5][2],0},{x[5][3],0},{x[5][4],0},{x[5][5],0},{x[5][6],0},{x[5][7],0},{x[6][0],0},{x[6][1],0},{x[6][2],1},{x[6][3],0},{x[6][4],0},{x[6][5],0},{x[6][6],0},{x[6][7],0},{x[7][0],0},{x[7][1],0},{x[7][2],0},{x[7][3],0},{x[7][4],0},{x[7][5],1},{x[7][6],0},{x[7][7],0}}
{tour[0],4},{tour[1],7},{tour[2],6},{tour[3],3},{tour[4],0},{tour[5],1},{tour[6],2},{tour[7],5}

変数の個数は64個，一次項の個数は64個，二次項の個数は1792個です．Easy Solverを用いることで，エネルギー値80の解が得られています．このエネルギー値はツアーの総移動距離を表します．得られたツアーは，次の順序でノードを巡るものです．ツアーは，4,7,6,3,0,1,2,5の順でノードを巡回します．このツアーを図示すると，最適解が得られていることが確認できます．

まとめ

• $n$ノードの巡回セールスマン問題をQUBO形式で解く場合，大きさ$n \times n$の二値変数行列を用います．
• このQUBO式は，各行および各列がone-hotベクトルとなる制約と，総移動距離を最小化する目的関数を組み合わせたものです．
• 変数の個数が多い場合，Exhaustive Solverでは現実的な時間内に探索が終了しないため，Easy Solverを使用します．

QUBO++情報

• QUBO＋＋ツールとABS2 QUBO Solverの詳細情報：QUBO++に関する公式サイトです．ツールの特徴や使用方法，ABS2 QUBO Solverとの連携について詳しく説明されています．
• 【研究成果】組合せ最適化問題を解くためのプログラミング用ツールQUBO＋＋を無償公開します：広島大学からの公式リリース情報です．QUBO++の開発背景や研究成果について記載されています．
• 利用条件：QUBO++は非商用利用や評価研究目的に限り無償で利用可能です．商用利用や定常的な業務での使用を希望する場合は，ライセンス契約が必要となります．