巡回セールスマン問題について

最近、仕事の都合で「数理最適化」というものに触れる機会に恵まれています。
折角なので、ちょっとだけプライベートで触れてみた結果を載せてみたいと思います。

やったことは、数理最適化の1例、「巡回セールスマン問題」
これを、なんと無料の環境pulpで実装できるということで、その環境を整えつつ、実装までやってみました。

やってみたい問題は以下のもの。

◆2.12 巡回セールスマン問題
https://www.msi.co.jp/nuopt/docs/v19/examples/html/02-12-00.html

以下に簡単に絵を描いて見ました。
A,B,C,Dが地図上のポイント。

Aを出発して、一筆書きをして帰ってくるとき、その距離の総和が一番小さいのはどのルートか???
という最適化問題のことを「巡回セールスマン問題」というらしいです。
NP困難と呼ばれて、結構難しいらしい。

しかし、python環境で触れるpulpというモデラーでのアプローチが無料で!!!簡単に出来るということで、以下の記事を参考にさせていただきます。

◆巡回セールスマン問題から始まる数理最適化
https://qiita.com/panchovie/items/6509fb54e3d53f4766aa

素晴らしいことに、pulpのインストールから、サンプルコードまで置いてあります。
早速、最初の設定をサンプルコードに入れてみます。
とりあえずコンパイルが通りそうな以下のようなコードでどうでしょうか。

import pulp as pp

# 最適化モデルの定義
mip_model = pp.LpProblem("tsp_mip", pp.LpMinimize)

N = 4
BigM = 10000

c_ = [
      [0,6,5,5], # pointA > [pointA,pointB,pointC,pointD]
      [6,0,7,4], # pointB > [pointA,pointB,pointC,pointD]
      [5,7,0,3], # pointC > [pointA,pointB,pointC,pointD]
      [5,4,3,0], # pointD > [pointA,pointB,pointC,pointD]
      ]

# 変数の定義
x = [[pp.LpVariable("x(%s,%s)"%(i, j), cat="Binary") for i in range(N)] for j in range(N)]
u = [pp.LpVariable("u(%s)"%(i), cat="Continuous", lowBound=1.0, upBound=(N)) for i in range(N)]

# 評価指標（式（１））の定義＆登録
objective = pp.lpSum(c_[i][j] * x[i][j] for i in range(N) for j in range(N) if i != j)
mip_model += objective

#　条件式(2)の登録
for i in range(N):
    mip_model += pp.lpSum(x[i][j] for j in range(N) if i != j) == 1

# 条件式(3)の登録
for i in range(N):
    mip_model += pp.lpSum(x[j][i] for j in range(N) if i != j) == 1

for i in range(N):
    mip_model += x[i][i] == 0
    
# 条件式(4) (MTZ制約)
for i in range(N):
    for j in range(N):
        if i != j:
            mip_model += u[i] + 1.0 - BigM * (1.0 - x[i][j]) <= u[j]


# 最適化の実行
status = mip_model.solve()

# 結果の把握
print("Status: {}".format(pp.LpStatus[status]))
print("Optimal Value [a.u.]: {}".format(objective.value()))

for i in range(N):
    for j in range(N):
        if i != j:
            print("x[%d][%d]:%f" % (i,j,x[i][j].value()))

for i in range(len(u)):
    print("u[%d] %f" % (i,u[i].value()))

この状態で実行すると･･･

Status: Infeasible
Optimal Value [a.u.]: 18.0
x[0][1]:0.000000
x[0][2]:1.000000
x[0][3]:0.000000
x[1][0]:0.999800
x[1][2]:0.000000
x[1][3]:0.000200
x[2][0]:0.000200
x[2][1]:0.000000
x[2][3]:0.999800
x[3][0]:0.000000
x[3][1]:1.000000
x[3][2]:0.000000

と出てきて、どうやら上手く計算できていないようです。
制約条件が上手く満たせていないような感じなので、色々と検討すると、どうやらこのままではMTZ制約の式がおかしくなるようでした。
最終状態のxからuについての式を4本全て書くと以下になります。

