TSP LP with scipy.optimize.linprog

今回はscipyを使ってTSP（巡回セールスマン問題）を解いてみました。前回までのPulpのような最適化ソルバー/モデラーであれば0-1整数計画法（混合整数計画法）として計算可能でしたが、scipyの線形計画法の関数であるlinprogは整数計画法まではカバーされていないため、離散値（整数値）ではなく連続値として計算しなければいけません。Pulpを使ったほうがコーディングも楽だし処理も早いと思いますが、今回は敢えてscipyで可能かどうかという試みです。アルゴリズム的には基本的に前回と同じですが、今回は連続値を扱うのとscipyのフォーマットで書き直すため微妙に異なる点があります。データは前回同様TSPLIBのアメリカ48都市を用いることにします。

目的関数とベースの制約：

最小化する式は以下のようになります。

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

s.t: sum(x(0, 1) + x(0, 2) + ... + x(0, 47)) = 2
     sum(x(0, 1) + x(1, 2) + ... + x(1, 47)) = 2
     ・・・
     sum(x(0, 47) + x(1, 47) + ... + x(46, 47)) = 2

とりあえず各都市において最低２本の経路が接続されていなければならないという条件を加えておきます。dist(i, j)は二都市間(i, j)における距離、x(i, j)は経路として採用された場合は1、不採用なら0になりますが、scipyの場合はこの値が連続値になるので0.0〜1.0の小数になります。

上記をscipy.optimize.linprogのフォーマットで書き直すと以下のようになります。

e = [(i, j) for j in range(num) for i in range(num) if i < j]
dist = [abs(XY[i] - XY[j]) for j in range(num) for i in range(num) if i < j]

A_eq = [[1 if n in (i, j) else 0 for (i, j) in e] for n in range(num)]
b_eq = [2] * num
bounds = [(0, 1)] * len(dist)
    
res = linprog(c=dist, A_eq=A_eq, b_eq=b_eq, bounds=bounds, method="revised simplex")

linprogの引数c=distは1次元配列でなければならないので、eとdistは二都市間(i, j)のマトリクスになりますが1次元配列として用意します。48*47/2=1128個の値を含んだ配列になります。

e：(i, j)の二都市を参照するための配列

dist：二都市間の距離を格納した配列

A_eq：各都市における二経路接続の制約マトリクス（二次元配列：1128x48）
b_eq：制約ベクター

bounds：x(i, j)の値の範囲（今回は0.0〜1.0の連続値）

この制約で一度答えを出してみます。methodはより正確なrevised simplexを用いることにします（そうしないとWarningが出る）。まず以下の方法でx(i, j)の値がそれぞれどのくらいになっているか見てみます。

plt.plot(list(range(len(dist))), res.x )

そうすると、以下のように1128通りあるx(i, j)は0か0.5か1.0（縦軸）になっています。

横軸はx(i, j)をx[n]のインデックス番号に変換した値です（n=0〜1127）。実際これを地図に落とし込むと以下のようになります。

赤い線がx(i, j)=0.5、グレーの線がx(i, j)=1.0です。今回はx(i, j)の値が0.0〜1.0の連続値になっているため、この段階では0.5と1.0の二種類の経路が混在しています（もしバイナリ設定が可能なPulpなどのライブラリを使えば0か1の表現になるのでもう少し簡単になるのですが）。

部分巡回路除去：

次の段階として、この状態から部分巡回路を選択して、その部分巡回路に追加の制約を与えます。まず目につくのが、右上にある7都市からなる部分巡回路と下部右寄りにある五角形の部分巡回路です。この段階で選択する部分巡回路の条件は、部分巡回路の都市数が全体の半数未満（48都市あるので24都市未満）でx(i, j)=1.0で構成されているもの（実際は0.5の経路群を含めても可）。複数ある場合は都市数が少ない方をターゲットとします。よって今回は下部右寄りにある五角形{2, 15, 21, 33, 40}を選択します。

ターゲットの部分巡回路の都市群をst_flat={2, 15, 21, 33, 40}として、以下の制約を追加して再度解きます。ST15[0]というのはこの都市群における経路群[(2, 21), (2, 33), (15, 21), (15, 40), (33, 40)]であり、st_flatというのは経路群を都市群に変換したものです。

st_flat = set(sum(ST15[0], ()))
A_ub = [[-1 if len({i, j} & st_flat) == 1 else 0 for (i, j) in e]]
b_ub = [-2]
res = linprog(c=dist, A_ub=A_ub, b_ub=b_ub, A_eq=A_eq, b_eq=b_eq, bounds=bounds, method="revised simplex")

