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

前回はPulp（オープンソース最適化ソルバー/モデラー）を使ってTSP（巡回セールスマン問題）の厳密解をMTZ（Miller-Tucker-Zemlin）部分巡回路除去法で求めてみました。前々回に試したDP（動的計画法）よりはMTZのほうが高速でしたが、せいぜい数十ノード程度が限界でした。

Lazy Constraints（逐次組込制約）：
そして今回もLPで厳密解を求めますが、Lazy Constraintsという制約式を使ってみようと思います。Lazy Constraintsの日本語訳が検索してもあまりでてこないので（ここくらい）、逐次組込制約と呼ぶことにしました。あるいは切除平面法と呼ばれているのかもしれません。
Lazy Constraintsの特長として、最初に基本的な制約だけを与えて一度結果を出し、その結果からその都度必要な制約を追加していき、何回かのループで答えを出します。予め全ての組み合わせの計算をしない分早く収束してくれるときがあります。必要なら計算するという方法が意外に効果的でこれまでの結果を圧倒してくれました。最適化ソルバーPulpのおかげかもしれませんが、これまでのアルゴリズムではせいぜい20ノードが限界でしたが、今回は100ノード以上でも短時間で厳密解を求めることができました。
実際、厳密解かどうかは不安なのでPyconcorde（TSP専用ソルバー）で求めた解と比較することにしました。


これまでの厳密解アルゴリズム比較：（Python 3.6、Jupyter Notebook使用）

全探索法	：20ノード未満で限界
DP（functools.lru_cache）	：20ノード：50.9s
DP（自前のmemoization）	：20ノード：57.7s
DP（ビットDP）	：20ノード：55.8s
LP（MTZ部分巡回路除去）	：20ノード：12.9s
LP（Lazy Constraints）	：20ノード：0.73s、100ノード：23.9s
Pyconcorde（TSPソルバー）	：20ノード：0.23s、100ノード：0.597s

20ノード以上になると極端に時間がかかっていましたが、今回のLazy Constraintsでは100ノードでも数十秒で計算可能でした。もちろんTSP専用ソルバーのPyconcordeなら100ノードでも1秒もかかりません。
ただ、Lazy Constraintsだと状況に応じてその都度制約を追加していくので、早く解決するときもあれば多少時間がかかるときもあり時間にばらつきがあります。今回200ノードまでは数分で厳密解を求めることができましたが、200ノードを越えるとRAMが足りないためなのか途端に処理が遅くなってしまいました。


Lazy Constraintsアルゴリズム：
アルゴリズムにおいては、前回のMTZと比較してもそれほど複雑なことをしているという感じでもありませんでした（ここを参考）。
基本はLPなので、

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

# node i to node j: distance between i and j
dist = dict({(i, j): (abs(XY[i] - XY[j])) for i in range(num) for j in range(num) if i != j})

#  1: the tour includes the edge from i to j, 0: if not included
x = LpVariable.dicts('x', dist, 0, 1, LpBinary)

# the objective to be minimumized
prob += lpSum([dist[(i,j)] * x[(i,j)] for i, j in dist]) 

前回同様、各二点間(i, j)をディクショナリーのキーとして登録しておき、それぞれ値に距離を入れておきます。各ノードの座標は事前にランダム生成してあるので、その座標から距離を求めることができます。x(i, j)は0か1のバイナリで、最終的な巡回路にiからjへの経路が採用されれば1、不採用なら0になります。つまり採用された(i, j)の組だけが距離distとして合算されていきます。
そして最小化する式、

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

を解いていきます。

ベースの制約として、これも前回同様iを基準としてiからjへの移動x(i, j)、あるいはその反対であるjからiへの移動における組み合わせではそれぞれ一通りずつなので、以下の制約式を加えておきます。

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

前回はこのあとMTZの部分巡回路除去の制約を加えていましたが、今回は一旦ここでこの問題を解いてみます。当然中途半端な答えがでてくるので、その答えを見ながら次に必要な制約式を加えていきます。ここがLazyと言われる部分なのでしょう。

この状態で解かせると、大体は以下のように二点間だけの巡回路がたくさんできてしまいます。「0→23→0」など単なる二点間の往復経路ばかり。

100ノードの場合。

最初の二つの制約は、(i, j)マトリクスで考えれば横軸上での組み合わせの中から一つ選び、縦軸上の組み合わせから一つ選ぶという二重制約です。

