今回は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：制約ベクター
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を用いて試してみたという感じです。
以下が全体のコード：