閉じている部分巡回路に対して部分巡回路外の都市と接続するように仕向けるための制約です。x(i, j)のiかjの一方が部分巡回路に属し、もう一方が属していない都市となるx(i, j)の組み合わせが2個以上という制約です。Pulpなどでの制約式としては以下のようになりますが（前回参照）、

sum([x(i, j) for (i, j) in e if len({i, j} & st_flat) == 1]) >= 2

scipy.optimize.linprogでは基本的に「<=」が使用されるため「>=」の場合はをマイナス反転するようです。よってA_ubにおいてx(i, j)が条件に従う場合は-1でそうでない場合は0という配列にしておきます（この部分が多少面倒）。またb_ubのほうにもマイナスをかけて「-2」にしておきます。

この結果、五角形の部分巡回路は以下のようにつながります。

次のターゲットは右上にある7都市の部分巡回路です。同様に制約を追加して解いていきます。

部分巡回路がない場合（接続数１）：

部分巡回路除去の手順を繰り返すと徐々に部分巡回路はつながっていき、以下のような状態になります。

この状態では分離した部分巡回路はありません。このようなひとつながりでありつつ巡回路としては不完全な状態のときには、条件を変えて次のターゲットを選ぶことにします。

今回の場合であれば、左側にある辛うじて一本の経路でつながっている7都市からなる都市群{3, 9, 23, 25, 34, 41, 44}を選択します。この経路はx(i, j)=0.5である赤い経路の都市群{9, 23, 25, 41, 44}とx(i, j)=1.0である{3, 25, 34, 44}と{9, 23}で構成されています。{41}の１箇所で全体に対してつながっているので、この都市群に最低でも二本の経路でつながるような制約を追加します。

似たような箇所として右上にある{5, 16, 18, 26, 27, 29, 35, 36, 42}も{35}の１箇所で辛うじてつながっています。このような１箇所でしか接続されていない都市群を探し出して二経路接続の制約を追加するというのが今回の手段になります。もちろんこれらの都市群の数は全体の都市数の半分未満の数でなければならないという条件と、もし複数ある場合は合計都市数の少ない方を選択するということにしておきます。

条件としては、

・x(i, j)=0.5だけで構成される部分巡回路を探す

・x(i, j)=0.5の部分巡回路に2箇所以上接続しているx(i, j)=1.0の部分巡回路/部分経路を足し合わせる

・足し合わせた都市群は全体都市数の半数未満でなければならない

という感じです。

この部分はかなり面倒な手続きになりますが、そのような都市群を探し出す関数を用意することにしました（あとで登場）。

よって、st_flat={5, 16, 18, 26, 27, 29, 35, 36, 42}に制約を与えると以下のようになります。

左側は部分巡回路として分離しましたがx(i, j)=0.5の赤い経路は消えてくれました。

以後は、これまでと同様に

・分離した部分巡回路（全体の半分未満の都市数からなる）

・分離した部分巡回路がない場合、１箇所でつながっている都市群

という条件で制約を追加していきます。

部分巡回路がない場合（接続数３以上）：

前述の手順を繰り返していると、以下のようなパターンもでてきます。

全体としてはつながっており、１箇所だけでつながっている部分もありません。次の条件としては、x(i, j)=0.5で構成される赤い部分巡回路で他の都市群に接続している部分が３箇所以上あるものを選び、他の都市群への接続数を２箇所以下に制約します。先ほどの方法同様、赤い部分巡回路に対して２箇所以上で接続しているx(i, j)=1.0の部分経路は赤い巡回路に足し合わせます。

具体的には、赤い巡回路{2, 13, 15, 22, 28, 33, 40}に着目します。x(i, j)=1.0である{33, 40}と{2, 15, 21}は、赤い巡回路{2, 13, 15, 22, 28, 33, 40}に対して２箇所で接続しているため、足し合わせて{2, 13, 15, 21, 22, 28, 33, 40}を今回のターゲットとします。

この都市群は{13}と{22}と{28}の3箇所において他の都市群とつながっているので、この接続数３を２以下に制約します。以下がこの部分のコード。

A_ub.append([1 if len({i, j} & st_flat) == 1 else 0 for (i, j) in e])
b_ub.append(2)

この制約で一旦解いてみると、以下のようになります。

