最近、仕事の都合で「数理最適化」というものに触れる機会に恵まれています。
折角なので、ちょっとだけプライベートで触れてみた結果を載せてみたいと思います。
やったことは、数理最適化の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]
ソースコードをまとめると以下になります。
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で最短らしいので、いくつか答えがあるということになるのですね。
数理最適化は、今まで触れてこなかった分野だけに、とても新鮮です。
気軽に試せる環境があるので、これを機に色々と勉強してみようと思います。