前回書いた以下のマトリクス（6x6の場合）、

------- 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(1, 2)を選ぶと、同じ横軸上のx(1, 0)、x(1, 3)、x(1, 4)、x(1, 5)は選べなくなり、同様に縦軸上のx(0, 2)、x(3, 2)、x(4, 2)、x(5, 2)も選べなくなるという制約。常に横軸上からひとつ、縦軸上から一つという条件を満たしていないといけません。

以下はシミューレータ（editor.p5js.orgで作成）。グリッド上の(i, j)のどれかを選ぶと次に選ぶことができない組み合わせができます。

グレー：選択可能な(i, j)
黒　　：選択不可能な(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)は選択不可。


選べない経路(i, j)の組み合わせを考える：

この条件だけでは、二点間最短距離のペアばかりできてしまいます。いわゆるグラフ理論における最小マッチングのようなものでしょうか。

これではちょっと先行き不安なので、もう一つだけ簡単な制約式を自分なりに考えて加えておきました。

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

これは、先ほどのマトリクス上で考えると、対角線上に線対称な組み合わせも選べないという制約です。たとえばx(1, 2)を選んだらその横軸と縦軸だけでなくx(2, 1)も選べないということです。これは先ほどの二点間のみの巡回路を制約しており、x(i, j)とx(j, i)は同時に存在しないので、x(i, j)+x(j, i)は1以下になるということです。

この制約を加えてまた一旦解かせてみると、

このように最低でも三点以上の巡回路になります。

少し改善されたように見えますが、実際効率が上がるのかは検証していないので不明。簡単な制約式なのでそれほど計算にも時間かからないのではと思って追加しただけです。かえって時間かかるようであれば、この制約がなくても解くことはできます。
ちなみに、三角形の経路を制約する場合、

lpSum(x[i, j] + x[j, k] + x[k, i]) <= 2

をさらに追加すると最低でも四角形の経路になりますが、3変数もあるのでかえって演算に時間がかかってしまいそうです。今回は2変数だけにして二点間往復路だけを制約してみました。

個人的に興味深いと思ったことは、道順と距離を求めることが目的なのにもかかわらず、それとは別に二点間のマトリクス上でのありえない組み合わせを探しているという部分です。TSPをそのような組み合わせパズルとして見たことはなかったので新鮮な視点が得られた感じがしました。これは一種の抽象的な思考方法なのかもしれませんが、こうやって別の角度から新たなルールを導き出しつつ解答に近くことができるのかもしれません。

ここまでが前準備で、そして以下からがLazy ConstraintsによるSubtour Elimination（部分巡回路除去）です。

Lazy Constraints Subtour Elimination：

制約式としては以下の部分だけです。

   for st in ST:
        nots = [j for j in path if j not in st]
        prob += lpSum(x[(i, j)] for i in st[:-1] for j in nots if (i, j) in x) >= 1
        prob += lpSum(x[(j, i)] for i in st[:-1] for j in nots if (i, j) in x) >= 1

先ほどの100ノードの場合であれば、部分巡回路は19個できてしまいました。この19個の部分巡回路のうち、全体のノード数の半分未満の部分巡回路に対して制約式を与えていきます。全体のノード数の半分以上の大きな部分巡回路を母体として残しておいて、全体のノード数の半分未満の小さな部分巡回路だけに変化（制約）を与えるということです。
今回は、全体のノード数の半分未満の小さな部分巡回路全てに制約を与えますが、代表的なひとつの部分巡回路だけに制約を与えていっても解けるようです。

制約の原理として、例えば「0→23→64→0」という三角形の部分巡回路がありますが、この部分巡回路に必要なのは最低でも一個以上のノードが他の部分巡回路のノードへつながらないといけない（そうしないと現状のまま閉じてしまうので）ということです。三角形内の辺x(0, 23)に対しては、i=0をそのままにしておいてj=23を他の巡回路へつないでもらうために、23が所属している部分巡回路(0, 23, 64)以外のノードをつなげたときに1になるような条件にします。notsという変数は、j=23であれば23が所属している部分巡回路(0, 23, 64)以外の他のノード全部ということです。

つまり、一度解いて出てきた結果x(0, 23)に対して、

x(0はそのまま, 23は{0, 23, 64}以外の全てのノード)

へつなげてみたらどうかという制約になります。