u[0] < u[2] ,\\
u[1] < u[0] ,\\
u[2] < u[3] ,\\
u[3] < u[1]

これを全部満たすとなると、

u[0] < u[2] < u[3] < u[1] < u[0] < ...

なんだか無限に増えていきますが、uの有効範囲は1~4の間のため矛盾が出てしまったようです。
そこで、最後にA地点に戻る変数についての制約を取り除くことで、以下のように式を構成しなおしてみます。

u[0] < u[2] ,\\
u[2] < u[3] ,\\
u[3] < u[1]

これによって、実装も変わり、MTZ条件の式は以下に変更となります。

# 条件式(4) (MTZ制約)
for i in range(N):
    for j in range(1,N):
        if i != j:
            mip_model += u[i] + 1.0 - BigM * (1.0 - x[i][j]) <= u[j]

ソースコードをまとめると以下になります。

sample_route.py

import pulp as pp

# 最適化モデルの定義
mip_model = pp.LpProblem("tsp_mip", pp.LpMinimize)

N = 4
BigM = 10000

c_ = [
      [0,6,5,5], # pointA > [pointA,pointB,pointC,pointD]
      [6,0,7,4], # pointB > [pointA,pointB,pointC,pointD]
      [5,7,0,3], # pointC > [pointA,pointB,pointC,pointD]
      [5,4,3,0], # pointD > [pointA,pointB,pointC,pointD]
      ]

# 変数の定義
x = [[pp.LpVariable("x(%s,%s)"%(i, j), cat="Binary") for i in range(N)] for j in range(N)]
u = [pp.LpVariable("u(%s)"%(i), cat="Continuous", lowBound=1.0, upBound=(N)) for i in range(N)]

# 評価指標（式（１））の定義＆登録
objective = pp.lpSum(c_[i][j] * x[i][j] for i in range(N) for j in range(N) if i != j)
mip_model += objective

#　条件式(2)の登録
for i in range(N):
    mip_model += pp.lpSum(x[i][j] for j in range(N) if i != j) == 1

# 条件式(3)の登録
for i in range(N):
    mip_model += pp.lpSum(x[j][i] for j in range(N) if i != j) == 1

for i in range(N):
    mip_model += x[i][i] == 0
    
# 条件式(4) (MTZ制約)
for i in range(N):
    for j in range(1,N):
        if i != j:
            mip_model += u[i] + 1.0 - BigM * (1.0 - x[i][j]) <= u[j]


# 最適化の実行
status = mip_model.solve()

# 結果の把握
print("Status: {}".format(pp.LpStatus[status]))
print("Optimal Value [a.u.]: {}".format(objective.value()))

for i in range(N):
    for j in range(N):
        if i != j:
            print("x[%d][%d]:%f" % (i,j,x[i][j].value()))

for i in range(len(u)):
    print("u[%d] %f" % (i,u[i].value()))

結果は以下になります。

Status: Optimal
Optimal Value [a.u.]: 18.0
x[0][1]:0.000000
x[0][2]:1.000000
x[0][3]:0.000000
x[1][0]:1.000000
x[1][2]:0.000000
x[1][3]:0.000000
x[2][0]:0.000000
x[2][1]:0.000000
x[2][3]:1.000000
x[3][0]:0.000000
x[3][1]:1.000000
x[3][2]:0.000000
u[0] 1.000000
u[1] 4.000000
u[2] 2.000000
u[3] 3.000000

無事に収束してくれました。
結果はグラフィカルに出ていませんが、

A地点　→　C地点　→　D地点　→　B地点　→　A地点

が最短らしく、その距離は18とのことでした。
問題が載っていたところには、

A地点　→　B地点　→　D地点　→　C地点　→　A地点

で、やはり、距離は18で最短らしいので、いくつか答えがあるということになるのですね。

数理最適化は、今まで触れてこなかった分野だけに、とても新鮮です。
気軽に試せる環境があるので、これを機に色々と勉強してみようと思います。