TSP LP (Linear Programming)：線形計画法による巡回セールスマン問題 / Subtour Elimination 部分巡回路除去

前回はDP（動的計画法）でTSP（巡回セールスマン問題）の厳密解を求めてみましたが、今回はLP（線形計画法）で厳密解を求めてみます。〜計画法と名前は似ていますが全く異なるアルゴリズムです。線形計画法は最大値や最小値を求める最適化計算のようで、大体は最適化ソルバーを使って計算するそうです。有償ソルバとしてIBMのCPLEXやGurobiなどありますが、今回はオープンソースのPulpを使ってみます（Python 3.6、Jupyter Notebook使用）。
追記：
IBMのCPLEXは「pip install cplex」でインストールできますが、無償版のためリミッターがついているらしく、ノード数を20以上にするとストップしてしまいました。
cvxpyというオープンソースのソルバーも使ってみました。cvxpyでもソルバーとしてCPLEXを選択できますが、やはりノード数を増やすとストップしてしまいます。
TSPの場合MIP（混合整数線形計画法）でなければいけないようで選択できるソルバーも限られているみたいです。その点PulpはデフォルトでもMIPに対応したCOINというソルバーを使っているようで（CPLEXにも変更可）問題なく解くことができました。


TSPの場合、経路合計距離が短いほどいいので最小化する式を考えなければいけません。
訪問先の都市（ノード）を

S = {0, 1, 2, 3, 4, 5}

とし二点間の経路を

x(i, j)

にしておきます。iもjも0〜5の値が入りi≠jとなります。
各経路には固有の距離があるので、その距離を

dist(i, j)

にしておきます（配列や辞書など使って事前に二点間の距離を記録しておきます）。
iとjの組み合わせは6*5=30通りありますが、最終的な経路として採用されるのは、例えば、

x(0, 1) → x(1, 2) → x(2, 3) → x(3, 4) → x(4, 5) → x(5, 0)

のような30通りあるうちの6つとなります。
そうすると、

x(i, j) = 1：経路として採用

x(i, j) = 0：経路として不採用

というように0と1で表し、それに距離を掛けて全体の経路長は以下のように表すことができます。

Σ dist(i, j) * x(i, j)

これを線形計画法で最小化していくようです。要は二点間(i, j)の距離に(i, j)が経路に含まれている場合は1を掛けて、含まれていない場合は0を掛けて全ての(i, j)の組み合わせを合算していき、最小となる組み合わせを探すというわけです。意外に簡単な式でTSPの特徴を表すことができます。

ただし、いくつか条件が必要になってきます。
例えば、0を起点に他の地点へ向かう組み合わせx(0, j)は５つあり、そのうちどれか一つが採用されることになります。採用された経路は1で残りは0となるので、

x(0, 1) + x(0, 2) + x(0, 3) + x(0, 4) + x(0, 5) = 1

という式が成り立ちます。
もしx(0, 1)とx(0, 2)が採用されてしまうと、0から1と2へ向かう分岐点が出来てしまうので矛盾してしまいます。よって、これらを足し合わせると必ず1になり、0から向かうことができる地点は1箇所しかないという制約になります。

他の地点に関しても同様で、

x(0, 1) + x(0, 2) + x(0, 3) + x(0, 4) + x(0, 5) = 1  (i≠j)

x(1, 0) + x(1, 2) + x(1, 3) + x(1, 4) + x(1, 5) = 1  (i≠j)

x(2, 0) + x(2, 1) + x(0, 3) + x(0, 4) + x(0, 5) = 1  (i≠j)

x(3, 0) + x(3, 1) + x(3, 2) + x(0, 4) + x(0, 5) = 1  (i≠j)

x(4, 0) + x(4, 1) + x(4, 2) + x(4, 1) + x(0, 5) = 1  (i≠j)

x(5, 0) + x(5, 1) + x(5, 2) + x(5, 3) + x(5, 4) = 1  (i≠j)

となり、for文を使えば、

for i in S:
    sum([x(i, j) for j in S if i != j]) == 1

となります。
上記はiからjへの移動に関する制約ですが、その反対であるjからiへの移動についても同じ制約があります。

x(1, 0) + x(2, 0) + x(3, 0) + x(4, 0) + x(5, 0) = 1  (i≠j)

x(0, 1) + x(2, 1) + x(3, 1) + x(4, 1) + x(5, 1) = 1  (i≠j)

x(0, 2) + x(1, 2) + x(3, 2) + x(4, 2) + x(5, 2) = 1  (i≠j)

x(0, 3) + x(1, 3) + x(2, 3) + x(4, 3) + x(5, 3) = 1  (i≠j)

x(0, 4) + x(1, 4) + x(2, 4) + x(3, 4) + x(5, 4) = 1  (i≠j)

x(0, 5) + x(1, 5) + x(2, 5) + x(3, 5) + x(4, 5) = 1  (i≠j)

同様にfor文で書けば、

for j in S:
    sum([x(i, j) for i in S if i != j]) == 1

となります。

二つの制約をiとjのマトリクス（以下）で考えれば、

------- x(0, 1) x(0, 2) x(0, 3) x(0, 4) x(0, 5)

x(1, 0) ------- x(1, 2) x(1, 3) x(1, 4) x(1, 5)

x(2, 0) x(2, 1) ------- x(2, 3) x(2, 4) x(2, 5)

x(3, 0) x(3, 1) x(3, 2) ------- x(3, 4) x(3, 5)

x(4, 0) x(4, 1) x(4, 2) x(4, 3) ------- x(4, 5)

x(5, 0) x(5, 1) x(5, 2) x(5, 3) x(5, 4) -------

