首都圏の主要駅１日で全部回れるか（競プロで得た技術で遊ぼう～その２～）

#はじめに
どのように回ったら最も早く回れるかというのを計算できるとうれしいことがたくさんあると思います。たくさん回る箇所がある場合どう回るのが最も早いかをテーマにしたものに有名な巡回セールスマン問題があります。今日はそれを紹介しつつ遊んでいきたいと思います。

#巡回セールスマン問題
wikipediaには「都市の集合と各2都市間の移動コスト（たとえば距離）が与えられたとき、全ての都市をちょうど一度ずつ巡り出発地に戻る巡回路のうちで総移動コストが最小のものを求める（セールスマンが所定の複数の都市を1回だけ巡回する場合の最短経路を求める）組合せ最適化問題である。」とあります。
 これを求めるに単純なものでは、各地点を結ぶコストを求め、たどる順番を全探索し、それぞれの総移動コストを求め最小となるものを出力すればいいことが分かります。しかし、この計算量はたどる地点数をnとするとn！通りのたどり方があるためO(n!)となり、現実的に解けるのは、n≤10とかそこらになります。これを少し改善したのがbitDPといわれる手法で計算量はO(n<sup>2</sup>2<sup>n</sup>)となります。

#bitDPについて
詳しくはこちらが分かりやすいです→https://algo-logic.info/bit-dp/
今までに通った地点の集合、今いる地点におけるコストの最小を記録していく感じです。

dp[今までに通った地点の集合][今いる地点] = 最小コスト

みたいな感じです。

#問題設定
東京駅を出発し、首都圏の主要駅を巡ったのち東京駅に帰る最も早いルートは？関東の主要駅はhttps://www.metro.tokyo.lg.jp/tosei/hodohappyo/press/2017/01/20/07_01.html
より19駅を採用しました。
（東京、品川、恵比寿、渋谷、新宿、池袋、巣鴨、上野、秋葉原、船橋、大宮、浦和、赤羽、桜木町、吉祥寺、三鷹、国分寺、八王子、横浜）です。なお、公共交通機関を使って移動することとします。

