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