同様に、iとjを反転した場合の制約も加えておきます。

この制約で一発で解決するわけではないのですが、繰り返し解いてはその結果から他のノードへつながるよう（部分巡回路として閉じてしまわないように）仕向けて、最終的に一つの大きな部分巡回路ができあがれば、それが答えとなるわけです。

ちなみにmatplotlibでその経過をアニメーションにしてみました。

Once Loop Reflect

100ノードの場合です。この場合だと6回のループで最適解が求まりました。普通のノートパソコンで23.9sかかりました。

動画が自動ループしているとどこが始まりかわかりにくいですが、最初は細かい三角形や四角形の部分巡回路ばかりで、少しずつ繋がってはまた分解し全体としては徐々に大きい経路へと成長していきます。

一気に解決するような強力なアルゴリズムではなく、必要に応じて徐々に解決していくという部分が面白いと思います。結果的には無駄なことをあまりしていないので早く解決できるのと、その都度計算も重くならないのか100ノードでも現実的な時間内で解決できます。Lazyなのにすごい。

繰り返しループ処理させる部分については、以下のような感じです。

ED = []
while True:
    # Find edges as x(i, j) = 1
    E = [(i, j) for (i, j) in x if x[(i, j)].varValue == 1]
    ED.append(E)
    # Connect the edges to form subtours
    ST = []
    for n in path:
        if n not in sum(ST, []):
            subtour = [n]
            for _ in path:
                for (i, j) in E:
                    if subtour[-1] == i and j not in subtour:
                        subtour.append(j)
                    elif j == subtour[-1] and i not in subtour:
                        subtour.append(i)
            subtour.append(subtour[0])
            ST.append(subtour)

    print('Subtours:', len(ST))
    print(ST)
    if len(ST) == 1:
        print('\nOPT TOUR:', ST[0])
        print('SOLVED')
        break
        
    # Lazy constraints: for the minor subtours(the number of nodes in the subtour is less than half of the whole nodes)
    # Either of (i, j) in the subtour has to connect to a node in the other subtour,
    # hence 'i' stays in the subtour and 'j' has to be changed to one in the other subtour, and vice versa for (j, i)
    for st in ST:
        if len(st) < num / 2: # case of being less than half of the whole node
            nots = [j for j in path if j not in st]
            prob += lpSum(x[(i, j)] for i in st[:-1] for j in nots if (i, j) in x) >= 1
            prob += lpSum(x[(j, i)] for i in st[:-1] for j in nots if (i, j) in x) >= 1
    prob.solve()

T = ST[0]

最初のEDという変数はアニメーション用の配列なので無視して構いません。

whileですべての部分巡回路がひとつながりになるまでループします。

・最初の最適化された結果からx(i, j)=1の辺を探します。

・そして、それらのばらばらの辺(i, j)をつなぎ合わせてすべての部分巡回路（subtour=[0, 23, 64, 0]など）を配列でつくります。

・さらに各部分巡回路を全体用の配列STに入れていきます。len(ST)で何個の部分巡回路があるのかわかります。最終的に1になれば一つの大きな部分巡回路=正解の巡回路というわけです。

・まだ複数の部分巡回路があるうちは、全ノード数の半分未満の小さな部分巡回路stに対してその内容にあわせてLasy Constraintsを与えていきます。各部分巡回路内のx(i, j)に対してiは現状と同じノードでjだけ現在の所属する部分巡回路以外のノードを入れてみるということを繰り返していきます。

・同様に反転したx(j, i)についても制約を与えていきます。

・1ループごとにprob.solve()で暫定的な答えを出してみます。

・全体の部分巡回路の数が一つになるまで上記の繰り返し。

最終的には最適解がでてきます。１ループごとに部分巡回路の状況を出力させています。それをもとにアニメーションをつくりました。

今回は冒頭でも書いたように、本当に最適化されているのかどうかをTSP専用ソルバーであるPyconcordeでも同じ条件で解かせてその結果を比較するということにしました。一応結果は一致したのでアルゴリズム的には機能しているようです。TSP専用ソルバーと比べればまだまだですが、これまでの方法よりも圧倒的に早くかつ解決できるノード数も上回ったので、とりあえずひと段落ついたところです。

次へ：TSP LP: US48 cities / 巡回セールスマン問題　線形計画法

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

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:巡回セールスマン問題について（まとめ）