#方針
###①各駅間のコスト決定
まずは各駅間での所要時間を求める必要があります。そこで今回用いたのがウェブスクレイピングという手法です。今回はyahoo乗り換え案内の方からデータを受け取ることとしました。詳しくはこちらの記事(https://selenium-rpa.net/selenium_auto_project/546/) に書いてあります。
###②最適化bitDP
コストを決めた後はbitDPで総コストが最小になるルートを求めればよいです。どのルートをたどったかについては経路復元の考え方を用いて、DPの更新と同時にどこから更新されたかを記録していけばよいです（DPの経路復元についてhttps://qiita.com/drken/items/0c7bab0384438f285f93 ).

#コード
①と②でそれぞれ結構時間がかかるので、分けて実装しました。
###①各駅間のコスト決定
エクセルファイルにはあらかじめたどる都市をいれておきます。
![cost.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/864771/e89adabe-1cbb-3c21-8f9a-ae02ed5974e1.jpeg)
（変更可能ですが今回は高速バス、新幹線は等は利用しないこととしました。）

```cost.py
import time
from selenium import webdriver
from selenium.webdriver.support.select import Select
import xlwings as xw

###各自で設定するところ(chromedriver.exeの場所,各地点間コストを記録するエクセルファイル*あらかじめ作製)
driver = webdriver.Chrome("C:/Users/user/Downloads/chromedriver_win32/chromedriver.exe") 
book_path = "C:/Users/user/Desktop/gamepython/timetable.xlsx"
###

wb1 = xw.Book(book_path)
ws1 = wb1.sheets[0]
dic = dict()
invdic = dict()
ii = 2
while ws1.cells(ii,1).value is not None:
    city = ws1.cells(ii,1).value
    dic[city] = ii-2
    invdic[ii-2] = city
    ii += 1
cityN = len(dic)
for i in range(cityN):
    for j in range(cityN):
        if i == j:
            ws1.cells(i+2,j+2).value = 0
            continue
        driver.get('https://transit.yahoo.co.jp/') #乗り換え案内
        if len(driver.find_elements_by_xpath('//a[contains(text(),"オフ")]')) == 1:
            driver.find_element_by_xpath('//a[contains(text(),"オフ")]').click()
        start_st = invdic[i]
        driver.find_element_by_xpath('//input[@id="sfrom"]').send_keys(start_st)
        end_st = invdic[j]
        driver.find_element_by_xpath('//input[@id="sto"]').send_keys(end_st)
        # 不変項目
        driver.find_element_by_xpath('//input[@id="tsAvr"]').click()
        select = Select(driver.find_element_by_xpath('//select[@name="ticket"]'))
        select.select_by_index(0)
        select = Select(driver.find_element_by_xpath('//select[@name="expkind"]'))
        select.select_by_index(0)
        select = Select(driver.find_element_by_xpath('//select[@name="s"]'))
        select.select_by_index(0)
        # 交通手段
        element = driver.find_element_by_xpath('//input[@id="air"]')
        if element.is_selected():
            element.click()
        element = driver.find_element_by_xpath('//input[@id="sexp"]')
        if element.is_selected():
            element.click()
        element = driver.find_element_by_xpath('//input[@id="exp"]')
        if element.is_selected():
            element.click()
        element = driver.find_element_by_xpath('//input[@id="hbus"]')
        if element.is_selected():
            element.click()
        element = driver.find_element_by_xpath('//input[@id="bus"]')
        if element.is_selected():
            element.click()
        element = driver.find_element_by_xpath('//input[@id="fer"]')
        if element.is_selected():
            element.click()
        driver.find_element_by_xpath('//input[@id="searchModuleSubmit"]').click()
        #time.sleep(1)    #Microsoft Edge
        if driver.current_url.find('https://transit.yahoo.co.jp/search/result') == 0:
                s = driver.find_element_by_xpath('//ul[@id="rsltlst"]/li[1]/dl/dd/ul/li[1]').text
                h = 0
                m = 0
                num = ''
                for k in s:
                    if k == '[' or k == '!' or k == ']':
                        continue
                    if k.isnumeric():
                        num += k
                        continue
                    if k == '時':
                        h = int(num)
                        num = ''
                    if k == '間':
                        continue
                    if k == '分':
                        m = int(num)
                mi = 60*h+m
                ws1.cells(i+2,j+2).value = mi
        else:
            ws1.cells(i+2, j+2).value = '検索エラー'
```
実行するとエクセルファイルに
![cost_result.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/864771/da518ed1-547e-db6b-5798-bf736ba154f0.jpeg)
と書き込まれます。
###②最適化bitDP
次にこのエクセルファイルを読み込んでbitDPを用いて最小コストを求めていきます。

```bitDP.py
import xlwings as xw
book_path = "C:/Users/user/Desktop/gamepython/timetable.xlsx"
wb1 = xw.Book(book_path)
ws1 = wb1.sheets[0]
dic = dict()
invdic = dict()
ii = 0
while ws1.cells(ii+2,1).value is not None:
    city = ws1.cells(ii+2,1).value
    dic[city] = ii
    invdic[ii] = city
    ii += 1
N = ii #都市の個数
lscost = [[0]*(N) for i in range(N)]
for i in range(N):
    for j in range(N):
        lscost[i][j] = int(ws1.cells(i+2,j+2).value)

dp = [[float('INF')]*N for i in range(2**N)]
comefrom = [[0]*N for i in range(2**N)]
dp[0][0] = 0
for i in range(2**N):
    for j in range(N):
        for k in range(N):
            if dp[i][k]+lscost[k][j] < dp[i|1<<j][j]:
                dp[i|1<<j][j] = dp[i][k]+lscost[k][j]
                comefrom[i|1<<j][j] = (i,k)
print('cost =',dp[-1][0],'min')
now = (-1,0)
root = [now]
while now != (0,0):
    a,b = now
    now = comefrom[a][b]
    root.append(now)
root.reverse()
ansls = []
for a,b in root:
    ansls.append(invdic[b])
print(*ansls,sep=' ')
```

実行結果

```result
cost = 378 min
東京 横浜 桜木町 品川 恵比寿 渋谷 新宿 吉祥寺 三鷹 国分寺 八王子 浦和 大宮 赤羽 池袋 巣鴨 上野 秋葉原 船橋 東京
```
というわけで６時間ちょっとで回ることができることが分かり、たどる順番も求めることができました。

#最後に
もっと地点数が多くなるとbitDPでは解くことが難しくなるので、数理モデルを用いたものや遺伝的アルゴリズムで最適に近いルートを割り出す手法が使われているようです（https://qiita.com/canonno/items/29eb369d32c030fea367 ）。巡回セールスマン問題は奥が深いなと思いました。

#参考記事一覧
https://ja.wikipedia.org/wiki/%E5%B7%A1%E5%9B%9E%E3%82%BB%E3%83%BC%E3%83%AB%E3%82%B9%E3%83%9E%E3%83%B3%E5%95%8F%E9%A1%8C

https://algo-logic.info/bit-dp/

https://www.metro.tokyo.lg.jp/tosei/hodohappyo/press/2017/01/20/07_01.html　

https://selenium-rpa.net/selenium_auto_project/546/

https://qiita.com/drken/items/0c7bab0384438f285f93

https://qiita.com/canonno/items/29eb369d32c030fea367