巡回セールスマン問題とは
ここでは，巡回路型の組合せ最適化問題の代表例である巡回セールスマン問題(traveling salesman problem)を考える．
巡回セールスマン問題は，$n$ 個の点（都市）から構成される無向グラフ $G=(V,E)$， 枝上の距離（重み，費用，移動時間）関数 $c: E \rightarrow \mathbf{R}$ が与えられたとき，すべての点をちょうど $1$ 回ずつ経由する巡回路 で，枝上の距離の合計（巡回路の長さ）を最小にするものを求める問題である．
上の定義のように，向きをもたないグラフ上（これを無向グラフと呼ぶ） で定義された問題を対称巡回セールスマン問題(symmetric traveling salesman problem)と呼ぶ． また，向きをもった（言い換えれば行きと帰りの距離が異なる）グラフ（これを有向グラフ と呼ぶ）上で定義される問題を，非対称巡回セールスマン問題(asymmetric traveling salesman problem)と呼ぶ． もちろん対称巡回セールスマン問題は，非対称巡回セールスマン問題の特殊形なので， (効率良く解けるかどうかは別にして）非対称な問題に対する定式化はそのまま対称な問題に対しても適用できる．
巡回セールスマン問題に対する動的最適化
ある点 $s(\in V)$ から出発し，点の部分集合 $S (\subseteq V)$ をすべて経由し， 点 $j (\in S)$ にいたる最短路長を $f(j,S)$ と書くことにする． このとき，初期条件 $f(j, \{j\})=c_{s j}$ および 次の再帰方程式によって計算された $f(s,V)$ が最短巡回路長である．
$$ f(j, S)= \min_{i \in S \setminus \{ j \}} \left\{ f(i,S \setminus \{ j \}) +c_{ij} \right\} $$点の数を $n$ とする． 始点を $s$ とし， 関数 $f(s,V)$ をよび出すことによって， 動的最適化アルゴリズムは最短巡回路長を出力する． 適当な情報を保存しておくことによって対応する巡回路も再現できる． このアルゴリズムの計算量は $O(n^2 2^n)$であり，必要な記憶容量は $O(n 2^n)$である．
以下に巡回セールスマン問題に対する動的最適化関数 tspdp のコードを示す．
引数：
- n: 点数
- c: 移動費用を表す辞書
- V: 解きたい点集合を表すリスト（Vの要素の対は，必ず辞書cのキーに含まれる必要がある．）
返値：
- 最適値と順回路のタプル
n = 17
random.seed(1)
V, c, x, y = make_data(n)
opt_val, tour = tspdp(n, c, V)
print("Optimal value=", opt_val)
print("Optimal tour=", tour)
G = nx.Graph()
for idx, i in enumerate(tour[:-1]):
G.add_edge(i, tour[idx + 1])
pos = {i: (x[i], y[i]) for i in range(n)}
nx.draw(G, pos=pos, node_size=1000 / n + 10, with_labels=False, node_color="blue");
部分巡回路除去定式化
巡回セールスマン問題を定式化するためには，何通りかの方法がある．まずは，対称巡回セールスマン問題に対する定式化を示す．
枝 $e \in E$ が巡回路に含まれるとき $1$，それ以外のとき $0$ を表す $0$-$1$ 変数 $x_e$ を導入する． 点の部分集合 $S$ に対して，両端点が $S$ に含まれる枝の集合を $E(S)$ と書く． 点の部分集合 $S$ に対して，$\delta(S)$ を端点の1つが $S$ に含まれ，もう1つの端点が $S$ に含まれない枝の集合とする．
巡回路であるためには，各点に接続する枝の本数が2本でなければならない． また，すべての点を通過する閉路以外は，禁止しなければならないので， 巡回路になるためには， 点集合 $V$ の位数 $2$以上の真部分集合 $S \subset V, |S| \geq 2$ に対して，$S$ に両端点が含まれる枝の本数は， 点の数 $|S|$ から $1$ を減じた値以下である必要がある．
上の議論から，以下の定式化を得る． $$ \begin{array}{lll} minimize & \sum_{e \in E} c_e x_e & \\ s.t. & \sum_{e \in \delta(\{i\})} x_e = 2 & \forall i \in V \\ & \sum_{e \in E(S)} x_e \leq |S|-1 & \forall S \subset V, |S| \geq 2 \\ & x_e \in \{0,1\} & \forall e \in E \end{array} $$
点に接続する枝の本数を次数と呼ぶので，最初の制約式は次数制約(degree constraint)と呼ばれる． 2番目の制約は，部分巡回路（すべての点を通らず点の部分集合を巡回する閉路）を除くので， 部分巡回路除去制約(subtour elimination constraint)と呼ばれる．
部分順回路除去制約は，問題の入力サイズの指数オーダーの本数をもつ．よって，必要に応じて現在の解を破っている制約だけを追加する分枝カット法を用いる．
def tsp(V, c):
"""tsp -- model for solving the traveling salesman problem with callbacks
- start with assignment model
- add cuts until there are no sub-cycles
Parameters:
- V: set/list of nodes in the graph
- c[i,j]: cost for traversing edge (i,j)
Returns the optimum objective value and the list of edges used.
"""
EPS = 1.0e-6
def tsp_callback(model, where):
if where != GRB.Callback.MIPSOL:
return
edges = []
for (i, j) in x:
if model.cbGetSolution(x[i, j]) > EPS:
edges.append((i, j))
G = nx.Graph()
G.add_edges_from(edges)
Components = list(nx.connected_components(G))
if len(Components) == 1:
return
for S in Components:
model.cbLazy(quicksum(x[i, j] for i in S for j in S if j > i) <= len(S) - 1)
print ("cut: (%s) <= %s" % (S,len(S)-1) )
return
model = Model("tsp")
# model.Params.OutputFlag = 0 # silent/verbose mode
x = {}
for i in V:
for j in V:
if j > i:
x[i, j] = model.addVar(vtype="B", name="x(%s,%s)" % (i, j))
model.update()
for i in V:
model.addConstr(
quicksum(x[j, i] for j in V if j < i)
+ quicksum(x[i, j] for j in V if j > i)
== 2,
"Degree(%s)" % i,
)
model.setObjective(
quicksum(c[i, j] * x[i, j] for i in V for j in V if j > i), GRB.MINIMIZE
)
model.update()
model.__data = x
return model, tsp_callback
def solve_tsp(V, c):
model, tsp_callback = tsp(V, c)
model.params.DualReductions = 0
model.params.LazyConstraints = 1
model.optimize(tsp_callback)
x = model.__data
EPS = 1.0e-6
edges = []
for (i, j) in x:
if x[i, j].X > EPS:
edges.append((i, j))
return model.ObjVal, edges
obj, edges = solve_tsp(V, c)
print("Optimal tour:", edges)
print("Optimal cost:", obj)
G = nx.Graph()
for (i,j) in edges:
G.add_edge(i, j)
pos = {i: (x[i], y[i]) for i in range(n)}
nx.draw(G, pos=pos, node_size=1000 / n + 10, with_labels=False, node_color="blue");
ポテンシャル定式化
今度は，多項式オーダーの本数の制約をもつ定式化を考えていこう．
非対称巡回セールスマン問題を考える．グラフは有向グラフ $G=(V,A)$ であり， 点集合 $V$，有向枝集合 $A$，枝上の距離関数 $c: A \rightarrow \mathbf{R}$ を与えたとき，最短距離の巡回路を 求めることが目的である．
点 $i$ の次に点 $j$ を訪問するとき $1$，それ以外のとき $0$ になる $0$-$1$ 変数 $x_{ij}$ と点 $i$ の訪問順序を表す実数変数 $u_i$ を導入する． 出発点 $1$ における $u_1$ を $0$ と解釈しておく（実際には $u_1$ は定式化の中に含める必要はない）． 点 $i$ の次に点 $j$ を訪問するときに，$u_j= u_i+1$ になるように制約を付加すると， 点 $1$ 以外の点 $i$ に対しては，$u_i$ は $1,2,\ldots,n-1$ のいずれかの整数の値をとることになる． これらの変数を用いると，非対称巡回セールスマン問題は，以下のように定式化できる．
$$ \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 \\ & u_i + 1 - (n-1) (1-x_{ij}) +(n-3) x_{ji} \leq u_j & \forall i=1,2,\ldots,n, j=2,3,\ldots,n, i \neq j \\ & 1 +(1-x_{1i}) +(n-3) x_{i1} \leq u_{i} & \forall i=2,3,\ldots,n \\ & u_{i} \leq (n-1) - (1-x_{i1}) -(n-3) x_{1i} & \forall i=2,3,\ldots,n \\ & x_{ij} \in \{0,1\} & \forall i \neq j \end{array} $$最初の制約は出次数制約と呼ばれるものであり， 点 $i$ から出る枝がちょうど1 本であることを規定する． 2 番目の制約は入次数制約とよばれ， 点 $i$ に入る枝がちょうど1本であることを規定する． 3番目の制約は，$u_i$ は点 $i$ のポテンシャルと解釈できることから，ポテンシャル制約と呼ばれるものの強化版である． 同様に，上下限制約も強化してある． ポテンシャル制約を用いた定式化は，部分巡回路除去制約を用いた定式化を非対称に拡張したものとくらべると，はるかに弱い． これは，$x_{ij}=1$ になったときのみ，$u_j= u_i+1$ を強制するための制約において， $(1-x_{ij})$ の項の係数 $n-1$ が，非常に大きな数を表す Big M と同じ働きをするためである．
以下に，強化していない定式化 mtz と強化版の定式化 mtz_strong を示す．
n = 17
random.seed(1)
V, c, x, y = make_data(n, start_index=1)
model = mtz(n, c)
model.optimize()
cost = model.ObjVal
print("Opt.value=", cost)
x, u = model.__data
arcs = [(i, j) for (i, j) in x if x[i, j].X > 0.5]
sol = sequence(arcs)
print(sol)
model = mtz_strong(n, c)
model.optimize()
cost = model.ObjVal
print("Opt.value=", cost)
x, u = model.__data
arcs = [(i, j) for (i, j) in x if x[i, j].X > 0.5]
sol = sequence(arcs)
print(sol)
AMPLによる求解
AMPLモデル(mtz.mod)
param n;
set V := {1..n};
param c {i in V, j in V};
var x {i in V, j in V: i != j} binary;
var u {V} >= 0, <= n-1;
minimize z: sum {i in V, j in V: i != j} c[i,j] * x[i,j];
subject to
OutDegree {i in V}:
sum {j in V: j != i} x[i,j] = 1;
InDegree {i in V}:
sum {j in V: j != i} x[j,i] = 1;
Potential {i in 1..n, j in 2..n : i != j}:
u[i] + 1 - (n-1) * (1-x[i,j]) <= u[j];AMPLモデル(mtz_lift.mod)
param n;
set V := {1..n};
param c {i in V, j in V};
var x {i in V, j in V: i != j} binary;
var u {V} >= 0, <= n-1;
minimize z: sum {i in V, j in V: i != j} c[i,j] * x[i,j];
subject to
OutDegree {i in V}:
sum {j in V: j != i} x[i,j] = 1;
InDegree {i in V}:
sum {j in V: j != i} x[j,i] = 1;
Potential {i in 1..n, j in 2..n : i != j}:
u[i] - u[j] + (n-1)*x[i,j] + (n-3)*x[j,i] <= n-2;
LiftedBndOut {i in 2..n}:
-x[1,i] - u[i] + (n-3)*x[i,1] <= -2;
LiftedBndIn {i in 2..n}:
-x[i,1] + u[i] + (n-3)*x[1,i] <= n-2;ampl = AMPL(Environment(AMPL_PATH))
ampl.option['solver'] = 'highs'
# ampl.read("mtz.mod")
ampl.read("../ampl/mtz_lift.mod")
ampl.param['n'] = n
ampl.param['c'] = c
ampl.solve()
z = ampl.obj['z']
x = ampl.var['x']
u = ampl.var['u']
edges = [(i,j) for i in V for j in V if i != j and x[i,j].value() > 0.5]
print("\n\ntsp optimum:", z.value())
print("optimal tour:", edges)
tmp = [(int(u[i].value()+.5), i) for i in V]
tmp.sort()
print("tour =", [i for (u_,i) in tmp])
単一品種フロー定式化
ここでは，「もの」の流れ（フロー）の概念を用いた定式化を紹介する． ここで示すのは単一品種フロー定式化(single commodity flow formulation)と呼ばれる．
いま，特定の点（$1$）に $n-1$単位の「もの」が置いてあり， これを他のすべての点に対してセールスマンによって運んでもらうことを考える． （当然，セールスマンは点 $1$ を出発するものと仮定する．） 点 $1$ からは $n-1$単位の「もの」が出て行き，各点では $1$単位ずつ消費される． また，セールスマンが移動しなた枝にだけ「もの」を流すことができるものとする．
セールスマンが枝 $(i,j)$ を通過することを表す $0$-$1$ 変数 $x_{ij}$の他に，枝 $(i,j)$ を通過する「もの」（品種）の量を表す連続変数を $f_{ij}$ 導入する． これらの記号を用いると，単一品種フロー定式化は以下のように書ける．
$$ \begin{array}{lll} minimize & \sum_{i \neq j} c_{ij} x_{ij} & \\ s.t. & \sum_{j: j \neq i} x_{ij} = 1 & \forall i=0,1,\ldots,n \\ & \sum_{j: j \neq i} x_{ji} = 1 & \forall i=0,1,\ldots,n \\ & \sum_{j} f_{1j} = n-1 & \\ & \sum_{j} f_{ji} -\sum_{j} f_{ij} = 1 & \forall i=2,3,\ldots,n \\ & f_{1j} \leq (n-1) x_{1j} & \forall j \neq 1 \\ & f_{ij} \leq (n-2) x_{ij} & \forall i \neq j, i \neq 1, j \neq 1 \\ & x_{ij} \in \{0,1\} & \forall i \neq j \\ & f_{ij} \geq 0 & \forall i \neq j \end{array} $$ここで最初の2つの制約は次数制約であり，各点に入る枝と出る枝がちょうど1本であることを規定する． 3番目の制約は，最初の点0から $n-1$単位の「もの」が出荷されることを表し， 4番目の制約は「もの」が各点で1ずつ消費されることを表す． 5番目と6番目の制約は容量制約であり，セールスマンが移動しない枝に「もの」が流れないことを表す． ただし，点$0$に接続する枝 $(0,j)$ に対しては最大 $n-1$ 単位の「もの」が流れ， それ以外の枝に対しては最大 $n-2$ 単位の「もの」が流れることを規定している．
def scf(n, c):
"""scf: single-commodity flow formulation for the (asymmetric) traveling salesman problem
Parameters:
- n: number of nodes
- c[i,j]: cost for traversing arc (i,j)
Returns a model, ready to be solved.
"""
model = Model("atsp - scf")
x, f = {}, {}
for i in range(1, n + 1):
for j in range(1, n + 1):
if i != j:
x[i, j] = model.addVar(vtype="B", name="x(%s,%s)" % (i, j))
if i == 1:
f[i, j] = model.addVar(
lb=0, ub=n - 1, vtype="C", name="f(%s,%s)" % (i, j)
)
else:
f[i, j] = model.addVar(
lb=0, ub=n - 2, vtype="C", name="f(%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
)
model.addConstr(quicksum(f[1, j] for j in range(2, n + 1)) == n - 1, "FlowOut")
for i in range(2, n + 1):
model.addConstr(
quicksum(f[j, i] for j in range(1, n + 1) if j != i)
- quicksum(f[i, j] for j in range(1, n + 1) if j != i)
== 1,
"FlowCons(%s)" % i,
)
for j in range(2, n + 1):
model.addConstr(f[1, j] <= (n - 1) * x[1, j], "FlowUB(%s,%s)" % (1, j))
for i in range(2, n + 1):
if i != j:
model.addConstr(f[i, j] <= (n - 2) * x[i, j], "FlowUB(%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, f
return model
model = scf(n, c)
model.optimize()
cost = model.ObjVal
print("Opt.value=", cost)
x, u = model.__data
arcs = [(i, j) for (i, j) in x if x[i, j].X > 0.5]
sol = sequence(arcs)
print(sol)
CONCORDE
巡回セールスマン問題に対する厳密解法と近似解法を実装したものとして CONCORDE があり，以下のサイトからダウンロードできる．
https://www.math.uwaterloo.ca/tsp/concorde/downloads/downloads.htm
アカデミック・非商用の利用は無料であるが，商用の場合には作者のWilliam Cookに連絡をする必要がある． 厳密解法は，巡回セールスマン問題に特化した分枝カット法であり，$85900$点の問題例を解くことに成功している． 近似解法は，反復Lin-Kernighan法 linkern であり，大規模な問題例でも高速に良好な解を算出することができる．
以下に，ベンチマーク問題例をlinkernを用いて求解するための関数を示す．
引数：
- file_name: ファイル名
- folder: ベンチマーク問題例を入れたディレクトリ名
返値：
- total: 近似値
- fig: Plotlyの図オブジェクト
name = "att532"
total, fig = linkern(name + ".tsp")
print("Cost=", total)
plotly.offline.plot(fig);
反復Lin-Kernighan法 linkernを用いて，ベンチマーク問題例以外のデータの近似解を算出する関数 tsplk を作成する． ただし，linkernは8点以下だと計算ができない仕様になっている．そのような場合には，上述の動的最適化を利用すれば良い．
引数：
- n: 点数
- c: 移動費用を表す辞書．キーは必ず点の番号のタプル $i,j ~(i=0,\ldots,n-1, j=0,\ldots,n-1)$ であり，値は非負整数とする
返値：
- total: 近似値
- tour: 順回路のリスト
- G: 順回路を表すnetworkXのグラフオブジェクト
n = 17
random.seed(1)
V, c, x, y = make_data(n)
cc = {}
for (i,j) in c:
cc[i,j] = int(c[i,j]*100000)
total, route, G = tsplk(n, cc)
print("Optimal value=", total)
print("Optimal tour=", route)
pos = {i: (x[i], y[i]) for i in range(n)}
nx.draw(G, pos=pos, node_size=1000 / n + 10, with_labels=False, node_color="blue");
LKH
また，HelsgaunによるLin-Kernighan法(LKH法)の実装 のソースコードが，以下のサイトからダウンロードできる． こちらもアカデミック・非商用のみ無料である．
http://webhotel4.ruc.dk/~keld/research/LKH-3/
ダウンロードしてコンパイルすると，LKH という実行ファイルが生成される（Windowsの場合にはバイナリが提供されている）． LKHは，PyLKHというパッケージ経由で呼び出すことができる．
LKH-3（バージョン3以上）を用いると，巡回セールスマン問題だけでなく，容量制約付き配送計画問題をはじめとした様々な順回路型のベンチーマーク問題を解くことができる． 複数の運搬車は，デポをコピーすることによって巡回セールスマン問題に変換され，容量制約などの付加条件はペナルティとして評価される．
現在，対応している問題は以下の通り．
- ACVRP: Asymmetric capacitated vehicle routing problem（非対称容量制約付き配送計画問題）
- ADCVRP: Asymmetric distance constrained vehicle routing problem（非対称距離制約付き配送計画問題）
- BWTSP: Black and white traveling salesman problem （白黒巡回セールスマン問題；点に白と黒の属性を付加して，連続する黒の点の長さや位数に対する制約が付加された問題）
- CCVRP: Cumulative capacitated vehicle routing problem（累積容量制約付き配送計画問題；各顧客への到着時刻の和を最小化する問題）
- CTSP: Colored traveling salesman problem（色付き巡回セールスマン問題；複数のセールスマンは色をもち，訪問する点は色の部分集合をもつ． 点は指定したいずれかの色をもつセールスマンによって処理されるという条件が付加された問題）
- CVRP: Capacitated vehicle routing problem（容量制約付き配送計画問題）
- CVRPTW: Capacitated vehicle routing problem with time windows（時間枠付き容量制約付き配送計画問題）
- DCVRP: Distance constrained capacitated vehicle routing problem（距離・容量制約付き配送計画問題）
- 1-PDTSP: One-commodity pickup-and-delivery traveling salesman problem（1品種・積み込み・積み降ろし巡回セールスマン問題）
- m-PDTSP: Multi-commodity pickup-and-delivery traveling salesman problem（複数品種・積み込み・積み降ろし巡回セールスマン問題）
- m1-PDTSP: Multi-commodity one-to-one pickup-and-delivery traveling salesman problem （複数品種・1対1・積み込み・積み降ろし巡回セールスマン問題）
- MLP: Minimum latency problem（最小待ち時間問題；顧客の待ち時間の和を最小化する巡回セールスマン問題）
- MTRP: Multiple traveling repairman problem（複数巡回修理人問題；複数の修理人が顧客を訪問する際の，顧客の待ち時間の和を最小化する問題）
- MTRPD: Multiple traveling repairman problem with distance constraints（距離制約付き複数巡回セールスマン問題）
- mTSP: Multiple traveling salesmen problem（複数巡回セールスマン問題）
- OCMTSP: Open close multiple traveling salesman problem（パス型閉路型混在巡回セールスマン問題；あるセールスマンはデポに戻り，他のセーするマンはデポに戻る必要がないパスで良いと仮定した問題）
- OVRP: Open vehicle routing problem（パス型配送計画問題；デポに戻る必要がないと仮定した配送計画問題）
- PDPTW: Pickup-and-delivery problem with time windows（時間枠付き積み込み・積み降ろし型配送計画問題）
- PDTSP: Pickup-and-delivery traveling salesman problem（積み込み・積み降ろし巡回セールスマン問題）
- PDTSPF: Pickup-and-delivery traveling salesman problem with FIFO loading（先入れ先出し型積み込み・積み降ろし巡回セールスマン問題）
- PDTSPL: Pickup-and-delivery traveling salesman problem with LIFO loading（後入れ先出し型積み込み・積み降ろし巡回セールスマン問題）
- RCTVRP: Risk-constrained cash-in-transit vehicle routing problem（リスク制約付き現金輸送配送計画問題）
- RCTVRPTW: Risk-constrained cash-in-transit vehicle routing with time windows（時間枠付きリスク制約付き現金輸送配送計画問題）
- SOP: Sequential ordering problem（先行制約付き巡回セールスマン問題）
- STTSP: Steiner traveling salesman problem（Steinr巡回セールスマン問題）
- TRP: Traveling repairman problem（巡回修理人問題；顧客の待ち時間の和を最小化する問題）
- TSPDL: Traveling salesman problem with draft limits（喫水制限付き巡回セールスマン問題）
- TSPPD: Traveling salesman problem with pickups and deliveries（喫水制限付き積み込み・積み降ろし巡回セールスマン問題）
- TSPTW: Traveling salesman problem with time windows（時間枠付き巡回セールスマン問題）
- VRPB: Vehicle routing problem with backhauls（帰り荷を考慮した配送計画問題）
- VRPBTW: Vehicle routing problem with backhauls and time windows（帰り荷を考慮した時間枠付き配送計画問題）
- VRPMPD: Vehicle routing problem with mixed pickup and delivery（配送と集荷を考慮した配送計画問題）
- VRPMPDTW: Vehicle routing problem with mixed pickup and delivery and time windows（配送と集荷を考慮した時間枠付き配送計画問題）
- VRPSPD: Vehicle routing problem with simultaneous pickup and delivery（同時配送集荷を考慮した配送計画問題）
- VRPSPDTW: Vehicle routing problem with simultaneous pickup-delivery and time windows（同時配送集荷を考慮した時間枠付き配送計画問題）
以下のコードでは，配送計画問題をベンチマークサイトからダウンロードして解いている．
problem_str = requests.get('http://vrp.atd-lab.inf.puc-rio.br/media/com_vrp/instances/E/E-n51-k5.vrp').text
problem = lkh.LKHProblem.parse(problem_str)
solver_path = './LKH'
tour = lkh.solve(solver_path, problem=problem)
folder = "../data/cvrp/"
file_name = "E-n51-k5.vrp"
m = 5 #運搬車の台数；ベンチマーク問題例のk?の?を代入
f = open(folder + file_name)
data = f.readlines()
f.close()
n = int(data[3].split()[-1])
Q = int(data[5].split()[-1])
print("n=",n, "Q=",Q)
x, y= {}, {}
pos = {}
for i, row in enumerate(data[7 : 7 + n]):
id_no, x[i], y[i] = list(map(int, row.split()))
pos[i] = x[i], y[i]
G = nx.Graph()
i = 1
for j in tour[0][1:]:
if j>n:
j=1 #depot
G.add_edge(i-1,j-1)
i = j
G.add_edge(j-1,0)
nx.draw(G, pos=pos, node_size=1000 / n + 10, with_labels=False, node_color="blue")
plt.show()
def tsp_length(d, # Distance matrix
tour): # Order of the cities
n = len(tour)
length = d[tour[n - 1]][tour[0]]
for i in range(n - 1):
length += d[tour[i]][tour[i + 1]]
return length
######### Build solution representation by successors
def tsp_tour_to_succ(tour):
n = len(tour)
succ = [-1] * n
for i in range(n):
succ[tour[i]] = tour[(i+1)%n]
return succ
######### Build solution representation by predeccessors
def tsp_succ_to_pred(succ):
n = len(succ)
pred = [-1] * n
for i in range(n):
pred[succ[i]] = i
return pred
######### Convert solution from successor of each city to city order
def tsp_succ_to_tour(succ):
n = len(succ)
tour = [-1] * n
j = 0
for i in range(n):
tour[i] = j
j = succ[j]
return tour
######### Convert a solution given by 2-opt data structure to a standard tour
def tsp_2opt_data_structure_to_tour(t):
n = int(len(t)/2 + 0.5)
tour = [-1] * n
j = 0
for i in range(n):
tour[i] = j>>1
j = t[j]
return tour
######### Compare 2 directed tours; returns the number of different arcs
def tsp_compare(succ_a, succ_b):
n = len(succ_a)
count = 0
for i in range(n):
if succ_a[i] != succ_b[i]:
count += 1
return count
######### Basic Lin & Kernighan improvement procedure for the TSP
def tsp_LK(D, # Distance matrix
tour, # Solution
length): # Tour length
n = len(tour)
succ = tsp_tour_to_succ(tour)
for i in range(n): succ[tour[i]] = tour[(i + 1) % n]
tabu = [[0 for _ in range(n)] for _ in range(n)] #Can edge i-j be removed ?
iteration = 0 # Outermost loop counter to identify tabu condition
last_a, a = 0, 0 # Initiate ejection chain from city a = 0
improved = True
while a != last_a or improved:
improved = False
iteration += 1
b = succ[a]
path_length = length - D[a][b]
path_modified = True
while path_modified: # Identify best ref. struct. with edge a-b removed
path_modified = False
ref_struct_cost = length # Cost of reference structure retained
c = best_c = succ[b]
while succ[c] != a: # Ejection can be propagated
d = succ[c]
if path_length - D[c][d] + D[c][a] + D[b][d] < length:
best_c = c # An improving solution is identified
ref_struct_cost = path_length - D[c][d] + D[c][a] + D[b][d]
break # Change improving solution immediately
if tabu[c][d] != iteration and \
path_length + D[b][d] < ref_struct_cost:
ref_struct_cost = path_length + D[b][d]
best_c = c
c = d # Next value for c and d
if ref_struct_cost < length: # Admissible reference structure found
path_modified = True
c, d = best_c, succ[best_c] # Update reference structure
tabu[c][d] = tabu[d][c] = iteration#Don't remove again edge c-d
path_length += (D[b][d] - D[c][d])
i, si, succ[b] = b, succ[b], d # Reverse path b -> c
while i != c:
succ[si], i, si = i, si, succ[si]
b = c
if path_length + D[a][b] < length: # A better solution is found
length = path_length + D[a][b]
succ[a], last_a = b, b
improved = True
tour = tsp_succ_to_tour(succ)
succ = tsp_tour_to_succ(tour)
a = succ[a]
return tour, length
賞金収集巡回セールスマン問題とその変形
巡回セールスマン問題の一般化として，賞金収集巡回セールスマン問題 (prize collecting traveling salesman problem) がある． これは，部分順回路問題の一種であり，点上に定義された賞金を収集することから，この名がつけられた．もともとは，（両者とも研究が進んでいる）ナップサック多面体と巡回セールスマン多面体を合わせるとどうなるかという理論的な興味から始まったが， 近年では，活発に研究されている問題である．
賞金収集巡回セールスマン問題は， $n$ 個の点（都市）から構成されるグラフ $G=(V,E)$， 枝上の費用関数 $c: E \rightarrow \mathbf{R}$，点上の賞金関数 $p: V \rightarrow \mathbf{R}$， 賞金収集額の下限 $b$ が与えられたとき， 点の部分集合をちょうど $1$ 回ずつ経由する巡回路で，部分集合内の賞金の合計が $b$ 以上で， 順回路の費用の合計を最小にするものを求める問題である.
賞金収集巡回セールスマン問題が，賞金の下限を収集する部分順回路を求めるのに対して，オリエンテーリング問題 (orienteering problem)は，与えられた巡回費用の上限以下で，収集した賞金の合計を最大化する． これは，その名前の由来になったオリエンテーリングという競技（制限時間内に，地点に付与されたスコアを収集し，その合計を競う）の他に，様々な応用をもつ．
オリエンテーリング問題は，$n$ 個の点（都市）から構成されるグラフ $G=(V,E)$， 枝上の移動時間関数 $c: E \rightarrow \mathbf{R}$，点上の賞金（スコア）関数 $p: V \rightarrow \mathbf{R}$， 移動時間の上限$T$が与えられたとき， 点の部分集合をちょうど $1$ 回ずつ経由する巡回路で，順回路の移動時間の合計が$T$以下で，部分集合内の賞金の合計を最大にするものを求める問題である.
定式化はほぼ同じなので，オリエンテーリング問題のものだけを示す．
部分巡回路除去定式化
まずは，対称オリエンテーリング問題に対する定式化を示す． ただし，2点から成る順回路を除いた解の中に最適解があると仮定して定式化を行う． （2点から成る順回路の数は $O(n^2)$ なので，事前にチェックしておけば良い．）
枝 $e \in E$ が巡回路に含まれるとき $1$，それ以外のとき $0$ を表す $0$-$1$ 変数 $x_e$ と，点を巡回するか否かを表す $0$-$1$ 変数 $y_i$ を導入する． 点の部分集合 $S$ に対して，$\delta(S)$ を端点の1つが $S$ に含まれ，もう1つの端点が $S$ に含まれない枝の集合とする．
$$ \begin{array}{lll} maximize & \sum_{i \in V} p_i y_i & \\ s.t. & \sum_{e \in \delta(\{i\})} x_e = 2 y_i & \forall i \in V \\ & \sum_{e \in \delta(S)} x_e \geq 2(y_i+y_j-1) & \forall S \subset V, i \in S, j \in V \setminus S, |S| \geq 2 \\ & \sum_{e \in E} c_e x_e \leq T & \\ & x_e \in \{0,1\} & \forall e \in E \\ & y_i \in \{0,1\} & \forall i \in V \end{array} $$部分順回路制約は，カットセット型で記述されている．これは，点の部分集合 $S$ 内の点 $i$ と $S$ 外の点 $j$ が部分順回路に含まれているときには， $S$ と $S$ 以外の間に2本以上の枝があることを意味する．
ポテンシャル定式化
非対称巡オリエンテーリング問題を考え，多項式オーダーの本数の制約をもつ定式化を示す．
点 $i$ の次に点 $j$ を訪問するとき $1$，それ以外のとき $0$ になる $0$-$1$ 変数 $x_{ij}$，点を巡回するか否かを表す $0$-$1$ 変数 $y_i$， ならびに 点 $i$ の訪問順序を表す実数変数 $u_i$ を導入する．
出発点 $1$ を出発点と考え， $u_1$ を $0$ と解釈しておく（実際には $u_1$ は定式化の中に含める必要はない）． 点 $i$ の次に点 $j$ を訪問するときに，$u_j= u_i+1$ になるように制約を付加する．
非対称オリエンテーリング問題は，以下のように定式化できる．
$$ \begin{array}{lll} maximize & \sum_{i \in V} p_i y_i & \\ s.t. & \sum_{j: j \neq i} x_{ij} = y_i & \forall i=1,2,\ldots,n \\ & \sum_{j: j \neq i} x_{ji} = y_i & \forall i=1,2,\ldots,n \\ & u_i + 1 - (n-1) (1-x_{ij}) \leq u_j & \forall i=1,2,\ldots,n； j=2,3,\ldots,n, i \neq j \\ & 1 \leq u_{i} \leq (n-1) & \forall i=2,3,\ldots,n \\ & \sum_{i \neq j} c_{ij} x_{ij} \leq T & \\ & x_{ij} \in \{0,1\} & \forall i \neq j \\ & y_i \in \{0,1\} & \forall i \in V \end{array} $$ベンチマーク問題例は，以下のサイトからダウンロードできる．
https://www.mech.kuleuven.be/en/cib/op
オリエンテーリング問題は，制限時間内になるべく価値の高い仕事をこなすという自然なモデルのため， 以下のような様々な分野に応用をもつ．
- モバイル・クラウドソーシング（時間枠）
- 旅程計画（時間枠と時刻依存移動時間）
- テーマパークのアトラクション巡回順（時刻依存移動時間）
- 在庫配送計画
folder = "../data/orienteering/"
fn ="tsiligirides_problem_1_budget_20.txt"
f = open(folder+fn)
data = f.readlines()
f.close()
Tmax, num_paths = tuple(map(int, data[0][:-1].split( "\t" )))
#print(Tmax)
pos, x, y = {},{},{}
prize ={}
for i, row in enumerate(data[1:]):
try:
x[i+1], y[i+1], p = list(map(float, row.split( "\t")))
pos[i+1] = (x[i+1], y[i+1])
prize[i+1] = p
except:
pass
n = len(prize)
c ={}
for i in range(1,n+1):
for j in range(1,n+1):
c[i,j] = distance(x[i], y[i], x[j], y[j])
model = Model("orienteering - mtz")
x, y, u = {}, {}, {}
for i in range(1, n + 1):
u[i] = model.addVar(lb=0, ub=n - 1, vtype="C", name= f"u({i})")
y[i] = model.addVar(vtype="B", name= f"y({i})")
for j in range(1, n + 1):
if i != j:
x[i, j] = model.addVar(vtype="B", name= f"x({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) == y[i], f"Out({i})"
)
model.addConstr(
quicksum(x[j, i] for j in range(1, n + 1) if j != i) == y[i], f"In({i})"
)
for i in range(1, n + 1):
for j in range(2, n + 1):
if i != j:
model.addConstr(
u[i] - u[j] + (n - 1) * x[i, j] + (n - 3) * x[j, i] <= n - 2,
"LiftedMTZ(%s,%s)" % (i, j),
)
for i in range(2, n + 1):
model.addConstr(
-x[i, 1] + u[i] + (n - 3) * x[1, i] <= n - 2, name="LiftedUB(%s)" % i
)
model.addConstr( quicksum(c[i, j] * x[i, j] for (i, j) in x) <= Tmax, "Time Constraint" )
model.setObjective(quicksum(prize[i] * y[i] for (i) in y), GRB.MAXIMIZE)
model.optimize()
cost = model.ObjVal
print("Opt.value=", cost)
arcs = [(i, j) for (i, j) in x if x[i, j].X > 0.5]
G = nx.Graph()
G.add_edges_from(arcs)
G.add_nodes_from(list(range(1,n+1)))
nx.draw(G, pos=pos, node_size=1000 / n + 10, with_labels=False, node_color="blue");
階層的巡回セールスマン問題
2段階の巡回セールスマン問題を考える．この問題は，都市内物流や緊急物資の配送への応用をもつ．
- $P$: 中継地点およびデポの集合
- $N$: 顧客の集合
- $E_1$: 階層1の枝（中継地点およびデポ間）の集合 $E_1=\{(i,j)\ |\ i,j \in P, i < j\}$
- $E_2$: 顧客間の枝の集合 $E_2=\{(i,j) | i,j \in N, i < j \}$
$E_3$: 顧客と中継地点との間の枝の集合 $E_3=\{(i,p)\ |\ i \in N, p \in P\}$
$x_e$: 配送車が枝$e$上を移動する場合1,そうでないとき0をとる $0$-$1$ 変数．
- $X_e$: 台車が枝$e$上を移動する場合1,そうでないとき0をとる$0$-$1$ 変数．
- $y_p$: 配送車が中継地点$p$に停車する場合1,そうでないとき0をとる$0$-$1$ 変数．
- $z_{ip}$: 顧客$i$に対し，中継地点$p$から往復輸送が行われる場合に1,そうでないとき0をとる$0$-$1$ 変数．
点の部分集合$S (\subset N \cup P)$について，以下を定義する．
$$ \delta_k \left( S \right) = \left\{ (i,j)\ |\ (i,j) \in E_k, (i \in S, j \not\in S\ \text{or}\ i \not\in S, j \in S) \right\} \quad (k=1,2,3)$$
また枝の集合$E$について，以下を定義する．
$$ x(E) = \sum_{e \in E} x_e, \quad X(E) = \sum_{e \in E} X_e, \quad z(E) = \sum_{e \in E} z_e$$
定式化 $$ \begin{array}{llll} minimize & { \displaystyle\sum_{e\in E_1} c_e^1x_e + \displaystyle\sum_{e\in E_2\cup E_3} c_e^2X_e + \displaystyle\sum_{e\in E_3} 2c_e^2 z_e }&& \\ s.t. & X\left(\delta_2\left(\{i\}\right)\cup\delta_3\left(\{i\}\right)\right) + 2z\left(\delta_3\left(\{i\}\right)\right) = 2 & \forall i\in N & (1) \\ & X\left(\delta_2\left(\{p\}\right)\cup\delta_3\left(\{p\}\right)\right) + 2z\left(\delta_3\left(\{p\}\right)\right) = 2y_p & \forall p \in P & (2) \\ & x\left(\delta_1\left(\{p\}\right)\right) = 2y_p & \forall p \in P & (3) \\ & X\left(\delta_2 (S) \cup \delta_3 (S)\right) + 2 z\left(\delta_3 ( S )\right) \geq 2 & \forall S\subset N, |S|\geq 3 & (4) \\ & x\left(\delta_1 (S)\right) \geq 2 (y_p+y_q-1) & \forall S\subset P, |S| \geq 3, p\in S, q \in P\setminus S & (5) \\ & \displaystyle\sum_{(k,p) \in \delta_3( \{k\} ), p \in I} X_{kp} + \displaystyle\sum_{(\ell,q) \in \delta_3(\{\ell\}), q \in (P \setminus I)} X_{\ell q}+\displaystyle\sum_{i,j \in S \cup \{k,\ell\}, i < j} X_{ij} \leq |S| +2 & \forall (S\cup\{k,\ell\})\subset N, k\ne \ell, I\subset P & (6)\\ & x_e \in \{0,1\} & \forall e \in E_1 & \\ & X_e \in \{0,1\} & \forall e \in E_2 \cup E_3 & \\ & z_e \in \{0,1\} & \forall e \in E_3 & \\ & y_p \in \{0,1\} & \forall p \in P & \end{array} $$
制約式(1),(2),(3)はいずれも次数制約である．式(1)は顧客$i (\in N)$から伸びる枝は2本であることを表す． 式(2)は中継地点$p (\in P)$を用いる場合，その点から階層2に向かって伸びる枝は2本であることを表す． 式(3)は中継地点$p (\in P)$を用いる場合，その点から階層1に向かって伸びる枝は2本であることを表す．
制約式(4)は，顧客の集合$S$に対するカットセット制約である．
制約式(5)は，中継地点の集合$S$に対するカット制約である． 2点$p (\in S)$，$q (\not\in S)$を訪れる場合，集合$S$と$P \setminus S$の間に2本以上の枝が存在しなければならないことを意味する．
制約式(6)は，階層2において中継地点を出発して，異なる中継地点へ向かうことを禁じるパス除去制約である．
カットセット制約ならびにパス除去制約の数は非常に多いので，実装の際には切除平面法（分枝カット法）を用いる必要がある．
階層的巡回セールスマン問題に対するルート先・クラスター後法
ルート先・クラスター後法（route-first/cluster-second method）では，はじめにすべての点（顧客およびデポ）を通過する 巡回路を（たとえば巡回セールスマン問題を解くことによって）作成し，その後でそれをクラスターに分けることによってルートを生成する．
階層的巡回セールスマン問題に対して，すべての点を通過する巡回路を，その巡回路の順番を崩さないように「最適」に分割する方法を考える．
点（顧客およびデポ）の集合 $N$ をちょうど1回ずつ通過する巡回路を表す順列を $\rho$とする．
$\rho (i)$ は $i$ 番目に通過する点の番号であり， $\rho(0)$ はデポ $(0)$ である． $C_{ij}^p$ を $i+1$ 番目から $j$ 番目の顧客を $\rho$ で定義される順に，中継地点 $p$ から巡回したときのルートの費用と定義する．
ただし，ルートに何らかの制限がないと，階層的な解にならない．1つのルートに含まれる顧客数の上限，もしくはルート長に制約を付加する必要がある． 以下の例題では，階層2の1つのルートに含まれる顧客数が $Q$ 以下であるという制約を付加する．
すべての中継地点を使う場合は，点集合 $P$ に対する最適な順回路を計算し，その後，顧客の中継地点への割り当ては，動的最適化で行う．
$i+1$ 番目から $j$ 番目の顧客を巡回するための最適な中継地点は， 以下のように計算できる．
$$ p^*_{ij} = \arg \min_{p} C_{ij}^p $$$C_{ij}^{p^*_{ij}}$ を枝の費用としたとき，点 $0$ から $n=|N|$ までの（有向閉路をもたないグラフ上での）最短路は，動的最適化で計算できる． 最短路に対応する巡回路 $\rho$ が，最適な分割になる．
$j$番目の点までの最適値を $F_j$ とする． $F_0 = 0$ の初期条件の下で，以下の再帰方程式によって最適値を得ることができる．
$$ F_j = \min \{ F_i +C_{ij}^{p^*_{ij}} \} \ \ j=1,2,\ldots,n $$ランダムな問題例を作成する．中継地点数は $30$，顧客数は $100$ とし，中央にデポがあるものとする．デポの番号は $0$ とする． 顧客とデポに対して，tsplk関数を用いて巡回セールスマン問題の近似解を算出し，それを順列 $\rho$ （プログラムではroute2）とする．
Q = 5 # capacity
n1, n2 = 30, 100
x1, y1, x2, y2 = {}, {}, {}, {}
for i in range(1, n1):
x1[i] = random.random() * 100.0
y1[i] = random.random() * 100.0
x1[0], y1[0] = 50.0, 50.0 # depot
for i in range(1, n2):
x2[i] = random.random() * 100.0
y2[i] = random.random() * 100.0
x2[0], y2[0] = 50.0, 50.0 # depot
def distance(x1, y1, x2, y2):
"""distance: euclidean distance between (x1,y1) and (x2,y2)"""
return int(math.sqrt((x2 - x1) ** 2 + (y2 - y1) ** 2) + 0.5)
c2 = {}
pos2 = {}
for i in range(n2):
c2[i, i] = 0
pos2[i] = x2[i], y2[i]
for j in range(i, n2):
c2[j, i] = c2[i, j] = distance(x2[i], y2[i], x2[j], y2[j])
total2, route2, G2 = tsplk(n2, c2)
nx.draw(G2, pos=pos2, with_labels=True)
print(total2, route2)
動的最適化の準備として， $i+1$ 番目から $j$ 番目の顧客を $\rho$ で定義される順に，最適な中継地点 $p$ から巡回したときのルートの費用 $C_{ij}^{p^*_{ij}}$ を計算しておく．
C = {}
pstar = {}
for i in range(n2 - 1):
cum = 0
for j in range(i + 1, n2 - 1):
if j - i > Q:
break
C[i, j] = cum
cum += c2[route2[j], route2[j + 1]]
# find the nearest parking point
min_ = np.inf
for p in range(n1):
dis = distance(
x2[route2[i + 1]], y2[route2[i + 1]], x1[p], y1[p]
) + distance(x2[route2[j]], y2[route2[j]], x1[p], y1[p])
if dis < min_:
min_ = dis
pstar[i, j] = p
C[i, j] += min_
if n2 - 1 - i <= Q:
C[i, n2 - 1] = cum
min_ = np.inf
for p in range(n1):
dis = distance(
x2[route2[i + 1]], y2[route2[i + 1]], x1[p], y1[p]
) + distance(x2[route2[n2 - 1]], y2[route2[n2 - 1]], x1[p], y1[p])
if dis < min_:
min_ = dis
pstar[i, n2 - 1] = p
C[i, n2 - 1] += min_
動的最適化の再帰方程式によって，最適なルートの分割を求め，使用した中継地点をParkに保管する．
F = {}
prev = {}
F[0] = 0
for j in range(1, n2):
min_ = np.inf
for i in range(j - 1, -1, -1):
if (i, j) not in C:
break
if F[i] + C[i, j] < min_:
min_ = F[i] + C[i, j]
prev[j] = i
F[j] = min_
j = n2 - 1
edges = []
Park = set([])
while 1:
i = prev[j]
for k in range(i + 1, j):
edges.append((route2[k], route2[k + 1]))
edges.append((n2 + pstar[i, j], route2[i + 1]))
edges.append((n2 + pstar[i, j], route2[j]))
Park.add(pstar[i, j])
if i == 0:
break
j = i
print("Parking Points=", Park)
使用する中継地点に対する巡回セールスマン問題を解き，第１階層の順回路を求め，結果を描画する．
c1 = {}
pos1 = {}
n1 = len(Park)
for i, p in enumerate(Park):
c1[i, i] = 0
pos1[i] = x1[p], y1[p]
for j, q in enumerate(Park):
c1[i, j] = distance(x1[p], y1[p], x1[q], y1[q])
if n1 <= 10:
V = list(range(n1))
total1, route1 = tspdp(n1, c1, V)
else:
total1, route1, G1 = tsplk(n1, c1)
print(total1, route1)
G = nx.Graph()
tour = route1
for idx, i in enumerate(tour[:-1]):
G.add_edge(i, tour[idx + 1])
G.add_edge(tour[-1], tour[0])
nx.draw(
G,
pos=pos1,
node_size=1000 / n + 10,
with_labels=False,
width=5,
edge_color="orange"
)
G3 = nx.Graph()
G3.add_nodes_from(G2.nodes)
for i in x1: #中継地点の追加
G3.add_node(n2 + i)
pos2[n2 + i] = x1[i], y1[i]
G3.add_edges_from(edges)
nx.draw(G3, pos=pos2, with_labels=False, node_size=10, node_color="y")
plt.show()