分離した部分巡回路ができあがりましたが、赤い経路は消えてx(i, j)=1.0の部分巡回路だけになりました。

次に、最初の方法で右側の部分巡回路に二経路接続の制約を与えると、

部分巡回路がなくなってひとまとまりになります。これで最適解が得られました。

制約を与える手順：

今回は0.5の経路と1.0の経路（実際は0の経路もあるけれども表示されない）があるので、制約を与えるターゲットをどのように選ぶかが複雑でしたが、要約すれば以下のような感じになります。

・ターゲットとなる部分巡回路もしくは足し合わせた都市群の数は全体都市数の半数未満とする

・分離した部分巡回路あるいは都市群がある場合は、その都市群がターゲット

・分離した部分巡回路がない場合は、0.5の部分巡回路とそれに付随する1.0の部分経路を足し合わせた都市群をターゲットとする

・足し合わせた都市群が他の都市群と接続している箇所が、２未満の場合は２以上の制約をあたえる

・足し合わせた都市群が他の都市群と接続している箇所が、３以上の場合は２以下の制約をあたえる

この条件分けで繰り返しループ処理させるコードが以下になります。

def st_elim(ST10, ST05, ST15):
    st_flat = None
    if len(ST15) > 1:
        st_flat = set(sum(ST15[0], ()))
    else:
        if len(ST15[0]) == num:#  and is_loop(ST15[0]):
            st_flat = set(sum(ST15[0], ()))
            print('Optimal')
        else:
            if len(ST05) > 0:
                for st05 in ST05:
                    st05f = set(sum(st05, ()))
                    for st10 in ST10:
                        st10f = set(sum(st10, ()))
                        if len(st05f & st10f) >= 2 and len(st05f | st10f) < num / 2:
                            st05f |= st10f
                            st_flat = st05f

    print('st:',st_flat)
    return st_flat

while True:
    E10 = [e[i] for i, x in enumerate(res.x) if x > 0.99]
    E05 = [e[i] for i, x in enumerate(res.x) if x > 0.49 and x <= 0.51]
    E15 = E10 + E05
    
    ST10 = subtour_list(E10)
    ST05 = subtour_list(E05)
    ST15 = subtour_list(E15)
    
    st_flat = st_elim(ST10, ST05, ST15)
    
    if st_flat == set(range(num)):
        print('Optimal')
        draw_plot(False, E10, E05)
        break
         
    if st_flat is None:
        print('st_flat is None')
        break
    else:            
        cnt = 0
        for i, j in E10:
            if len({i, j} & st_flat) == 1:
                cnt += 1
            
        if cnt > 2:
            A_ub.append([1 if len({i, j} & st_flat) == 1 else 0 for (i, j) in e])
            b_ub.append(2)
            #print('cnt:', cnt)
        else:
            A_ub.append([-1 if len({i, j} & st_flat) == 1 else 0 for (i, j) in e])
            b_ub.append(-2)
        res = linprog(c=dist, A_ub=A_ub, b_ub=b_ub, A_eq=A_eq, b_eq=b_eq, bounds=bounds, method="revised simplex")

最初のst_elim()というのが、ターゲットとなる都市群を探し出す関数です。そして続くwhileループ内でターゲットとなる都市群に条件（他都市との接続数）に応じて制約を追加していきます。今回の場合1ループで約10秒、48都市を解決するのに約2分かかりました。

他のサンプルでは試していないので、まだ完璧ではないかもしれませんが（これから検証）、一応scipyを用いて試してみたという感じです。

関連：TSP: Traveling Salesman Problem / 巡回セールスマン問題について

以下が全体のコード：

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

今回はTSPLIBにある「att48.tsp」というアメリカ48都市のデータを用いて前回の方法（LP:Lazy constraints）を試してみました。前回同様Pulpを使ってLPで解くことにします。TSPLIBには厳密解の巡回路データ「att48.opt.tour」も用意されているので結果を照らし合わせることができます。前回、線形計画法によるTSPの解法がある程度理解できたのでアルゴリズムを整理し直してみました。


事前準備：
まず、ダウンロードした「att48.tsp」をPandasで読み込みます。

data = pd.read_csv('att48.tsp', header=5)
df = data['NODE_COORD_SECTION'].str.split(' ', expand=True)
X = df[1].values[:-1].astype('int32')
Y = df[2].values[:-1].astype('int32')
XY = X + Y * 1j
num = len(X)

