容量制約付き配送計画問題(capacitated vehicle routing problem, CVRP)に関して、以下のURLを参考にして定式化およびpythonでのシミュレーションを試みました。
参考文献:配送計画問題をpythonで最適化する
URL:https://www.letsopt.com/entry/2020/08/30/180258
作成日:2020/8/30
作成者:cresselia2012
しかし、試みる過程で、以下のMTZ制約に関する制約式とpythonのプログラムがなぜ一致しているのか理解できませんでした。
参考にさせていただいたURL内のpythonのプログラム、特にMTZ制約が正しい理由をご教授いただけますと大変ありがたいです。
理解できなかった理由としては2点あります。
- 参考文献に記載されている本来のMTZ制約の制約式ではu[j]とdemand[j]に入力される顧客jは同様のはずですが、参考文献に記載されているpythonのプログラムではu[j-1]とdemand[j]と記述されており、一見異なる顧客を入力しているように思えます。
### MTZ制約
for i in range(1,num_nodes):
for j in range(1,num_nodes):
if i != j and demand[i] + demand[j] <= capacity:
problem += u[i-1] - u[j-1] + capacity * x[i][j] <= capacity - demand[j]
- 参考文献に記載されている変数u[i]の入力iは0を含まないと定義されています。
しかし、MTZ制約に関するpythonのプログラムではi=1のとき、u[0]が生み出されてしまい、矛盾が起きてしまうのかなと疑問に思いました。
## 変数を定義
u = [ pulp.LpVariable( 'u_{}'.format( i ), demand[i], capacity, cat="Integer" ) \
for i in range(1,num_nodes) ]
### MTZ制約
for i in range(1,num_nodes):
for j in range(1,num_nodes):
if i != j and demand[i] + demand[j] <= capacity:
problem += u[i-1] - u[j-1] + capacity * x[i][j] <= capacity - demand[j]
ちなみに、以下のように本来のMTZ制約の制約式のままでpythonのプログラムを記述して、シミュレーションをしたところ、
for i in range(1,num_nodes):
for j in range(1,num_nodes):
if i != j and demand[i] + demand[j] <= capacity:
problem += u[i] - u[j] + capacity * x[i][j] <= capacity - demand[j]
結果としては以下のようなエラーが起き、プログラムが正しくないことがわかりました。
u[i-1] - u[j-1]をu[i] - u[j]に変更したことでエラーが起こったのは明らかなのですが、具体的な原因はわかりませんでした。
Traceback (most recent call last):
File "C:\Users\Control-PC\Documents\easy_vrp5_compver.py", line 62, in <module>
problem += u[i] - u[j] + capacity * x[i][j] <= capacity - demand[j]
IndexError: list index out of range
最後に私がURLを参考にして記述したpythonのプログラムを記載させていただきます。
以下のプログラムのMTZ制約は参考文献と同様のプログラムを使用しており、シミュレーションを行うとエラーは起きません。
import math
def makeCVRP():
num_nodes = 7
capacity = 30
demand = [ 0, 9, 11, 13, 7,\
19, 17 ]
coordinate = []
coordinate.append( ( 190, 3 ) )
coordinate.append( ( 98, 184 ) )
coordinate.append( ( 5, 42 ) )
coordinate.append( ( 117, 89 ) )
coordinate.append( ( 61, 162 ) )
coordinate.append( ( 9, 97 ) )
coordinate.append( ( 80, 15 ) )
def computeDistance( c1, c2 ):
return math.sqrt( pow( c2[0] - c1[0], 2 ) + pow( c2[1] - c1[1], 2 ) )
distance = [ [ round(computeDistance( c1, c2 )) for c1 in coordinate ] \
for c2 in coordinate ]
return num_nodes, capacity, demand, distance, coordinate
# make problem
num_nodes, capacity, demand, distance, coordinate = makeCVRP()
import pulp
# 最適化問題を定義
problem = pulp.LpProblem( "CVRP", pulp.LpMinimize )
## 変数を定義
x = [ [ pulp.LpVariable( 'x_{}_{}'.format( i, j ), cat="Binary" ) \
if i != j else None for j in range(num_nodes) ] \
for i in range(num_nodes) ]
u = [ pulp.LpVariable( 'u_{}'.format( i ), demand[i], capacity, cat="Integer" ) \
for i in range(1,num_nodes) ]
num_vehicle= [ pulp.LpVariable( 'num_vehicle', 1, num_nodes, cat="Integer" ) ]
## 最小化したい関数を定義
problem += pulp.lpSum( distance[i][j] * x[i][j] for i in range(num_nodes) \
for j in range(num_nodes) if i != j )
## 制約条件を定義
### \sum x_ij = 1
for j in range(1,num_nodes):
problem += pulp.lpSum( x[i][j] for i in range(num_nodes) if i != j ) == 1
for i in range(1,num_nodes):
problem += pulp.lpSum( x[i][j] for j in range(num_nodes) if i != j ) == 1
### \sum x_0j = |K|
problem += pulp.lpSum(x[0][j] for j in range(1,num_nodes)) == num_vehicle
### MTZ制約
for i in range(1,num_nodes):
for j in range(1,num_nodes):
if i != j:
problem += u[i-1] - u[j-1] + capacity * x[i][j] <= capacity - demand[j]
# solve
result = problem.solve(pulp.CPLEX_CMD())
print("objective value = {}".format(pulp.value(problem.objective)))
# pulp(CBC)の結果から辺を定義
edges = [ ( i, j ) for i in range(num_nodes) for j in range(num_nodes)
if i != j and pulp.value(x[i][j]) == 1 ]
# edgesから経路を取得
paths = []
for i,j in edges:
if i == 0:
path = [ i, j ]
while path[-1] != 0:
for v, u in edges:
if v == path[-1]:
path.append(u)
break
paths.append(path)
import networkx as nx
import matplotlib.pyplot as plt
# 有向グラフの作成
G = nx.DiGraph()
G.add_edges_from( edges )
color = [ "r", "g", "y", "m", "c" ]
edge_color = []
# 経路毎に色をセット
for i,j in G.edges:
for t,path in enumerate(paths):
if i in path and j in path:
edge_color.append( color[t] )
break
assert len(edges) == len(edge_color)
# グラフの描画
pos = { i : coordinate[i] for i in range(num_nodes) }
fig = plt.figure()
nx.draw_networkx( G, pos, edge_color=edge_color, alpha=0.5)
# 画像を保存
plt.axis("off")
fig.savefig("test.png")
print(plt.show())
初めての質問であり、質問の仕方が間違っていたら申し訳ありません。
何卒よろしくお願い申し上げます。