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