データは6行目から各都市の「番号 X座標 Y座標」となっており、主にX座標とY座標だけ抜き出します。最後の行には「EOF（End Of File）」があるので、その部分を含めないようにします。

厳密解データ（解答）も読み込んでおきます。厳密解は巡回順の整数値が並んでいるだけです。

T_opt = pd.read_table('att48.opt.tour', header=4)
T_opt = list(T_opt.values[:-2, 0].astype('int32') - 1) + [0]
draw_plot(T_opt)

表示用draw_plot()関数を用意しておき厳密解を表示させます。

実際このような巡回路となれば正解というわけです。48都市あるので各都市のインデックス番号を0〜47としておきます。TSPLIBでは1〜48の番号になっていますが、とりあえず今回は0〜47で計算することにします。


LPアルゴリズム：
48都市あるので、２都市間経路の組み合わせは48*47/2=1128となります。この場合都市0から都市1への移動となる(0, 1)という組み合わせに対して、その反転した組み合わせ(1, 0)は同じものとみなします（順列ではなく組み合わせのほうが数が減るので）。

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

dist = {(i, j): abs(XY[i] - XY[j]) for i in range(num) for j in range(num) if i < j}
x = LpVariable.dicts('x', dist, 0, 1, LpBinary)

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

for n in range(num):
    prob += lpSum([x[(i, j)] for i, j in dist if n in (i, j)]) == 2
    
result = prob.solve()
print(LpStatus[result])

ED = [(i, j) for (i, j) in x if x[(i, j)].varValue == 1]
draw_plot(False, ED)

前回同様Pulpを使って最小化する式を用意します。ここで注意すべきことは、i<jという大小関係になっていることと、変数xがバイナリLpBinaryになっているということです。変数x[(i, j)]はその組み合わせが経路として採用されれば1で不採用なら0となります。dist[(i, j)]*x[(i, j)]を積算しその値が最小値になるようにLPで解いていきます。
ただし、最低限の制約として、

for n in range(num):
    prob += lpSum([x[(i, j)] for i, j in dist if n in (i, j)]) == 2

があります。この制約の部分を前回よりも簡潔に書き直してみました。巡回路になるためには各都市において必ず二つの経路が接続されていなければならないので、都市1であれば、(i, 1)もしくは(1, j)が二つ必要ということになります。以下のような感じです。
(0, 1), (0, 2), (0, 3)...(0, 46), (0, 47) から二つ採用（都市0において）
(0, 1), (1, 2), (1, 3)...(1, 46), (1, 47) から二つ採用（都市1において）
(0, 2), (1, 2), (2, 3)...(2, 46), (2, 47) から二つ採用（都市2において）
(0, 3), (1, 3), (2, 3)...(3, 46), (3, 47) から二つ採用（都市3において）
...
(0, 46), (1, 46), (2, 46)...(45, 46), (46, 47) から二つ採用（都市46において）
(0, 47), (1, 47), (2, 47)...(45, 47), (46, 47) から二つ採用（都市47において）

一旦この制約でsolve()し、varValue=1となる経路（採用された経路）を表示させます。

そうするとこのような複数の部分巡回路が出来上がります。この段階では各辺がバラバラなので部分巡回路ごとにグループ分けする関数を用意して仕分けします。

def subtour_list(ED):
    E = []
    for (i, j) in ED:
        E_flat = sum(E, [])
        if (i, j) not in E_flat:
            S = [(i, j)]
            while True:
                SC = S.copy()
                for (ii, jj) in ED:
                    if (ii, jj) not in S and (ii, jj) not in E_flat:
                        if {ii, jj} & set(sum(S, ())):
                            S.append((ii, jj))
                if SC == S:
                    E.append(S)
                    break

    EL = [len(e) for e in E]
    AS = np.argsort(EL)
    ST = [E[a] for a in AS]
    return ST

# output subtours
ST = subtour_list(ED)
for st in ST: 
    print(len(st), st)

部分巡回路ごとにグループ分けすると、

4 [(6, 17), (6, 27), (17, 35), (27, 35)]
5 [(2, 21), (2, 33), (15, 21), (15, 40), (33, 40)]
7 [(5, 29), (5, 36), (18, 36), (29, 42), (16, 42), (18, 26), (16, 26)]
13 [(1, 28), (1, 41), (4, 28), (4, 47), (25, 41), (38, 47), (3, 25), (3, 34), (31, 38), (34, 44), (9, 44), (23, 31), (9, 23)]
19 [(0, 7), (0, 8), (7, 37), (8, 39), (14, 39), (30, 37), (30, 43), (43, 45), (11, 14), (32, 45), (10, 11), (10, 22), (13, 22), (13, 24), (19, 32), (19, 46), (20, 46), (12, 20), (12, 24)]

