Watsonで巡回セールスマン問題を解く

# はじめに
[Watsonで数独を解く! Decision Optimizerを使ってみた](https://qiita.com/makaishi2/items/d1cc9a3f49f640a3b649)に引き続きDecision Optimizerシリーズ第二段です。
今回はあまりにも有名な問題「**巡回セールスマン問題**」にチャレンジしてみます。

[2020-03-16 githubのリポジトリ移動]

# 巡回セールスマン問題とは
　一人のセールスマンが、N箇所の場所を一筆書きで回りたい。この場合に最短時間で回れるコースをどうやってみつけるか？
という問題です。
数学的には「NP困難」と呼ばれる領域の問題で、少しNが大きくなると、調べるべき組み合わせの数が爆発的に増大し、完全解を求められないことがわかっています。より詳しい話は、下記のWikipediaの記事を参照してください。

[巡回セールスマン問題](https://ja.wikipedia.org/wiki/巡回セールスマン問題)

# CPLEXの2つのライブラリ
これから、この問題を**CPLEX**を使って解いていくのですが、その前にCPLEXの2種類のライブラリについて説明しておきます。
**Mathematical Programming Modeling** (MP) と呼ばれるライブラリと、**Constraint Programming Modeling** (CP) と呼ばれるライブラリがその2つです。
前回紹介した「**数独**」は、CPで解く種類の問題です。
前回の記事の最後に紹介した「**原油のブレンド問題**」は、MPで解く問題となります。違いは、MPの場合は**目的関数**と呼ばれる関数があり、この関数の値を最大にする組み合わせをみつける問題である点となります。今回解く、巡回セールスマン問題も、**移動距離合計**という目的関数が存在するのでMPで解くべき問題ということになります。

(注)一点補足すると、上で説明した「NP困難」の数学的事実はCPLEXでも解決することはできません。
CPLEXは「NP困難な問題に対しても実用上問題のない次善の解(最適解ではない)を、実用になる有限時間で見つけるツール」と理解してもらえればいいです。
もう一点補足すると、最近注目を浴びている「量子コンピュータ」は、まったく新しい技術要素を用いて、この問題に対する最適解を有限時間で求めようとしているアプローチということになります。

# 数学的原理
CPLEXで問題を解く場合、制約や目的関数など、数学的な原理を事前に明確にしておく必要があります。
前回紹介した数独では、この原理はほぼ自明なものだったのですが、今回は結構数式としての記述が難しい問題になります。
この点に関しては、下記リンク先の情報を参考にさせてもらいました。

http://satemochi.blog.fc2.com/blog-entry-24.html

ただ、条件式に関しては、コードで表現しやすいように、一部同値な式で書き直しています。
($x_{ii} = 0$という制約を付け加えることで、条件式を簡略化した)

最終的な数式を書き直すと、次のようになります。
この式の中で$x_{ij}とu_{i}$がどういう意味を持つかについては、下のPythonコードのコメントのところに書いておきましたので、そちらを参照して下さい。
なんとなく$u_iはいらなくてx_{ij}$だけでいいのではという気がするのですが、この変数は、2個以上の小さなループの集合になる解を排除して、全体で一つの大きなループにするために必要な変数です。このことを記述しているのが、制約式の一番最後にある謎の式で、どうしてこれでいいかについては、上のリンク先記事を参照して下さい。

$ x_{ii} = 0 $  　$(i = 0, 1, .. N-1)$
$ 0 \leqq  x_{ij} \leqq 1$　　$i\ne j$　　$(i, j = 0, 1, .., N-1)$

$\displaystyle\sum_{i=0}^{N-1} x_{ij} = 1$　　$(j = 0, 1, .. N-1)$
$\displaystyle\sum_{j=0}^{N-1} x_{ij} = 1$　　$(i = 0, 1, .. N-1)$

$ u_0 = 0$
$ 1 \leqq u_i \leqq N-1$　　$(i = 1, 2, .., N-1)$

$u_i - u_j + N \cdot x_{ij} \leqq N-1$　　$(i, j = 0, 1, .., N-1)$

また、目的関数は次の式になります。

$ cost = \displaystyle\sum_{i=0}^{N-1}\displaystyle\sum_{j=0}^{N-1}
c_{ij}x_{ij}
$

# Jupyter Notebookコード
ここまでの準備ができれば、後は上の数式をPython APIで表現するだけです。
実装コードの全体は、
https://github.com/makaishi2/sample-notebooks/blob/master/cplex/TSP.ipynb
にアップしてあります。

また、前の記事のようにURL指定でWatson StudioのJupyter Notebookにロードする場合は、次のURLを指定して下さい。

https://github.com/makaishi2/sample-notebooks/raw/master/cplex/TSP.ipynb

# 結果概要

オリジナルデータは48件分あります。
Notebookのコードでは、変数Nの設定で、このうち何件分のデータで最適化を行うか、指定できるようになっています。
かかった時間は25件の場合で約10秒、30件の場合で約4分弱でした。

30件の場合の元データ散布図と、それに対してCPLEXが示した解によるルート図を以下に示します。

![スクリーンショット 2019-08-16 19.02.02.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/206130/7fb52afc-9de1-13f1-f22a-8734929058a0.png)


![スクリーンショット 2019-08-16 19.02.12.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/206130/1c7d52cc-05d8-d900-9590-4c40ae8640e0.png)


# コード解説
いつものようにコードの簡単な解説を以下に示します。

##  問題データのロード

### ライブラリロード
```py3
# ライブラリのロード

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from IPython.display import display
```

### データのロード

オリジナル問題の48件分のデータを、githubからロードする作りにしました。

```py3
# データのロード
# 48件は巡回セールスマン問題のオリジナルデータです

url = 'https://raw.githubusercontent.com/makaishi2/sample-data/master/data/att48.csv'
df = pd.read_csv(url)
display(df.head())
print('size = ', df.values.shape)
```

下の図のような結果が返ってくるはずです。

![スクリーンショット 2019-08-16 19.09.41.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/206130/8f1e6beb-d686-f374-aeb4-9e17d71fad66.png)

### 対象項目絞り込み

Nの変数設定で、分析対象項目の絞り込みを行います。絞り込んだ結果に関しては散布図表示しています。

```py3
# 対象項目数
# N = 25 は10秒程度で解決
# N = 25

# N=30 は4分弱かかる
N = 30

# これ以上の値での結果は未確認

# 絞り込んだデータを使って、初期化と散布図表示

data = df[:N]
data_np = np.array(data)
plt.figure(figsize=(6,6))
plt.scatter(data_np[:,0], data_np[:,1])
plt.grid()
plt.show()
```

### コスト配列の作成
任意のノード間のコストの配列(c)を座標値から計算します。
下の例では便宜上ユークリッド距離にしていますが、この距離の値をGoogle APIなど利用して2点間の実距離に差し替えれば、そのままリアルな最適化問題として使えるかと思います。

```py3
# 距離関数の定義
# 便宜上ユークリッド距離にしています

def distance(i, j):
    return np.sqrt(np.sum((data_np[i,:] - data_np[j,:])** 2))

# コスト配列の定義
c = np.zeros((N,N))

for i in range(N):
    for j in range(N):
        c[i, j] = distance(i, j)
```

## CPLEX APIによる問題定義
ここから先は、CPLEX APIを利用した、問題定義となります。

### 環境の確認
最初に次のAPI呼出しで、環境を確認してみます。

```py3
# CPLEX環境の確認

from docplex.mp.environment import Environment
env = Environment()
env.print_information()
```

手元の環境では、次のような結果になりました。

```
* system is: Linux 64bit
* Python version 3.6.8, located at: /opt/conda/envs/Python36/bin/python
* docplex is present, version is (2, 10, 150)
* CPLEX library is present, version is 12.9.0.1, located at: /opt/conda/envs/Python36/lib/python3.6/site-packages
* pandas is present, version is 0.24.1
```

### モデルの定義
最初にモデルインスタンスの定義を行います。TSPは「Traveling Salesman Problem」の略です。

```py3
# モデルの宣言

from docplex.mp.model import Model
mdl = Model(name="TSP")
```

### モデルパラメータの設定
モデルに対しては必要に応じてパラメータ設定ができるようです。下記のサンプルではスレッド数(MAXは2のようです)と最大時間数を設定しています。
(他にどんなことができるかは現在調査中)

```py3
# モデルのパラメータ定義

# スレッド数: 2
mdl.parameters.threads = 2
# 最大時間数: 600秒
mdl.parameters.timelimit = 600
```
### 決定変数の宣言
次は決定変数の宣言です。
決定変数が行列型のxと配列型のuの二種類を宣言します。それぞれの意味・目的はコメントに記載しておきました。

```py3
# 決定変数の宣言

# x : 移動matrix
# i番目のノードからj番目のノードに移動する時 x_ij = 1
# それ以外の場合 x_ij = 0
x = mdl.integer_var_matrix(N, N)

# u: 順序変数
# 0番目のノードの次に移動するノードがiの場合
# u[i] = 1
# その次に移動するノードがjの場合
# u[j] = 2
# .. のように定義する
u = mdl.integer_var_list(N)
```

### 制約の定義
次に制約の定義を行います。
元になる数式は、前の節に記載しておきましたので、そちらを参照して下さい。


```py3
# u(順序変数)の制約
# 最初のノードは0番目 : u[0] = 0
# それ以外の順序変数は1以上 N-1以下

mdl.add_constraint(u[0] == 0)
for i in range(1, N):
    mdl.add_constraint(u[i] >= 1)
    mdl.add_constraint(u[i] <= N-1)

# x(移動matrix)の制約
# 自分から自分への移動はないので、u_ii = 0
# それ以外の場合は u_ij は0か1

for i in range(N):
    for j in range(N):
        if i == j:
            mdl.add_constraint(x[i, j] == 0)
        else:
            mdl.add_constraint(x[i, j] >= 0)
            mdl.add_constraint(x[i, j] <= 1)

# x(移動matirx)に関する縦の制約

for i in range(N):
    mdl.add_constraint(mdl.sum(x[i, j] for j in range(N)) == 1)

# x(移動matirx)に関する横の制約

for j in range(N):
    mdl.add_constraint(mdl.sum(x[i, j] for i in range(N)) == 1)

# 部分路制約
# ループが全体で1つであるための条件
# 参考リンク http://satemochi.blog.fc2.com/blog-entry-24.html

for i in range(1, N):
    for j in range(1, N):
        mdl.add_constraint(u[i] - u[j] + N * x[i,j] <= N -1)
```

### 最適化関数の定義
最後に最適化関数の定義を行います。

```py3
# 総移動コスト

total_cost = mdl.sum(c[i, j] * x[i, j] for i in range(N) for j in range(N))

# 最適化の定義

mdl.minimize(total_cost)
```

### 設定の確認
これでモデルの記述は全部終わりました。
最後に``print_information``関数で、記述のサマリーを確認します。

```py3
# 制約設定の確認

mdl.print_information()
```

下記のような結果が返ってくるはずです。

```
Model: TSP
 - number of variables: 930
   - binary=0, integer=930, continuous=0
 - number of constraints: 2730
   - linear=2730
 - parameters: defaults
 - problem type is: MILP
```

元のノード数が30個の場合で、決定変数の数が930個、制約の数は2730個になっています。
これはとても人間の力では解けないので、ツールにお任せするしかなさそうです。

## 問題を解く
いよいよ``solve``関数を呼び出して、問題を解きます。
今回は、``log_output``のオプションも付けて、途中経過の表示も行ってみましょう。

```py3
# 問題を解く
# ログ出力をONにして、詳細情報も表示します

mdl.solve(log_output = True)
mdl.report()
```

下記のような結果が返ってくるはずです。

```
WARNING: Number of workers has been reduced to 2 to comply with platform limitations.
CPXPARAM_Read_DataCheck                          1
CPXPARAM_Threads                                 2
Tried aggregator 1 time.
MIP Presolve eliminated 1858 rows and 31 columns.
MIP Presolve modified 1624 coefficients.
Reduced MIP has 872 rows, 899 columns, and 4176 nonzeros.
Reduced MIP has 870 binaries, 29 generals, 0 SOSs, and 0 indicators.
Presolve time = 0.01 sec. (3.03 ticks)
Probing time = 0.01 sec. (2.90 ticks)
Tried aggregator 1 time.
Reduced MIP has 872 rows, 899 columns, and 4176 nonzeros.
Reduced MIP has 870 binaries, 29 generals, 0 SOSs, and 0 indicators.
Presolve time = 0.01 sec. (3.18 ticks)
Probing time = 0.01 sec. (2.89 ticks)
Clique table members: 466.
MIP emphasis: balance optimality and feasibility.
MIP search method: dynamic search.
Parallel mode: deterministic, using up to 2 threads.
Root relaxation solution time = 0.00 sec. (1.41 ticks)

        Nodes                                         Cuts/
   Node  Left     Objective  IInf  Best Integer    Best Bound    ItCnt     Gap

*     0+    0                        86255.7609        0.0000           100.00%
      0     0    18184.3571    62    86255.7609    18184.3571       51   78.92%
      0     0    20769.9760    72    86255.7609      Cuts: 93       94   75.92%
      0     0    20848.4887    72    86255.7609      Cuts: 45      126   75.83%
      0     0    20848.4887    72    86255.7609      Cuts: 43      141   75.83%
      0     0    20848.4887    71    86255.7609      Cuts: 15      152   75.83%
*     0+    0                        41843.4590    20848.4887            50.18%
      0     2    20848.4887    71    41843.4590    20848.4887      152   50.18%
Elapsed time = 0.18 sec. (85.35 ticks, tree = 0.02 MB, solutions = 2)
*   141   109      integral     0    32550.0557    20852.8478     1158   35.94%
*   239+  162                        32483.1076    20859.3823            35.78%
*   246   187      integral     0    32154.2342    20859.3823     1971   35.13%

(途中略)

325448 20165        cutoff          25354.1056    25157.6747  2674217    0.77%
 334254 15022        cutoff          25354.1056    25192.3348  2742881    0.64%
 343461  8762        cutoff          25354.1056    25244.3208  2809424    0.43%
Elapsed time = 221.50 sec. (173320.15 ticks, tree = 7.03 MB, solutions = 19)

GUB cover cuts applied:  22
Clique cuts applied:  12
Cover cuts applied:  22
Implied bound cuts applied:  15
Flow cuts applied:  1
Mixed integer rounding cuts applied:  102
Zero-half cuts applied:  41
Multi commodity flow cuts applied:  1
Lift and project cuts applied:  21
Gomory fractional cuts applied:  22

Root node processing (before b&c):
  Real time             =    0.19 sec. (85.04 ticks)
Parallel b&c, 2 threads:
  Real time             =  225.75 sec. (176705.75 ticks)
  Sync time (average)   =    3.91 sec.
  Wait time (average)   =    0.12 sec.
                          ------------
Total (root+branch&cut) =  225.94 sec. (176790.79 ticks)
* model TSP solved with objective = 25354.106
```

### 詳細の確認
``get_solve_details``関数で、解の詳細を確認してみます。

```py3
# 詳細情報の表示

print(mdl.get_solve_details())
```

こんな結果が返ってくるはずです。

```
status  = integer optimal, tolerance
time    = 225.95 s.
problem = MILP
gap     = 0.00996211%
```

## 結果のグラフ表示
最後に見つかった解を散布図上に表示してみましょう。

```
# 結果の取得

indexes = [u[i].solution_value for i in range(N)]
matrix = [[x[i, j].solution_value for j in range(N)] for i in range(N)]
ar = np.array(matrix)

# 結果の散布図上での表示

data_np = np.array(data)
plt.figure(figsize=(6,6))
plt.scatter(data_np[:,0], data_np[:,1])
for i in range(N):
    for j in range(N):
        if ar[i, j] == 1:
            plt.plot(data_np[[i,j],0], data_np[[i,j],1], c='b')
plt.grid()            
plt.show()
```

下記のようなグラフが表示されれば成功です。お疲れ様でした。

![スクリーンショット 2019-08-16 19.40.06.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/206130/2100bdfd-7f83-8418-64c4-8f2e756abc74.png)

# 追記 Nをもっと大きな数にすると
今までの記事で、Nをもっと大きくすると、どうなるのか知りたいと思います。
私も知りたいのでやってみました。
驚いたことに、N=35とN=40の場合は、N=30よりもっと短い時間で終わりました。
(N=35だと4秒、N=40だと42秒)
N=48のときは、全然終わりそうにない感じだったので、記事で紹介した時間制限を1800秒に設定しました。

※備考 Decisin OptimizerのNotebookを利用する場合、21CUH (1時間で21ユニット)の課金を使います。
Liteアカウントで利用可能な課金額は1か月で50ユニットですので、この処理を何度かやるとそれだけで課金を使いきることになります。Liteアカウントで試される場合はこの点にご注意下さい。

設定のためのコードを抜き出すと下記のとおりです。

```py3
# モデルのパラメータ定義

# スレッド数: 2
mdl.parameters.threads = 2
# 最大時間数: 1800秒
mdl.parameters.timelimit = 1800
```

この場合の結果を、以下に記載します。

### solve関数実行中のログ

```
CPXPARAM_Read_DataCheck                          1
CPXPARAM_Threads                                 2
CPXPARAM_TimeLimit                               1800
Tried aggregator 1 time.
MIP Presolve eliminated 4702 rows and 49 columns.
MIP Presolve modified 4324 coefficients.
Reduced MIP has 2258 rows, 2303 columns, and 10998 nonzeros.
Reduced MIP has 2256 binaries, 47 generals, 0 SOSs, and 0 indicators.
Presolve time = 0.01 sec. (7.30 ticks)
Probing time = 0.02 sec. (4.13 ticks)
Tried aggregator 1 time.
Reduced MIP has 2258 rows, 2303 columns, and 10998 nonzeros.
Reduced MIP has 2256 binaries, 47 generals, 0 SOSs, and 0 indicators.
Presolve time = 0.01 sec. (7.74 ticks)
Probing time = 0.01 sec. (3.58 ticks)
Clique table members: 1149.
MIP emphasis: balance optimality and feasibility.
MIP search method: dynamic search.
Parallel mode: deterministic, using up to 2 threads.
Root relaxation solution time = 0.00 sec. (3.68 ticks)

        Nodes                                         Cuts/
   Node  Left     Objective  IInf  Best Integer    Best Bound    ItCnt     Gap

      0     0    26762.8312    97                  26762.8312       94         
      0     0    31288.3138   122                   Cuts: 190      175         
      0     0    31520.2607   121                    Cuts: 84      229         
      0     0    31520.2607   116                    Cuts: 40      266         
      0     0    31520.2607   114                    Cuts: 32      299         
      0     0    31520.2607   117                    Cuts: 17      313         
*     0+    0                        44059.9234    31520.2607            28.46%
*     0+    0                        42371.0430    31520.2607            25.61%
      0     2    31520.2607   117    42371.0430    31520.2607      313   25.61%
Elapsed time = 1.34 sec. (1352.12 ticks, tree = 0.02 MB, solutions = 2)
    120   109    32236.4597    44    42371.0430    31520.2607     1982   25.61%
    240   222    34976.1678   105    42371.0430    31520.2607     3650   25.61%
    500   414    31536.0224   104    42371.0430    31535.0426     6334   25.57%
    783   701    41326.1454    47    42371.0430    31535.0426    10594   25.57%
*   873+  528                        41411.3945    31941.8932            22.87%
*   873+  352                        41082.6911    31941.8932            22.25%
    873   353    31941.8932   116    41082.6911    31941.8932    12784   22.25%
    875   355    31961.3763   126    41082.6911    31962.3013    12809   22.20%
    898   366    32044.5635   122    41082.6911    31963.5683    12886   22.20%
*   900+  226                        41075.0765    31963.5683            22.18%
*   908+  147                        40864.7463    31963.5683            21.78%
*   919+   94                        39048.8864    31963.5683            18.14%
    922   141    32138.2066   121    39048.8864    31963.5683    12966   18.14%
    940   110    32343.3044   125    39048.8864    31969.4607    13043   18.13%
   1152   219    34356.5554    97    39048.8864    31969.4607    14091   18.13%
Elapsed time = 11.57 sec. (8748.51 tick

(途中略)

396345 244626    32853.2816   122    33701.2618    32765.2102  5661464    2.78%
 396975 244777    33546.2942    86    33701.2618    32766.1908  5667055    2.77%
 397710 245350    33283.6134   119    33701.2618    32766.6977  5683128    2.77%
 398375 245730    33650.2717   115    33701.2618    32767.1283  5695264    2.77%
 399332 246349    33077.0530   142    33701.2618    32767.5542  5714250    2.77%
 400233 246907    33018.6804   143    33701.2618    32768.2769  5727993    2.77%
 401141 247321    32789.1947   103    33701.2618    32768.7445  5742685    2.77%
 402336 247850    33147.3776   102    33701.2618    32769.1675  5755769    2.77%
Elapsed time = 1782.07 sec. (1378944.05 ticks, tree = 306.56 MB, solutions = 24)
 403351 248512    33376.6836    96    33701.2618    32769.7014  5777789    2.76%
 404310 249171    33303.6713   101    33701.2618    32770.4043  5800246    2.76%
 405129 249660    33256.2355   120    33701.2618    32771.4467  5818181    2.76%

GUB cover cuts applied:  385
Clique cuts applied:  20
Cover cuts applied:  108
Implied bound cuts applied:  19
Flow cuts applied:  3
Mixed integer rounding cuts applied:  1147
Zero-half cuts applied:  38
Lift and project cuts applied:  4
Gomory fractional cuts applied:  31

Root node processing (before b&c):
  Real time             =    1.34 sec. (1351.27 ticks)
Parallel b&c, 2 threads:
  Real time             = 1798.67 sec. (1391470.14 ticks)
  Sync time (average)   =   61.04 sec.
  Wait time (average)   =    0.11 sec.
                          ------------
Total (root+branch&cut) = 1800.01 sec. (1392821.41 ticks)
* model TSP solved with objective = 33701.262
```

### 問題の散布図表示

![スクリーンショット 2019-08-17 13.37.33.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/206130/ede41eb3-d21a-ca05-bfd5-0bf30bbfbac0.png)

### 最終的な解のグラフ表示

![スクリーンショット 2019-08-17 13.37.40.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/206130/30ffc866-172f-17fd-0187-e01c1d35fd91.png)