時間枠付き巡回セールスマン問題に対する定式化とアルゴリズム
from gurobipy import Model, quicksum, GRB
import math
import networkx as nx
時間枠付き巡回セールスマン問題
ここでは，巡回セールスマン問題に時間枠を追加した 時間枠付き巡回セールスマン問題(traveling salesman problem with time windows)を考える．
この問題は，特定の点 $0$ を時刻 $0$ に出発すると仮定し， 点間の移動距離 $c_{ij}$ を移動時間とみなし， さらに点 $i$ に対する出発時刻が最早時刻 $e_i$ と最遅時刻 $\ell_i$ の間でなければならないという制約を課した問題である． ただし，時刻 $e_i$ より早く点 $i$ に到着した場合には，点 $i$ 上で時刻 $e_i$ まで待つことができるものとする．
ポテンシャル定式化
巡回セールスマン問題に対するポテンシャル制約の拡張を考える．
点 $i$ を出発する時刻を表す変数 $t_i$ を導入する． $t_i$ は以下の制約を満たす必要がある． $$ e_i \leq t_i \leq \ell_i \ \ \ \forall i=1,2,\ldots,n $$ ただし， $e_1=0, \ell_1=\infty$ と仮定する．
点 $i$ の次に点 $j$ を訪問する $(x_{ij}=1)$ ときには， 点 $j$ を出発する時刻 $t_j$ は，点 $i$ を出発する時刻に移動時間 $c_{ij}$ を加えた値以上であることから， 以下の式を得る． $$ t_i + c_{ij} - M (1-x_{ij}) \leq t_j \ \ \ \forall i,j, j \neq 1, i \neq j $$ ここで，$M$ は大きな数を表す定数である． なお，移動時間 $c_{ij}$ は正の数と仮定する．$c_{ij}$ が $0$ だと $t_i=t_j$ になる可能性があり， 部分巡回路ができてしまう．これを避けるためには，巡回セールスマン問題と同様の制約を付加する必要があるが， $c_{ij}>0$ の仮定の下では，上の制約によって部分巡回路を除去することができる．
このような大きな数Big Mを含んだ定式化はあまり実用的ではないので，時間枠を用いて強化したものを示す．
$$ \begin{array}{lll} minimize & \sum_{i \neq j} c_{ij} x_{ij} & \\ s.t. & \sum_{j: j \neq i} x_{ij} = 1 & \forall i=1,2,\ldots,n \\ & \sum_{j: j \neq i} x_{ji} = 1 & \forall i=1,2,\ldots,n \\ & t_i + c_{ij} - [\ell_i +c_{ij}-e_j]^+ (1-x_{ij}) \leq t_j & \forall i,j,j \neq 1, i \neq j \\ & x_{ij} \in \{0,1\} & \forall i,j, i \neq j \\ & e_i \leq t_{i} \leq \ell_i & \forall i=1,2,\ldots,n \end{array} $$巡回セールスマン問題のときと同様に，ポテンシャル制約と上下限制約は， 持ち上げ操作によってさらに以下のように強化できる．
$$ t_i + c_{ij} - [\ell_i +c_{ij}-e_j]^+ (1-x_{ij}) + [\ell_i-e_j+\min \{-c_{ji}, e_j-e_i\} ]^+ x_{ji} \leq t_j \ \ \ \forall i,j, j \neq1, i \neq j $$$$ e_i + \sum_{j \neq i} [e_j+ c_{ji} -e_i]^{+} x_{ji} \leq t_i \ \ \ \forall i=2,\cdots,n $$$$ t_i \leq \ell_i -\sum_{j \neq 1,i} [\ell_i -\ell_j+c_{ij}]^+ x_{ij} \ \ \ \forall i=2,\cdots,n $$以下に，強化していない標準定式化 mtztw と強化した定式化 mtz2tw のコードを示す．
def mtztw(n, c, e, l):
"""mtzts: model for the traveling salesman problem with time windows
(based on Miller-Tucker-Zemlin's one-index potential formulation)
Parameters:
- n: number of nodes
- c[i,j]: cost for traversing arc (i,j)
- e[i]: earliest date for visiting node i
- l[i]: latest date for visiting node i
Returns a model, ready to be solved.
"""
model = Model("tsptw - mtz")
x, u, v = {}, {}, {}
for i in range(1, n + 1):
u[i] = model.addVar(lb=e[i], ub=l[i], vtype="C", name="u(%s)" % i)
v[i] = model.addVar(lb=0, ub=n - 1, vtype="C", name="v(%s)" % i)
for j in range(1, n + 1):
if i != j:
x[i, j] = model.addVar(vtype="B", name="x(%s,%s)" % (i, j))
model.update()
for i in range(1, n + 1):
model.addConstr(
quicksum(x[i, j] for j in range(1, n + 1) if j != i) == 1, "Out(%s)" % i
)
model.addConstr(
quicksum(x[j, i] for j in range(1, n + 1) if j != i) == 1, "In(%s)" % i
)
for i in range(1, n + 1):
for j in range(2, n + 1):
if i != j:
M = max(l[i] + c[i, j] - e[j], 0)
model.addConstr(
u[i] - u[j] + M * x[i, j] <= M - c[i, j], "MTZ(%s,%s)" % (i, j)
)
model.addConstr(
v[i] - v[j] + (n - 1) * x[i, j] <= n - 2, "MTZorig(%s,%s)" % (i, j)
)
model.setObjective(quicksum(c[i, j] * x[i, j] for (i, j) in x), GRB.MINIMIZE)
model.update()
model.__data = x, u
return model
def mtz2tw(n, c, e, l):
"""mtz: model for the traveling salesman problem with time windows
(based on Miller-Tucker-Zemlin's one-index potential formulation, stronger constraints)
Parameters:
- n: number of nodes
- c[i,j]: cost for traversing arc (i,j)
- e[i]: earliest date for visiting node i
- l[i]: latest date for visiting node i
Returns a model, ready to be solved.
"""
model = Model("tsptw - mtz-strong")
x, u = {}, {}
for i in range(1, n + 1):
u[i] = model.addVar(lb=e[i], ub=l[i], vtype="C", name="u(%s)" % i)
for j in range(1, n + 1):
if i != j:
x[i, j] = model.addVar(vtype="B", name="x(%s,%s)" % (i, j))
model.update()
for i in range(1, n + 1):
model.addConstr(
quicksum(x[i, j] for j in range(1, n + 1) if j != i) == 1, "Out(%s)" % i
)
model.addConstr(
quicksum(x[j, i] for j in range(1, n + 1) if j != i) == 1, "In(%s)" % i
)
for j in range(2, n + 1):
if i != j:
M1 = max(l[i] + c[i, j] - e[j], 0)
M2 = max(l[i] + min(-c[j, i], e[j] - e[i]) - e[j], 0)
model.addConstr(
u[i] + c[i, j] - M1 * (1 - x[i, j]) + M2 * x[j, i] <= u[j],
"LiftedMTZ(%s,%s)" % (i, j),
)
for i in range(2, n + 1):
model.addConstr(
e[i]
+ quicksum(
max(e[j] + c[j, i] - e[i], 0) * x[j, i]
for j in range(1, n + 1)
if i != j
)
<= u[i],
"LiftedLB(%s)" % i,
)
model.addConstr(
u[i]
<= l[i]
- quicksum(
max(l[i] - l[j] + c[i, j], 0) * x[i, j]
for j in range(2, n + 1)
if i != j
),
"LiftedUB(%s)" % i,
)
model.setObjective(quicksum(c[i, j] * x[i, j] for (i, j) in x), GRB.MINIMIZE)
model.update()
model.__data = x, u
return model
ベンチマーク問題例の読み込み
時間枠付き巡回セールスマン問題のベンチマーク問題例は，以下のサイトから入手できる．
https://lopez-ibanez.eu/tsptw-instances
読み込んでデータを準備する．
folder = "../data/tsptw/"
fn = "n60w20.001.txt"
# fn ="n200w20.002.txt"
# fn ="n100w40.002.txt"
# fn ="n100w60.002.txt"
f = open(folder + fn)
data = f.readlines()
f.close()
x_, y_, early, late = [], [], [], []
for row in data:
try:
L = list(map(float, row.split()))
if L[0] >= 900:
continue
x_.append(L[1])
y_.append(L[2])
early.append(L[4])
late.append(L[5])
except:
pass
n = len(x_)
e, l = {}, {}
for i in range(n):
e[i + 1] = int(early[i])
l[i + 1] = int(late[i])
c = {}
for i in range(n):
for j in range(n):
c[i + 1, j + 1] = int(math.sqrt((x_[i] - x_[j]) ** 2 + (y_[i] - y_[j]) ** 2))
for i in range(1, n + 1):
for j in range(1, n + 1):
for k in range(1, n + 1):
c[i, j] = min(c[i, j], c[i, k] + c[k, j])
model = mtztw(n, c, e, l)
model.optimize()
x, u = model.__data
sol = [i for (v, i) in sorted([(u[i].X, i) for i in u])]
print("Opt.value =", model.ObjVal, sol)
次に，時間枠を用いて制約を強化した定式化を用いて求解する．
model = mtz2tw(n, c, e, l)
model.optimize()
x, u = model.__data
sol = [i for (v, i) in sorted([(u[i].X, i) for i in u])]
print("Opt.value =", model.ObjVal, sol)
G = nx.Graph()
for idx, i in enumerate(sol[:-1]):
G.add_edge(i, sol[idx + 1])
G.add_edge(sol[-1],sol[1])
pos = {i: (x_[i-1], y_[i-1]) for i in range(1,n+1)}
nx.draw(G, pos=pos, node_size=1000 / n + 10, with_labels=False, node_color="blue");
