##巡回セールスマン問題にアニーリングマシンで挑む

一つ前の記事で、隣り合う都道府県が異なる色になるように、日本地図を塗り分ける問題を、アニーリングマシンで解く方法を紹介しました
計算機に高度な塗り絵をさせる話

アニーリングマシンで解く、ということは、つまりアニーリング方式量子コンピュータで解くことができる形にすることを目論んでいます。
もちろん、今私の手元の計算機で (Google Colaboratoryを使用) 解ける計算量のものなら良いのですが、実世界の問題を解くことを考えると、計算量が半端ないので、もっと速い量子コンピュータを用いなくてはいけません。
量子コンピュータはまだ実用段階ではありませんが、アニーリング(焼きなましの方法)方式で問題を表現する方法は考えられています。
私は、量子コンピュータを専門にされる先生のゼミに参加して、その方法を学んだので、その一部分をQiitaにまとめています。

それでは行きましょう！巡回セールスマン問題！！

##巡回セールスマン問題

巡回セールスマン問題とは、各都市をセールスマンが巡回し、出発地点に戻る最短経路を求める問題のことです。

この黄色い点を効率よく通るにはどうしますか？？

わからない！！！
計算機に任せましょう。
まずは、この問題を計算機に投げるために、QUBO形式で表現することを考えます。
###QUBO形式とは

QUBOとは、Quadratic Unconstrained Binary Optimization (制約なし二次形式二値変数最適化)の略である。
0または1, あるいは +1または-1といった二値のみを取れる被最適化変数xiに対して、目的関数が変数の2次の項までで記述ができ（(xi)(xj)(xk)などの3次以上の項を含まない）、
かつ、∑Ni=1xi<Cのような明示的な制約を含まない問題に対する最適化のことを意味している。

(引用元：東北大学量子アニーリング研究開発センターT-QARDが運営するD-Waveマシンのユースケースナレッジベースサイト)

前回は省略してしまったので、今回は巡回セールスマン問題をQUBO形式で表した数式を置いておきます。

##計算機に託そう
###下準備
まずはパッケージの読み込みから。

# イジング模型およびQUBO形式の問題を解くためのパッケージ
!pip install dimod

# イジング模型およびQUBO形式の問題を解くためのパッケージ
import dimod
# NumPy
import numpy as np
# グラフ（ネットワーク）用のパッケージ
import networkx as nx
# 図の描画の関係パッケージ
%matplotlib inline
import matplotlib.pyplot as put

この問題では、都市の配置をランダムに設定してプロットし、それらの効果的な回り方を算出する。
つまり、回り方を、有向グラフを用いて表す。
そのために、そのグラフをプロットするための関数を用意しておく。

def plot_city(cities, sample = {}):
    N = len(cities)
    cities_dict = dict(cities)
    G = nx.DiGraph() # 有向グラフ
    # グラフの頂点（都市）を設定
    for city in cities_dict:
        G.add_node(city)
    
    # 経路
    if sample:
        city_order = [-1]*N
        # 都市(i)と順番(a)を解の情報からとりだす
        for key, val in sample.items():
          i = key // N
          a = key % N
          if val == 1:
            city_order[i] = a
        # 経路を設定する
        for i in range(N):
          a = city_order[i]
          if a < 0: # 都市iを訪問しない場合はとばす
            continue
          for j in range(N):
            if i == j:
              continue
            if city_order[j] == (a + 1) % N: # 都市jが次の訪問地なら辺を設定
              G.add_edge(cities[i][0], cities[j][0])
              
    # グラフの描画
    plt.figure(figsize=(3,3))
    #pos = nx.spring_layout(G)
    nx.draw_networkx(G, cities_dict, node_color='orange')
    plt.axis("off")
    plt.show()

そして、以下のように、頂点位置を設定した5つの都市をプロットしてみる。
(もちろん、Nの値を変えて、cities配列の中の頂点を設定したり、削除して都市の数を調整することは可能です)

N = 5
cities = [
    ("A", (0, 0)),
    ("B", (0.3, 0.9)),
    ("C", (0.8, 0.6)),
    ("D", (0.5, 0.3)),
    ("E", (0.1, 0.4))
]
plot_city(cities)

プロットすると下のようになる。

まずは、上で設定した都市間の距離を、二点間距離を求める方法で調べる。

dij = np.zeros((N, N))
for i in range(N):
  for j in range(N):
    pos_i = cities[i][1]
    pos_j = cities[j][1]
    dij[i, j] = np.sqrt((pos_i[0] - pos_j[0])**2 + (pos_i[1] - pos_j[1])**2)

いよいよ準備が整ったので、最短経路を求めていきます！
###数式のお話
アニーリングマシンのソルバーでは、

という形で2次の項の情報 Jkl と1次の項の情報 hk を与えることになる。
よって、この記事の冒頭にQUBO形式で書いたハミルトニアン(下の式)を、上の形に変形しなくてはならない。

これは数式変換の部分なので、説明の資料を掲載するに留めます。
興味がある人は追ってみてください。

###数式をコードにする

def Hamiltonian(N, Jij, alpha = 1.0):
  """ハミルトニアンを計算してQUBO形式で返す
  N: 都市の数
  Jij: 点i,j間の相互作用の強さ（今回は2都市間距離）
  alpha: 制約の強さ 
  """
  A = np.zeros((N*N, N*N)) # N都市
  # 第1項：問題を表す
  for i in range(N): # 0 <= i < N
    for j in range(N): # 0 <= j < N
      if i == j:
        continue # i = j を回避
      for a in range(N-1): # 巡る順番
        k1 = N*i + a
        k2 = N*j + a + 1
        A[k1, k2] += Jij[i, j]
      # a = N-1 のときの処理：出発地に戻る
      k1 = N*i + N - 1
      k2 = N*j
      A[k1, k2] += Jij[i, j]
  # 第2項：各都市を1回ずつ訪問する制約条件
  for i in range(N):
    for a in range(N):
      k1 = N*i + a
      for b in range(N):
        k2 = N*i + b
        if k1 == k2: # 1次の項は対角成分に入れる
          A[k1, k2] -= alpha
        else:
          A[k1, k2] += alpha
  # 第3項：同時に複数の都市を訪問できない制約条件
  for a in range(N):
    for i in range(N):
      k1 = N*i + a
      for j in range(N):
        k2 = N*j + a
        if k1 == k2: # 1次の項は対角成分に入れる
          A[k1, k2] -= alpha
        else:
          A[k1, k2] += alpha
  # dimodソルバー用のQUBO形式に変換
  Q = {}
  for k1 in range(N*N):
    for k2 in range(N*N):
      if A[k1, k2] != 0:
        Q[(k1, k2)] = A[k1, k2]
  return Q

###実行してみよう
準備が少し長かったが、各関数の準備ができたので、上で設定した(頂)点を最短で回るためのルートを求めてみましょう

# ハミルトニアンの設定
Q = Hamiltonian(N, dij, alpha = 1)
# ソルバーの設定：シミュレーテッド・アニーリング
solver = dimod.SimulatedAnnealingSampler()
# ソルバーによる求解
response = solver.sample_qubo(Q, num_reads = 10)
# 解とエネルギーを取り出す
sample0, energy0 = next(response.data(fields=['sample', 'energy']))
# 結果のプロット
plot_city(cities, sample0)

おお…！予想通りすぎてあまり面白くないので、点を10個にしましょう。

wow
計算機様様ですね！