巡回セールスマン問題としての「道の駅スタンプラリー」

#地点数128の巡回セールスマン問題
北海道にある道の駅を巡る最短ルートを探してみました．いわゆる巡回セールスマン問題．
Python使って解かせてみたけど，

* PuLP → 計算終わらず
* OR-Tools → 非現実的なルート選択
* 深層強化学習 → 学習が進まずお手上げ

という結果になったので，その失敗談を残しておきます:sob:．
128を分割して幾つかのグループに分けてPuLPを適用することで，一応は`それらしいルート`を得ましたけど...悔しいです！

##北海道 道の駅スタンプラリーを最速で制覇するには

　「北海道はでっかいどう！」よく耳にするフレーズです．私の息子も幼少時よく大きな声で言ってました．どれだけ北海道がデッカいのか，それを身をもって体験できるのが「北海道 道の駅スタンプラリー」です．私は，北海道に来る前は東北に住んでまして，東北 道の駅スタンプラリーを制覇した経験があります．その自信もあって，北海道でも制覇してやる！と意気込んで参加しましたが，結果は….はい，無理でした（涙）．北の大地は広すぎました．そこで考えました．神様にお願いしても北海道のサイズは変わりません．せめてラリーの走行距離を短くできないかしら．これは…最短ルートを見つけるしかないっ！
　前置きが長くなりましたが，このような理由から「北海道内にある全ての道の駅[^1]を巡る最短ルート」を見つけ出すことに取り組みました．そして，なんとか得られたルートを以下に示します．ここでは「なとわ・えさん」から出発した場合を示していますが，巡回路（一筆書きで書けるルート）ですから，どこの道の駅をスタートにしても，このルートに沿って走れば最短距離でスタートした駅まで戻ってこれます．総移動距離は3917.4 kmです．

| Start↓ | つづき1 | つづき2 | つづき3 | つづき4 |
|:------------:|:------------:|:------------:|:------------:|:------------:|
| なとわ・えさん  | 自然体感しむかっぷ  | 流氷街道網走  | てしお  | 北欧の風　道の駅とうべつ  |
| 縄文ロマン　南かやべ  | 南ふらの  | メルヘンの丘めまんべつ  | なかがわ  | 望羊中山  |
| しかべ間歇泉公園  | しかおい  |  ノンキーランド　ひがしもこと | えんべつ富士見  | 名水の郷きょうごく  |
| つど〜る・プラザ・さわら  | うりまく  |  ぐるっとパノラマ美幌峠 | ☆ロマン街道しょさんべつ  | あかいがわ  |
| YOU・遊・もり  | ピア21しほろ  | 摩周温泉  | ほっと♡はぼろ  | スペース・アップルよいち  |
| くろまつない  | しほろ温泉  | あいおい  |  風Wとままえ | オスコイ！かもえない  |
| らんこし・ふるさとの丘  | かみしほろ  |  オーロラタウン93りくべつ |  るもい |  いわない |
| ニセコビュープラザ  | 足寄湖  |  おんねゆ温泉 | おびら鰊番屋  | シェルプラザ・港  |
| 真狩フラワーセンター  | あしょろ銀河ホール21  |  サロマ湖 | 森と湖の里ほろかない  |  みなとま〜れ寿都 |
| 230ルスツ  |  ステラ★ほんべつ | 愛ランド湧別  | 絵本の里けんぶち  |  よってけ！島牧 |
| とうや湖  | ガーデンスパ十勝川温泉  | かみゆうべつ温泉チューリップの湯  |  とうま |  てっくいランド大成 |
| とようら  | おとふけ  |  遠軽　森のオホーツク |  ひがしかわ「道草館」 | ルート229元和台  |
| あぷた  | なかさつない  | まるせっぷ  | びえい「丘のくら」 | あっさぶ  |
| みたら室蘭  | さらべつ  | しらたき  |  びえい「白金ビルケ」 | 江差  |
| だて歴史の杜  | 忠類  | 香りの里たきのうえ  | あさひかわ  | 上ノ国もんじゅ  |
| そうべつ情報館i   | コスモール大樹  | オホーツク紋別  | ライスランドふかがわ  |  北前船　松前 |
| フォーレスト276大滝  |  うらほろ | おこっぺ  |  鐘のなるまち・ちっぷべつ | 横綱の里ふくしま  |
| サーモンパーク千歳  | しらぬか恋問  | おうむ  | サンフラワー北竜  | しりうち  |
| 花ロードえにわ  | 阿寒丹頂の里  | にしおこっぺ花夢  |  田園の里うりゅう |  みそぎの郷　きこない |
| マオいの丘公園  | 厚岸グルメパーク  |  もち米の里☆なよろ | たきかわ  | なないろ・ななえ  |
| ウトナイ湖  | スワン44ねむろ  | びふか  | スタープラザ芦別  |  なとわ・えさん |
| サラブレッドロード新冠  |  おだいとう |  おといねっぷ |  うたしないチロルの湯 | **Goal !** :smiley: |
| みついし  | 知床・らうす  | ピンネシリ  | つるぬま  | - |
| むかわ四季の館  | うとろ・シリエトク  | マリーンアイランド岡島  | ハウスヤルビ奈井江  | - |
| あびらD51ステーション  | しゃり  | 北オホーツクはまとんべつ  |  三笠 | - |
| 夕張メロード  | パパスランドさっつる  |  さるふつ公園 | しんしのつ  | - |
| 樹海ロード日高  | はなやか（葉菜野花）小清水  | わっかない  | 石狩「あいろーど厚田」  | - |
| **つづき1へ ↗** | **つづき2へ ↗** | **つづき3へ ↗** | **つづき4へ ↗** | - |