になります。一応、辺数の少ない順に並べ替えています。先頭の値は辺数。


Lasy constraints/部分巡回路除去：
次にこれら部分巡回路のひとつをターゲットとしLazy constraintsを追加し部分巡回路除去を行います。今回の場合は辺数の一番少ない[(6, 17), (6, 27), (17, 35), (27, 35)]に制約を与えることにしました（前回は複数の部分巡回路に制約を与えていましたが）。
部分巡回路は完結した閉路なので、最低でも二本の経路を外につなげるように仕向けます。それが以下の制約になります。

st_flat = set(sum(ST[0], ()))
prob += lpSum([x[(i, j)] for (i, j) in dist if len({i, j} & st_flat) == 1]) >= 2
prob.solve()

(i, j)の組み合わせのうち、iかjのどちらか一方がその部分巡回路に含まれて、他方が含まれていない組み合わせが二つ以上必要という条件にします。「if len({i, j} & st_flat) == 1」の部分が、(i, j)のどちらか一方が部分巡回路に含まれているかどうかという条件です。こうすることで、部分巡回路内にある最低でも二つの都市が他の部分巡回路の都市へつながろうとします。
この条件を与えたあと一旦solve()して結果を見ます。部分巡回路に多少変化が見られるので、何度か繰り返すと部分巡回路がまとまって一つの大きな巡回路になります。while文で最終的にひとつの巡回路になるまで繰り返すコードが以下です。

while len(ST) > 1:
    st_flat = set(sum(ST[0], ()))
    prob += lpSum([x[(i, j)] for (i, j) in dist if len({i, j} & st_flat) == 1]) >= 2
    prob.solve()
    ED = [(i, j) for (i, j) in x if x[(i, j)].varValue == 1]
    ST = subtour_list(ED)
    print('ST[0]:', ST[0])

今回は15ループ（3秒ほど）で48都市を最適化できました。

上記の状態だとまだバラバラの辺の集合なので、以下の処理で最終的な都市順路に並び替えます。

T = list(ED[0])
while len(T) < num:
    for (i, j) in ED:
        if i == T[-1] and j not in T:
            T.append(j)
        elif i not in T and j == T[-1]:
            T.append(i)

T += [T[0]]
print('Optimal Tour:', T)
draw_plot(T)

この結果、

Optimal Tour: [0, 7, 37, 30, 43, 17, 6, 27, 5, 36, 18, 26, 16, 42, 29, 35, 45, 32, 19, 46, 20, 31, 38, 47, 4, 41, 23, 9, 44, 34, 3, 25, 1, 28, 33, 40, 15, 21, 2, 22, 13, 24, 12, 10, 11, 14, 39, 8, 0]

という順路リストが得られます。逆順になっている場合もあるので必要に応じて反転させます。TSPLIBにあるデータは都市番号1から始まっているので、それぞれに+1すれば同じ値になります。


アニメーション：
また変化をアニメーション表示するには以下（前半の設定は省略）。

from matplotlib import animation
from IPython.display import HTML

fig = plt.figure()
ims = []

while True:
    # lazy constraints
    st_flat = set(sum(ST[0], ()))
    prob += lpSum([x[(i, j)] for (i, j) in dist if len({i, j} & st_flat) == 1]) >= 2
    prob.solve()
    
    img = draw_plot(False, ED, ST[0], ANIM=True)
    ims.append(img)
    
    ED = [(i, j) for (i, j) in x if x[(i, j)].varValue == 1]
    ST = subtour_list(ED)
    
    if len(ST) == 1:
        img = draw_plot(False, ED, ST[0], ANIM=True)
        ims.append(img)
        break

ani = animation.ArtistAnimation(fig, ims, interval=800)
#ani.save("tsp_48_1.gif", writer = "imagemagick") # save as animation gif
plt.close()
HTML(ani.to_jshtml())

以下が作成したアニメーション。赤い部分は制約を与える部分巡回路。

Once Loop Reflect

以下が全体のコード。

上記Gist内のコードには(i, j)の組み合わせが倍の2256通りあるパターンについても試しています。組み合わせが多いのでその分時間がかかりますが、同じアルゴリズムで解くことができます。