横軸で見るか縦軸で見るかの違いになります。つまり横軸の中からは一つしか選べず、同様に縦軸からも一つしか選べないということになります。

仮にx(0, 1)を選べば、次を選ぶ場合、横軸内x(0, 2), x(0, 3), x(0, 4), x(0, 5)はどれも選べず、同様に縦軸内x(2, 1), x(3, 1), x(4, 1), x(5, 1)も選べないことになります。よってこのような二重の制約内で経路を探すことになります。

追記：
以下はノード数10の場合、どの経路が選べるかというシミュレータです（editor.p5js.orgで作成）。

グレー：選択可能な(i, j)の組み合わせ
黒　　：選択不可能
ピンク：選択した(i, j)

条件：
・選択した(i, j)と同じ横軸上の(i, j')は選択不可。
　例えば(1, 2)を選択したら(1, 0)〜(1, 9)は選択不可。
・選択した(i, j)と同じ縦軸上の(i', j)は選択不可。
　例えば(1, 2)を選択したら(0, 2)〜(9, 2)は選択不可。
・選択した(i, j)を反転した(j, i)は選択不可。
　例えば(1, 2)を選択したら(2, 1)は選択不可。
・二つ以上（かつ全ノード数未満）選択して部分巡回路が{i, j, k}となったとき(k, i)は選択不可。
　例えば(1, 2)、(2, 4)を選択したら(4, 1)は選択不可。また(1, 2)、(2, 4)、(4, 5)を選択したら(5, 1)は選択不可。

Subtour Elimination（部分巡回路除去）：

しかしながら先ほどの二つの条件だけだと、ひとつながりの経路ではなく、

0 → 2 → 0

3 → 1 → 4 → 5 → 3

のような複数の巡回路ができてしまう場合があり、まだ制約としては十分ではないようです。

そこでSubtour elimination（部分巡回路除去）という制約を追加しなければいけません。これにはいくつかのアルゴリズムがあるようで、今回はまずMTZ(Miller-Tucker-Zemlin)の方法を試してみることにしました。

以下が追加の制約式です。一見わかりにくいですが、(i, j)の二点間における制約です。

ui - uj <= N * (1 - x(i, j)) - 1  (i≠0, j≠0, i≠j)

Nはすべてのノード数：N=|{0, 1, 2, 3, 4, 5}|=6。uiとujは、x(i, j)の２点間経路における経路順を表す値です。x(i, j)はiからjへ向かうという意味なので、仮にiが巡回路上の3番目（ui=3）にある地点だとすれば、jはその次の4番目（uj=4）となり、

ui - uj = 3 - 4 = -1

になります。iは常にjの一つ手前なのでui-ujが-1になるのは当然ですが、一周回って最後のときではiが6番目でjが１番目なので、

ui - uj = 6 - 1 = 5

になってしまいます。

5 <= 6 * (1 - 1) -1 = -1

となるので右辺と左辺の計算が合いません。
右辺の(1 - x(i, j))の部分は、巡回路に(i, j)が含まれれば1-1=0で、含まれなければ1-0=1となるという意味で、その結果にノード総数を掛けて、

N * (1 - x(i, j)) - 1 = -1

になります。もしx(i, j)が巡回路になければN-1=5となり、-1か5の値をとるということがわかります。

この制約式にはi≠0, j≠0, i≠j（あるいは1から始まるのであればi≠1, j≠1, i≠j）という条件があるので、スタート地点0を含まない方の部分巡回路で検証してみます。

3 → 1 → 4 → 5 → 3

に対してuiを割り当てると

u3,  u1,  u4,  u5, u3

となり、それぞれの部分巡回路内での道順の値を入れると、地点3が1番目なのでu3=1、地点1が2番目なのでu1=2という感じで（多分最初を0番目としても構わない）、

u3=1, u1=2, u4=3, u5=4, u3=1

となります。まず3から1への移動において制約式に代入してみます。x(i, j)の部分は経路に含まれている場合は1で含まれていない場合は0ですが、この部分巡回路には含まれているので1となります。

u3 - u1 + 6 * x(3, 1) <= 6 - 1

1 - 2 + 6 * 1 <= 6 - 1

5 <= 5

となり問題ありません。

次に、1から4へ移動するx(1, 4)を制約式に代入してみます。

u1 - u4 + 6 * x(1, 4) <= 6 - 1

2 - 3 + 6 * 1 <= 6 - 1

5 <= 5

これも問題ありません。同様に4から5へ移動するx(4, 5)を代入しても問題ありません。
しかし、最後の5から3への移動x(5, 3)を代入してみると、

u5 - u3 + 6 * x(5, 3) <= 6 - 1

4 - 1 + 6 * 1 <= 6 - 1

9 <= 5（不適合）

最後の不等式の辻褄が合わないためこの経路については無効となり除外されてしまいます。

もし仮に正解の巡回路が以下だとします。

0 → 2 → 3 → 1 → 4 → 5 → 0

最後の5 → 0の部分を検証してみると（スタート地点0を1番目、地点5を6番目とすると）

u5 = 6, u0 = 1

u5 - u0 + 6 * x(5, 0) = 5 - 1 + 6 * 1 = 10

となって、左辺と右辺を比較してみると、

10 < 6 - 1

となって適合しませんが、制約式の条件がi≠0, j≠0, i≠jなのでスタート地点を含んでいる場合は検証する必要がありません。どうやらこれはスタート地点を除いた他の地点に関する制約のようです。uiが巡回路に沿って追いかけていき、いずれ終点（スタート地点でもある）に到着したときに必ず制約式に違反してしまいます。しかしながら、スタート地点を含む巡回路に適用しないということで、その結果、他の部分巡回路を排除するという感じでしょうか。
このMTZの部分巡回路除去の制約式は便利そうですが、ノード数が増えるほどN^2-Nのバイナリ計算量が必要となるので、せいぜい数十個程度のものにしか適用できないようです。


Pulpで実装：
最適化ソルバーPulpを使ってこのアルゴリズムを書いてみました（ここを参考）。
Pulpではまず最適化（今回は最小化）するProblemを用意します。

prob = LpProblem(name='TSP_LP', sense=LpMinimize)

二点間の距離はマトリクスで表すよりも、二点間を辞書のキーに登録して、値にその距離を割り当てておきます。

dist = dict({(i, j): (abs(XY[i] - XY[j])) for i in range(num) for j in range(num) if i != j})

N*(N-1)通りの組み合わせがあります。
そして、最適化する変数Variable: xを用意しますが、これも辞書型を利用して下限値0、上限値1のバイナリなので、

x = LpVariable.dicts('x', dist, 0, 1, LpBinary)

になります。なんとなくTensorflowの書き方似ています。

次に、最適化（最小化）する式を用意します。

prob += lpSum([dist[(i,j)] * x[(i,j)] for i, j in dist])

二点間の距離とその経路の有無（1か0）を掛けて合計したものが最小になればいいということです。Loss関数のような感じです。

そしてベースとなる制約式です。

for i in path:
    prob += lpSum([x[(i, j)] for j in path if i != j]) == 1
    
for j in path:
    prob += lpSum([x[(i, j)] for i in path if i != j]) == 1

iから向かう他の地点は1箇所、同様に他からiに向かう地点も１箇所という制約です。
さらに、部分巡回路除去のMTZの制約を加えると、

u = LpVariable.dicts('u', path, 0, num)

for i in path:
    for j in path:
        if i != j and (i != 0 and j != 0) and (i, j) in x:
            prob += u[i] - u[j] <= num * (1 - x[(i, j)]) - 1

となります。前述したようにこの制約はスタート地点0には適用しないので、i!=0, j!=0になっています。
あとは、これをソルバーに解かせるだけです。

prob.solve()

問題がなければ、

LpStatus[prob.status]

で、

Optimal

が出力されます。

x[(i, j)].valValueでどの(i, j)の組み合わせが1になっているかわかります。1になっている(i, j)が今回経路として選ばれたことになるので、それを数珠つなぎに経路用配列Tへ入れていきます。

T = [0]
while len(T) < num + 1:
    for (i, j) in x:
        if x[(i, j)].varValue == 1 and i == T[-1]:
            T.append(j)

tour_dist = sum([dist[(T[i - 1], T[i])] for i in range(1, len(T))])

print('tour         :', T)
print('tour distance:', tour_dist)
plot_path(T)

全体の距離はdist[(i, j)]を積算すればでてきます。
plot_path()は事前に用意しておいた描画ファンクションです。

ノード数20で計算させてみました。5秒ほどかかりましたが、ソルバーのおかげなのか前回の動的計画法（DP）よりは速いです。
Subtour Elimination（部分巡回路削除）の方法にいくつかあるようなので、他のアルゴリズムも試してみようと思います。

続き：
TSP LP Lazy Subtour Elimination Constraints：巡回セールスマン問題 　逐次組込制約による部分巡回路除去

関連：Traveling Salesman Problem:巡回セールスマン問題について（まとめ）

TSP DP（その２） bit DP / 巡回セールスマン問題 / 動的計画法

前回のTSP DP（動的計画法による巡回セールスマン問題）では、訪問先を表現するためにset()関数を用いましたが、今回はbit DPで試してみました。
bit演算させたほうが高速になるのと、メモ化ライブラリであるfunctools.lru_cacheも使えるようになります。

ビット演算でどこを訪れていて、どこを訪れていないか（訪問フラグ）は以下のように表すようです。
まず訪問先を0, 1, 2, 3とすれば4桁の二進数になります。Pythonでは左側に0bをつければ二進数となります。

S = 0b0001

0ビット目（一番右側）を1にしておいて、ここが出発点となります。
つぎに2に訪れたら（右から0番目、1番目、2番目）、

S = 0b0101

となり、すでに訪れた場所の状態がわかります。
次に1に訪れて、

S = 0b0111

つぎに3に訪れて、

S = 0b1111

最終的にすべてが1になれば完了というわけです。
ビット演算用の記号を使いますが、実際は整数としても認識するので、

S = 1 | 1 << 2 = 0b0101

こんな感じになります。

ということから二進数を用いて前回のTSPを書き直してみました。


functools.lru_cacheを使った場合：（Python 3.6、Jupyter Notebook使用）
メモ化ライブラリであるfunctools.lru_cacheを使うと以下のようになります。
ちなみにこれは合計最短経路のみの出力です。尚、TSP/DPのアルゴリズムについては前回の方が詳しく説明してあるかと。

from functools import lru_cache

@lru_cache(maxsize=None)
def TSP_lru(S, pos):        
    if S == all_visited:
        return dist[pos][0]
    return min([dist[pos][i] + TSP_lru(S | (1 << i), i) for i in range(num) if S & (1 << i) == 0])

# @lru_cache(maxsize=None)
# def TSP_lru(S, pos):  # same as the above
#     if S == all_visited:
#         return dist[pos][0]
#
#     d_min = float('inf')
#     for i in range(num):
#         if S & (1 << i) == 0:     # if not visited 
#             d = dist[pos][i] + TSP_lru(S | (1 << i), i)
#             if d < d_min:
#                 d_min = d
#   return d_min

# Solve by TSP with functools.lru_cache
%time print('Total Distance:', TSP_lru(1, 0))
print(TSP_lru.cache_info())  # show cache info
TSP_lru.cache_clear()        # clear cache otherwise same outputs as before

メモ化されていないファンクションの上に@マークのデコレータをつけるだけでメモ化されます。
maxsize=Noneでメモリ制限なしにしています。デコレータなしだとメモ化しないためかなり遅くなります。
ただし、前回のようなset関数やリストを引数にしている場合はこのlru_cacheが使えません。タプルを引数にするなら大丈夫だと思います。今回は訪問先の集合をビットで表したので使うことができるというわけです。おそらくビット演算にするだけでも、少しは高速になっているのではないでしょうか。
それから、一度runするとcacheに内容が保持されたままになるので、再度TSR_lru(1, 0)をrunすると、前回と同じ答えを出してしまいます。ノード数を変えて試してみたい場合は、毎回cache_clear()する必要があります。何度も同じ答えを呼び出すには高速ですが、キャッシュクリアしないかぎりは同じ答えを出し続けるので、他のプログラムで使うときには注意したほうがいいと思います。


memoize関数を使った場合：
自前で書いたmemoize関数をデコレータにすることもできるようです（こちらを参考）。
合計最短経路のみの出力。

def memoize(f):
    cache = {}
    def func(*args):
        if not args in cache:
            cache[args] = f(*args)
        return cache[args]
    return func

@memoize
def TSP_memo(S, pos):  # re-run this function to clear cache otherwise same outputs as before
    if S == all_visited:
        return dist[pos][0]
    return min([dist[pos][i] + TSP_memo(S | (1 << i), i) for i in range(num) if S & (1 << i) == 0])

# Solve by TSP with memoization decorator
%time print('Total Distance:', TSP_memo(1, 0))

これも先ほどのlru_cacheと同じようなものです。
この場合もキャッシュクリアしない限りは同じ答えを出すので、TSP_memo(1, 0)をrunさせる場合は、def TSP_memo():も一緒に（jupyter上のこのセルごと）runさせたほうがいいです（このキャッシュクリアの方法がわからない）。


リストを使ったDPメモ化：
ビット化した以外はほぼ前回と同じです。
最後にメモ化した内容から経路順を割り出します。

def TSP_BDP(S, pos):
    # if visited before
    if dp[S][pos] >= 0:
        return dp[S][pos]
    
    # if all visited  
    if S == all_visited:         
        dp[S][pos] = dist[pos][0]
        return dp[S][pos]
    
    # if not visited yet
    dp[S][pos] = min([dist[pos][i] + TSP_BDP(S | (1 << i), i) for i in range(num) if S & (1 << i) == 0])
    
    return dp[S][pos]

# Solve TSP by bit DP (15 nodes: 0.75 sec, 18 nodes: 8.42 sec, 20 nodes: 44 sec)
dp = [[-1] * num for _ in range(1 << num)]    # initialize memoization
%time print('Total Distance:', TSP_BDP(1, 0))

基本的に、

dp[現在訪問状態][現在地] = dist[現在地][移動先] + TSP_BDP(次の訪問状態, 移動先)

を毎回記録/計算しています。forループのi（移動先）は、未訪問の地点があるかどうか確かめてから、あればそこへ移動する場合の距離を計算しmin()によって最短経路を選択します。
具体的な値を入れてみると、

dp[0b1011][3] = dist[3][2] + TSP_BDP(0b1011 | (1 << 2), 2)

Sについて2ビット目が0となっているので（右から数えて0,1,2番目）、まだ訪れていない地点は2ビット目だけということになり、forループでi=2のときこの処理がなされます。
そうすると、

dp[0b1011][3] = dist[3][2] + TSP_BDP(0b1111, 2)

右辺においてS=0b1111（全て訪問済み）となり、あとは現在地3から2までの距離と残りTSP_BDP(0b1111, 2)による距離を足して左辺のdp[0b1011][3] へ記録することになります。
しかし、dp[0b1011][3] へ記録する前に、TSP_BDP(0b1111, 2)を計算しなければいけないので、TSP_BDP(0b1111, 2)は次の再帰ループに突入してif S == all_visited:の条件分岐により、

dp[0b1111][2] = dist[2][0]

になり距離テーブルdist[2][0]から値を得てdp[0b1111][2]に記録します。このとき再帰ループで一つ先に移動しているので現在地pos=2、移動先は0になっています。距離が記録されたのちdp[0b1111][2]（= dist[2][0]）がリターンされたのでTSP_BDP(0b1111, 2)に戻ります。
この結果、

dp[0b1011][3] = dist[3][2] + TSP_BDP(0b1111, 2)

この式の右辺側にあるTSP_BDP(0b1111, 2)は、

TSP_BDP(0b1111, 2) = dist[2][0]

とわかったので、それを代入すると

dp[0b1011][3] = dist[3][2] + dist[2][0]

つまり、

未訪問2で現在地3 = 3から2までの距離 + 2から0（ゴール）までの距離

となってゴールから二歩手前までの距離が求まります。

訪問フラグと距離の関係（4ビットの場合、posはその時の現在地）：
dp[ob1111][pos]（4ビット1の場合）ゴールから一歩手前までの距離
dp[0b1011][pos]（3ビット1の場合）ゴールから二歩手前までの距離
dp[0b0011][pos]（2ビット1の場合）ゴールから三歩手前までの距離
dp[0b0001][pos]（1ビット1の場合）ゴールからスタートまでの距離

ということから、この仕組みを理解して経路の生成を以下で行いました。


経路の生成：
上で記録した二次元配列dpを元に経路の順番を生成します。

def gen_dp_tour():
    tour = [0]
    S = 1
    while len(tour) < num:
        d_min = np.inf
        pos = None
        for i in range(num):
            if S & (1 << i) == 0:
                d = dist[tour[-1]][i] + dp[S | (1 << i)][i]
                if d < d_min:
                    d_min = d
                    pos = i
        S |= (1 << pos)
        tour.append(pos)
    return tour + [0]

# Generate Tour
Tour = gen_dp_tour()
print('Tour:', Tour)
plot_path(Tour)

基本的には最短経路を求めるコードと似ています。前のコードのforループ内のTSP_BDP()の部分をdp[][]に置き換えて距離を計算させるかわりに直接メモリから呼び出します。再帰ループの代わりにwhileループにしています。最初出発点にいる状態から全てを訪問するまで、min関数で得た最短距離と同時に現在地も記録しておき、それをtour配列に入れていきます。
先ほどの訪問フラグと距離の関係から、dp[0b0001]にはスタートからゴールまでの合計最短距離が記録されています。これを一つずつ分解していく感じです。

dp[0b0001][pos] = dist[pos][i] + dp[S | (s << i)][i]

というのは、

dp[現在訪問状態][現在地] = dist[現在地][移動先] + dp[次の訪問状態][移動先]

なので、初期状態は0ビット目にいるため0b0001になり現在地pos=0になります。次の移動先iは複数の可能性があるのでforループで各パターンを確かめます。min()関数を使ってdist[][] + dp[][]の合計距離で最短の地点iを次の移動先にしますが、この時最短のiを経路用配列tourにappendしていきます。次のループでは再度訪問フラグをチェックして、複数ある未訪問の地点iについてまたdist[][] + dp[][]の合計距離を調べます。あとはそれの繰り返しとなります。尚、現在地posはその都度追加していった経路用配列tourの最後の値tour[-1]を参照すればいいことになります。初期状態では0から出発することになるので、tour=[0]にしています。

次へ：TSP LP (Linear Programming)：線形計画法による巡回セールスマン問題 / Subtour Elimination 部分巡回路除去

関連：Traveling Salesman Problem:巡回セールスマン問題について（まとめ）

TSP DP: 巡回セールスマン問題 / 動的計画法 / メモ化再帰

前回の再帰関数と動的計画法の続きとして、今回は巡回セールスマン問題（TSP/Travelling Salesman Problem）の厳密解を求める動的計画法（DP/Dynamic Programming）を試してみました。

再帰関数の例題によく出てくるフィボナッチ数などは単純なので分かりやすいのですが、TSPになると途端に難しくなり、中身を理解するのに予想以上に時間がかかってしまいました。

今回もまた理解した内容をできるだけ自前でプログラミングしてみたので、正式な動的計画法になっているかはわかりませんが一応機能します（DPといってもいろいろやり方があるようなので、どれがいいのかよく分からない）。

調べてみると、ビッド演算子を使うbit DPという方法が便利そうですが、今回はビッド演算子を使わず通常のforループとset()関数で試してみました。

驚きの数学 巡回セールスマン問題

posted with amazlet at 19.05.19

ウィリアム・J・クック
青土社
売り上げランキング: 100,713

Amazon.co.jpで詳細を見る

下準備：（Python 3.6、Jupyter Notebook使用）

まずノードの集合（簡単に4ノード）として、

S = {0, 1, 2, 3}

をPythonのセット関数（集合関数）で用意します。

0からスタートして残りの{1, 2, 3}を巡って最後にまた0に戻ってくることにします。同時に各ノード間の距離テーブルdist[a][b]も用意しておきます。

dist = [[0, 1, 5, 4]
　　　　[1, 0, 3, 6]
　　　　[5, 3, 0, 2]
　　　　[4, 6, 2, 0]]

そして距離を求める再帰関数として、

TSP(a, S-{a}, b)

を用意します。aからスタートして、ノード集合S-{a}を全て巡ってbに戻るという設定です。

実際の数値を代入すれば、

TSP(0, {1,2,3}, 0)

になります。

S-{a}というのは集合関数で、

S-{0} = {0, 1, 2, 3} - {0}　= {1, 2, 3}

となり、集合Sから要素{0}を除去するという意味です。ちなみに要素を付け足すにはプラス（+）ではなくS.add(0)になります。

TSP DPの考え方：

DPの特長として、全体を部分に切り分けて答えを求めていく「分割統治法」というのがあります。その際、部分的な答えを記憶保存しておく「メモ化」というテクニックを使います。
その際求めたい答え（TSPなら全体的な道順）を後回しにして、とりあえず目先の部分的な処理（次の行き先）だけを考えますが、この再帰的な仕組みが理解しにくい。必ずしも再帰ループを使わなくてもいいらしいですが、部分を解いて結果的に全体を解くという流れにすると再帰的になってしまいます。再帰的な途中処理の結果をメモ化で記録して再利用することで効率があがるようです。
まずは、最短経路の合計距離だけを求める式を考えてみます（道順についてはあとで求めます）。

先ほど用意した

TSP(0, {1,2,3}, 0)

というのは、左の0から出発し、次に順番はどうであれ中継地点{1, 2, 3}をそれぞれ一回ずつ通って目的地0に到着するということを表現しています。

TSP(出発地, {中継地点}, 目的地)

というわけです。最終的にはこの式から距離が求まりますが再帰的な処理をしていくので、いきなりは答えがでてきません。
{1, 2, 3}の道順の組み合わせは6通りありますが、とりあえず次に向かう先だけを考えます。そうすると1か2か3の3通りだけを考えることになります。残り2つは気になるけど後回しにしておきます。

次の手続きとして、もし中継地点のうち1に移動すると0から1までの距離dist[0][1]と、次の出発点1として残り{2, 3}が中継地点となります。

それを式に表すと、

TSP(0, {1, 2, 3}, 0) = dist[0][1] + TSP(1, {2, 3}, 0)

ということになります。全体の道順をいきなり解くというより、次の行き先だけを解いて残りは保留（後で解く）という感じです。次の行き先までの距離dist[0][1]は距離テーブルから求められますが、TSP(1, {2, 3}, 0)の距離はこのままでは計算できません。ちなみに、この段階では出発点0から最も近い中継地点に向かうことがベストとは限りません。出発点0から一番遠い地点であっても残りの経路TSP(1, {2, 3}, 0)で挽回することがありうるので、とりあえず式だけ組んでおくという感じです。
そして次の可能性を考えますが、もし1の次に2に向かうなら、

TSP(1, {2, 3}, 0) = dist[1][2]　+ TSP(2, {3}, 0)

と表すことができます。これは1から2へ向かって（距離は計算可能で確定）、2から3へ向かうということになります（予定）。このように確定と予定の二つで式をつくっていきます。あとは中継地点が{3}しかないので、3を次の出発点にして0へ向かうということになりますが、同様の手順を踏んでTSP(2, {3}, 0)の部分を、

TSP(2, {3}, 0) = dist[2][3] + TSP(3, {}, 0)

と表すことができ、さらに

TSP(3, {}, 0) = dist[3][0]

となり、中継地点はないので3から0へ直接移動するということになります。
この一連の流れでは、{}内の要素が一つずつ減っていき空集合になるまで進めていくと、最終的にはdist[a][b]として距離を求められます。また、中継地点のうちどれか一つが次のステップでの出発点となり、左側の変数aは変化し続けます。右側の変数bは最終目的地（出発地でもある）なのでずっと固定値です。
こうやって未定の部分を可能性として掘り下げていくと距離テーブルをつかって実際の値を求めていくことができます。{}内の要素がゼロになったとき、この再帰ループが終わるので、それを終了条件としておけばいいことになります。
あとで気づきましたが、TSP(2, {3}, 0)を求める段階でこの距離はdist[2][3] + dist[3][0]とわかるので、{}内の要素が一つのとき再帰ループの終了条件とし、それ以降のTSP(3, {}, 0)までループさせる必要がないかもしれません。

しかし、上記で行ったのは一つのパターンであって、実際はもっと多くの組み合わせがあります。最初のTSP(0, {1, 2, 3}, 0)のときは、0から次に向かう先は1か2か3の3通りあり、そのうちの最短経路を求めたいのでmin()関数を使って書き換えると、

TSP(0, {1, 2, 3}, 0) = min(dist[0][1] + TSP(1, {2, 3}, 0), 
                           dist[0][2] + TSP(2, {1, 3}, 0), 
                           dist[0][3] + TSP(3, {1, 2}, 0))

になります。これらのうちの最小のものを採用することになります。

しかし、このままではTSP()の部分が計算できないので、先ほどやったように以下のパターン全てを掘り下げて計算していきます。

TSP(1, {2, 3}, 0) = min(dist[1][2] + TSP(2, {3}, 0),
                        dist[1][3] + TSP(3, {2}, 0))

TSP(2, {1, 3}, 0) = min(dist[2][1] + TSP(1, {3}, 0),
                        dist[2][3] + TSP(3, {1}, 0))

TSP(3, {1, 2}, 0) = min(dist[3][1] + TSP(1, {2}, 0),
                        dist[3][2] + TSP(2, {1}, 0))

同様に右辺式後半の各TSP()部分が未確定なので、

TSP(2, {3}, 0) = dist[2][3] + TSP(3, {}, 0)
TSP(3, {2}, 0) = dist[3][2] + TSP(2, {}, 0)
TSP(1, {3}, 0) = dist[1][3] + TSP(3, {}, 0)
TSP(3, {1}, 0) = dist[3][1] + TSP(1, {}, 0)
TSP(1, {2}, 0) = dist[1][2] + TSP(2, {}, 0)
TSP(2, {1}, 0) = dist[2][1] + TSP(1, {}, 0)

となり、さらに

TSP(3, {}, 0) = dist[3][0]
TSP(2, {}, 0) = dist[2][0]
TSP(1, {}, 0) = dist[1][0]

となり、{}が空集合になるまで繰り返すと最後はdist[a][b]だけで値を求められます。あとはここから元に戻っていけば最初の答えがわかるというわけです。このときメモ化をつかって、各段階でTSP(a, S-{a}, b)を記録しておきます。そうすることで、一度計算したパターンについては、また掘り下げて計算しなくて済むようになります。その分計算量が減ってより高速になるというわけです。今回の例ではノード数が少ないため重複する部分も少ないですが（今回の場合は最後の3パターンが一つ前の段階で重複しているだけ）、ノード数が増えるほど重複計算も増えるので効果がでてきます。

実際は以下のようにfor文や再帰処理を使うので、ここまで書き出す必要はなくなります。

TSP（再帰）：
まずはメモ化なしで再帰的な処理だけで組んでみます。

def TSP(a, S, b):
    if len(S) == 0:
        return dist[a][b]

    d_min = float('inf')
    for s in S - {a}:
        d = dist[a][s] + TSP(s, S - {a, s}, b)
        if d < d_min:
            d_min = d

    return d_min

%time print('Total Length:', TSP(0, S, 0))

関数内の最初の処理としては、len(S)で中継地点の数をカウントし、もしゼロ（空集合）であれば、dist[a][b]の値を返すようにしてあります（再帰ループの終了条件）。ここで一旦リターンされますが、内部的にまだ再帰処理が残っている場合は、残りの処理を続けるためにまたループし始めます。
Sが空集合でない場合は、続くforループで複数のパターンを処理させます。
関数TSP(a, S, b)の部分は、最初TSP(0, {0,1,2,3}, 0)であり、関数内でS-{a}={0,1,2,3}-{0}={1,2,3}にしてからforループさせています。次の処理となるforループ内では、S-{a, s}={0,1,2,3}-{0,s}で前回の出発点0と次の目的地sを差し引いて処理し最小値を取り出しています。
もっと完結に書けば、

def TSP(a, S, b):
    if len(S) == 0:
        return dist[a][b]
    return min([dist[a][s] + TSP(s, S - {a, s}, b) for s in S - {a}])

%time TSP(0, S, 0)

このくらいシンプルになりますが、手順が分かりにくいので最初のほうが理解しやすいかと。
再帰処理の場合、終了条件（今回の場合は中継地点Sが空集合になるまで）がなければ無限ループになってしまうので、まずは終了条件を何にするのか、そして内部処理では、終了条件に向かうために何かを変化させていけなければいけないのですが、その部分を何にするかというのが最初は分かりにくいという感じです。終了条件さえ決められれば、while文でループさせてもいいかもしれませんが、自己代入で再帰的に処理するほうが余計な変数など減っていいのかもしれません。

ちなみに、このコードでは最短経路の合計距離（厳密解）だけを求めています。再帰的とはいえ結局はすべてのパターンを計算しています。

このままだと計算量は(n - 1)!なので、10ノードだと1秒くらいかかってしまいます（遅い）。ということから、前回参考にしたメモ化（Memoization）を加えて計算が重複しないように工夫してみました。

TSP DP（メモ化）：

先ほどのコードにメモ化を追加して、一度計算した内容は遡って計算しないように記録しておき、その都度直接呼び出せるようにします。これで一応動的計画法になるようです。
例えば、TSP(1, {2, 3}, 0)の値は、TSP(2, {3}, 0)とTSP(3, {2}, 0)、TSP(3, {}, 0)やTSP(2, {}, 0)、そしてdist[3][0]とdist[2][0]まで遡りつつ比較最小値をとらなければいけませんが、一度最小値の計算結果を求めてあれば、その値をmemoから直接一回で呼び出せます。

memo = {}

def TSP_DP(a, S, b):
    if len(S) == 0:
        memo[(a, tuple(S - {a}), b)] = dist[a][b]
        return dist[a][b]

    d_min = float('inf')
    for s in S - {a}:
        if (s, tuple(S - {a, s}), b) not in memo:
            memo[(s, tuple(S - {a, s}), b)] = TSP_DP(s, S - {a, s}, b)
        d = dist[a][s] + memo[(s, tuple(S - {a, s}), b)]  
        if d < d_min:
            d_min = d

    return d_min

%time print('Total Length:', TSP_DP(0, S, 0))

今回はPythonの辞書機能を使用して、キーをそのままTSP()の引数にしています。辞書dict()では、任意の数や文字を直接キーに割り当てられますが、リスト[]やセット{}は登録できないようです。タプル()なら登録できるので、セット{}をタプル()に変えて登録しています。

結果的には、10ノードで1秒もかかっていた処理速度が34.5 msまで短縮できました。理論上の計算量は2ˆn * nˆ2になるそうです。(n - 1)!よりは高速ですが、ノード数15を超えたあたりからきつくなってくるので、DPとはいえ数十、数百ノードのTSPを処理させることは不可能です。

追記：
その後少しメモ化の部分を変更してみました（以下）。

memo = {}

def TSP_DP(a, S, b):
    S = S - {a}
    if len(S) == 0:
        memo[(a, tuple(S), b)] = dist[a][b]
        return dist[a][b]
    
    if (a, tuple(S), b) in memo:
        return memo[(a, tuple(S), b)]
    
    d_min = min([dist[a][s] +  TSP_DP(s, S - {s}, b) for s in S])
    memo[(a, tuple(S), b)] = d_min

    return d_min

%time print('Total Distance:', TSP_DP(0, S, 0))

辞書機能を使っているのは同じですが、最小値を得たあとに記憶させておくことにしました。この方が効率がいいかもしれません。今回はmin()関数を使い一行forループにし、それとS = S - {a}を最初の方で書いておきもう少しすっきりさせました。

メモ化するためのライブラリfunctools.lru_cacheもありますが、集合関数set()や{}あるいはlist[]を含めてしまうと適用できないようでbit DPにして整数で集合を表現した方がよさそうです（その方がより高速になる）。

経路の出力：

先ほどのコードでは最短距離（厳密解）しか求められないので、経路も求められるコードを追加してみました。経路順はメモ化した内容から割り出せます。最短候補のノードを一つずつ拾い上げてつなげていく感じです。

P = [0]  # starting from to "0"

for i in range(num - 2, -1, -1):
    d_min = np.inf
    p_min = None
    for m in memo:
        if len(m[1]) == i and  set(P) | {m[0]} | set(m[1]) == S:
            d = dist[P[-1]][m[0]] + memo[m]
            if d < d_min:
                d_min = d
                p_min = m[0]
    P.append(p_min)

print('Tour:', P + [0])
plot_path(P + [0])

S={0,1,2,3}の場合のmemoの中身は、

{(1, (), 0): 0.5761741880674179,
 (1, (2,), 0): 0.6102626177520184,
 (1, (2, 3), 0): 1.8033372646495835,
 (1, (3,), 0): 1.5373375235172968,
 (2, (), 0): 0.3706271783270686,
 (2, (1,), 0): 0.8158096274923677,
 (2, (1, 3), 0): 1.7769729629422466,
 (2, (3,), 0): 1.5637018252246337,
 (3, (), 0): 0.7446829688310912,
 (3, (1,), 0): 1.3688287427536237,
 (3, (1, 2), 0): 1.4029171724382241,
 (3, (2,), 0): 1.189646034720611}

となっており、最初の状態(0, {1,2,3}, 0)以降の状態が記録されています。実際for m in memoで呼び出されるのは、(1, (2, 3), 0)というキーの部分だけですが、そのキーを入力すれば対応する距離（計算結果）を一発で呼び出せます。{}がタプルに変換されて()になっています。ちなみに各キーに対応する値（dist[a][b]による距離の計算結果）はnp.random.random()で生成した値を使っているので（以下のGistに記載してあります）このページ最初に設定した値とは異なっています。
最初の出発地点は0、中継地点は{1,2,3}になるので、次の移動先は{1,2,3}のどれかになります。
つまり次のメモ化内容は、出発地が{1,2,3}のどれかで、中継地点の数が2個のものになります。
中継地点{1,2,3}の3通り分だけforループ処理して、さらに検索条件として中継地点の数が2個len(S)==2のものを抜き出します。そうすると、

(1, (2, 3), 0): 1.8033372646495835
(2, (1, 3), 0): 1.7769729629422466
(3, (1, 2), 0): 1.4029171724382241

のうちどれかになります。
TSPの経路として最小値を選択するには、これらの距離だけでは判断できません。これらに一つ前の距離（開始点0からの距離）を足した合計距離で最短のものを選択します。よって式は、

dist[0][1] + memo[(1, (2, 3), 0)]: 1.8033372646495835
dist[0][2] + memo[(2, (1, 3), 0)]: 1.7769729629422466
dist[0][3] + memo[(3, (1, 2), 0)]: 1.4029171724382241

となり、dist[i][j]は距離テーブルから実際の値が求まるので、それにmemo[m]から呼び出した距離を合計して最短経路を選びます。
この結果、最小値のものが、

dist[0][3] + memo[(3, (1, 2), 0)]: 1.4029171724382241

となったとします。
memo[(出発点, (中継地点), 最終到着点)]であるので、出発点3は前回のステップから見れば0からの移動先です。この値3を配列Pにappendしておき最終的な経路を記録していきます。
ということから、さらに次の移動先は中継地点の数が1個で、もうすでに訪れた3が含まれていないということが条件になります。今回の検索条件は

if len(m[1]) == i and  set(P) | {m[0]} | set(m[1]) == S:

としています。mはforループで抜き出したmemo内のキーです。(m[0], m[1], 0)は(出発地, 中継地点, 到着地)に対応しているので、len(m[1])は中継地点の数となります。set(P) | {m[0]} | set(m[1]) == Sの部分は、格納した経路地点と出発地と中継地点を合わせれば常に{0,1,2,3} == Sになるという条件です。 そうすると、次のステップは中継地点数が1個のもので、最終経路の一部としてPに格納されている3が含まれていないもので、出発地が1か2のもの、中継地点も1か2のものとなります。そうすると、

 (1, (2,), 0): 0.6102626177520184
 (2, (1,), 0): 0.8158096274923677

次の移動先はこれら二つに絞られ、さらにこれらにdist[i][j]を加えてから最小値の方を選択します。よって距離の合計式は、

 dist[3][1] + memo[(1, (2,), 0)]: 0.6102626177520184
 dist[3][2] + memo[(2, (1,), 0)]: 0.8158096274923677

となって、

dist[3][1] + memo[(1, (2,), 0)]: 0.6102626177520184

が選ばれたとします。 そうすると1が次の移動先となり、1を配列Pにappendしておきます。いまのところ0→3→1となります。つぎは消去法で2が移動先になりますが、一応先ほどと同じ条件で検索してみます。 次の条件は中継地点数が0個のもので、Pに格納されている[0,3,1]以外のものとなります。そうすると

(2, (), 0): 0.8158096274923677

が抜き出され、Pには2がappendされて[0,3,1,2]になります。 最後に0を加えて[0,3,1,2] + [0] = [0,3,1,2,0]にして最終的な巡回経路にします。

以下はn=10で試したものです。一応厳密解の経路が形成されました。

画像上部のTour:が経路順です。厳密解を見つけることができるのでいいのですが、TSPの規模が大きくなると無理なので、今回はTSPを材料に動的計画法やメモ化について学んでみたという感じです。他のことに応用できるか分かりませんが、それなりの収穫があったと思います。

TSPに関してはこれまでいくつかのヒューリスティックな方法を試しましたが、次は整数計画法を試してみたいです。

プログラミングコンテストチャレンジブック [第2版]　～問題解決のアルゴリズム活用力とコーディングテクニックを鍛える～

posted with amazlet at 19.05.19

秋葉拓哉 岩田陽一 北川宜稔
マイナビ
売り上げランキング: 10,015

Amazon.co.jpで詳細を見る

この↑「プログラミングコンテストチャレンジブック（通称：蟻本）」には様々なアルゴリズムが書いてあります。TSP　DPについても書いてありますが、やはりbit DPを使うと便利らしい。
以下のGitの最後にはmemoizationとbit DPを使ったコードも書いておきました。

続き：TSP DP（その２）　bit DP / 巡回セールスマン問題 / 動的計画法
関連：Traveling Salesman Problem:巡回セールスマン問題について（まとめ）