Googleマップの著作権の関係から，ここに最短ルートを示した北海道全体のマップは載せられませんので，数が多くなりますが分割して示します．なお，有料道路は避けています．

https://goo.gl/maps/MXn6ZJU98GoNHkeT8

https://goo.gl/maps/4SiQVogLUEo4HxKe9

https://goo.gl/maps/o2YbPyLkNJWBRaR17

https://goo.gl/maps/9qVJShVhWFBpSMHC9

https://goo.gl/maps/aHa8kvCoRNjfmpxB9

https://goo.gl/maps/1w3DD1XbLAjSdtvN9

:point_up:**「知床・らうす」~「うとろ・シリエスク」間は要注意です！**
最短距離である国道334号線経由は4月下旬~10月しか通れません．

https://goo.gl/maps/bN2ppJ6FKYpdk5G79

https://goo.gl/maps/B244TjeYbehrD8bR8

https://goo.gl/maps/x1Ybc7fjMexCDqfw7

https://goo.gl/maps/sS1aCdQhxMpYtwYL8

https://goo.gl/maps/ehdCgxFR6XBqc1DG6

https://goo.gl/maps/cmPqJgS6YqPJgEdX7

https://goo.gl/maps/J988bLFCuZ7HFzyq6

https://goo.gl/maps/jzrYJnLVoyYbcuURA

https://goo.gl/maps/XN9NmBuDfxcLCb5C6

ところどころ，来た道を引き返していますが，「最短距離」にこだわると致し方ない結果でしょう．

[^1]: 令和2年7月末日時点の地図データを基にしているため，道の駅の数は128です．

##このルートを見つけるまでの試行錯誤

　今回の「こだわり」は，Google mapを使って得た正確な距離情報を基に，128という多くの訪問地点を対象に解析したことです．10~30の訪問地点間を単純に直線でつないで最短巡回ルートを求めた例[^2]はインターネット上で散見されますが，実際の道に沿った距離を使って北海道の全ての道の駅を対象にした巡回セールスマン問題を解きたかったのです．

[^2]:  たとえば https://qiita.com/Kanahiro/items/08df9b18d471c0ba678e

　まず，Google mapを使って各道の駅間の距離を調べて距離行列をつくります．ただ，これを素直に作るとなると ${}_{128}\mathrm{C}_2 $ 通りの距離を調べる必要があるため，めんどくさい:stuck_out_tongue:．そこで閾値[^3]を設けて，明らかにその値以上の距離にある場合は1000kmという大きな値を入れました．

[^3]:  最寄り「道の駅」同士の距離が大きい釧路・根室エリアを考慮して，80kmに設定しました.

後日談：以下の記事を読んで，Webスクレイピングを使えばGoogleマップから距離情報を自動で収集できるのかな～？と思ったりもして，まだ試していませんが，いずれチャレンジしようと思っています．

https://qiita.com/tomato1997/items/bb828396e8abfd767def#%E6%96%B9%E9%87%9D

###1. Pythonライブラリ PuLP

　PythonのライブラリであるPuLPを使って最短ルートを求める計算プログラムを作ります．

参考にしたページはこちら

https://qiita.com/panchovie/items/6509fb54e3d53f4766aa

このページ下部に訪問地点数と求解に要した時間の関係図が載っており，それを見ると地点数128の場合は途方もない時間がかかるようです．試しに，ノートPCのCPUで計算してみると案の定，1ヶ月かかっても計算が終わらない:sweat_smile:．そこで，北海道全体を幾つかのサブエリアに分けて，各サブエリアの中でPuLPを使って最短距離を見つけて，最後にそれらを連結する！という工夫をしました．この記事を書いている時に偶然見つけた論文 [Alhanjouri (2018)](http://ijiset.com/vol5/v5s6/IJISET_V5_I06_04.pdf) も同じような考え方してますね．

ということで，地道に作った距離行列を入力データとして，コンピュータにせっせと計算してもらった結果，上で示したような巡回ルートをゲットできました．

###2. Pythonライブラリ OR-Tools

　有名な数理最適化ライブラリであるOR-Toolsを使ってもトライしてみました．

参考にしたページはこちら

https://qiita.com/SaitoTsutomu/items/ab9d657c49879df69928

このアルゴリズムは近似解法です．`limit_time=10`として計算しました．
結果は...うまく行きませんでした．どうも十勝から亀田半島に飛びたがるのです．つまり非現実的なルートを選択します．OR-Toolsでは日高ロード沿いにある「むかわ」「新冠」「みついし」をUターンしてくるルートを捻り出すことができない模様．距離行列を作る際に，全ての距離を正確に入力しないとダメなのかもしれません（「遠すぎる場合は一様に1000km」ではダメということ）．

###3. 深層強化学習

　深層強化学習を使った求解にもトライしてみました．

参考にしたページはこちら．

https://qiita.com/panchovie/items/86323946cceca6695e91

関連する書籍も出ています → [現場で使える！Python深層強化学習入門](https://www.shoeisha.co.jp/book/detail/9784798159928 "6.2節 巡回セールスマン問題")
この書籍のwebサイトから，深層強化学習の実装例がダウンロードできます[^4]．アルゴリズムに関しては，どうやらarXivの論文 [Belloら (2017)](https://arxiv.org/abs/1611.09940 "原著論文") を参考にしているようです．Pointer Networkを利用した方策勾配ベース(Actor-Critic)のアルゴリズムとのこと．

[^4]: フレームワークがTensorflowの1系のため，初めてTensorflowを触る私にはコードの改変が大変でした．特に，for構文がサポートされていないため，二重ループをtf.while_loopで表現する方法を考えるのに苦労しました．

　私はディープラーニングの初学者かつ組合せ最適化問題の非専門家（いわゆる素人さん）のため，「巡回セールスマン問題に深層強化学習を適用するのは得策ではない」と仰る専門家の先生を尻目に，あえて興味本位でトライしてみました．

https://note.com/mikiokubo/n/n94853a581bb8

その結果は...やはり訪問地点数128ではダメですね（涙）．
下記の通り，（効率的に解を算出できるかは別として）地点数42だと正解にたどり着くことを確認しましたが，さすがに100オーバーだとキビしいようです．

>地点数42の例題は[こちら](https://people.sc.fsu.edu/~jburkardt/datasets/tsp/tsp.html "DANTZIG42")
ページ中程にある"DANTZIG42"　距離行列が既に準備されています．

>訪問地点の分布状態はこんな感じ↓
![DANTZIG42.JPG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1455270/2db22946-f2eb-4d25-7cb1-e836a3b70777.jpeg)

>深層強化学習により得られた最短巡回ルート↓
![Best_route.JPG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1455270/e13790b3-f9c4-bd61-b94c-0e666129b69b.jpeg)

>上記の例題ページに "The minimal tour has length 699" と書いてあるので，おそらくこれで合っていると思います．

　Belloさん達の原著論文では「訪問地点数100でも上手くいったよ」的な書きっぷりだったので，その論文からハイパーパラメタを拝借してミニバッチ数128，学習係数0.0001の5000ステップ毎に0.96倍，重みの初期値は[-0.08, 0.08]の一様乱数，閾値1の勾配クリッピングの導入，オプティマイザをAdamに変更....などを採用して計算を走らせるも，学習が途中で立ち往生して進まない:weary:．

このゴールデンウィーク中ずっと「ああでもないこうでもない」といろいろ試してみましたが，力尽きて諦めました．
Pythonを勉強して1年，Pytorch歴3ヶ月，Tensorflowに関しては「何ですか？それ」レベルの私にとっては，この取り組み自体が自分の勉強のために良い教材になりました．

##あとがき

　得られた最短ルートのマップを見て，「こんなの当たり前じゃ~ん」と言う人がいるかもしれません．しかし，北海道の白地図上に道の駅を点で示したものを我々人間が見て，パッと最短ルートを見つけることは難しいです．いや，無理です（実際にやってみればわかります）．いわゆる天才ならできるのかな？

[北海道開発局](https://www.hkd.mlit.go.jp/ky/kn/dou_kei/ud49g70000008gn6.html "北海道開発局のHP")　←ここに載っています．