●巡回セールスマン問題
都市が下図のように散らばっているとします。セールスマンは各都市を 1 回ずつもれなく訪問して帰ってこなくてはなりません。このとき、一番短い巡路 (全ての都市を含む単純な閉路) を見つけるのが「巡回セールスマン問題 (Traveling Salesperson Problem, TSP)」です。
巡路の長さは自由に定義してもかまいません。たとえば、時間であったり料金であってもいいのですが、今回はオーソドックスに「距離」と定義することにしましょう。
Ａ Ｂ
● ●
● ●
Ｅ Ｆ
● ●
Ｃ Ｄ
都市の配置
図 : 巡回セールスマン問題
上図の場合、都市の個数は 6 つあるので、巡路の総数は (6 - 1)! / 2 = 60 通りしかありません。これならば総当りで簡単に答えを求めることができます。ところが、都市の個数が多くなると簡単ではありません。n 個の都市の場合は (n - 1)! / 2 の巡路を調べなくてはいけません。巡路の総数が n の階乗に比例して増えるというのは、まさに爆発的に増えるので、n が増えると厳密解を求めるのは大変困難になります。
今回は PuLP を使って、巡回セールスマン問題を解いてみましょう。
●プログラムの作成
今回は拙作のページ「Algorithms with Python: 巡回セールスマン問題 [4]」で作成したプログラムを少し変更して、ファイルからデータを読み込むことにします。あとのプログラムは拙作のページ「ハミルトン閉路」で作成したものとほぼ同じですが、目的関数を設定して最小値を求めるように修正します。プログラムは次のようになります。
リスト : 巡回セールスマン問題
using JuMP, HiGHS # Cbc
using Plots
# ファイルからデータを読み込む
function read_data(filename)
buff = []
open(filename, "r") do fin
for s = eachline(fin)
push!(buff, [parse(Int, x) for x = split(s)])
end
end
buff
end
# 距離を求める
function distance(p1, p2)
dx = p1[1] - p2[1]
dy = p1[2] - p2[2]
sqrt(dx * dx + dy * dy)
end
# 隣接行列の生成
function make_matrix(ps)
size = length(ps)
xs = zeros(size, size)
for i = 1:size
for j = 1:size
xs[i, j] = distance(ps[i], ps[j])
end
end
xs
end
# 解を求める
function solver(ps)
size = length(ps)
ws = make_matrix(ps) # 隣接行列
model = Model(HiGHS.Optimizer)
# model = Model(Cbc.Optimizer)
# 変数の生成
@variable(model, xs[1:size, 1:size], Bin)
@variable(model, 0 <= us[1:size] <= size-1, Int)
# 目的関数
@objective(model, Min, sum(ws .* xs))
# 制約条件
for i = 1:size
@constraint(model, sum(xs[j, i] for j = 1:size) == 1)
@constraint(model, sum(xs[i, j] for j = 1:size) == 1)
end
# 自分自身への辺は選択しない
for i in 1:size
@constraint(model, xs[i, i] == 0)
end
# 頂点に対応する変数の制約
for i = 1:size
for j = 2:size
if i != j
@constraint(model, us[i] + 1 - (size - 1) * (1 - xs[i, j]) <= us[j])
end
end
end
@constraint(model, us[1] == 0) # スタート
for i = 2:size
@constraint(model, 1 <= us[i])
end
# 実行
# println(model)
set_silent(model)
optimize!(model)
# 結果の表示
println("結果: ", termination_status(model))
println("値 ", objective_value(model))
tour = [round(Int, x) for x = value.(us)]
println(tour)
# 描画
draw(ps, tour)
end
データは関数 read_data() で読み込み、それを関数 solver() に渡します。引数 ps から関数 make_matrix() で隣接行列を生成します。目的関数は拙作のページ「最短経路問題」と同じで、選択した辺の長さの合計値になります。これは sum(ws .* xs) で求めることができます。制約条件ですが、まずは「ハミルトン閉路」と同じ条件で試してみましょう。
テストデータは拙作のページ「Algorithms with Python: 巡回セールスマン問題 [4]」と同じものを使用します。
表 : テストデータ
r19.txt r20.txt r21.txt r22.txt r23.txt
--------------------------------------------
422 295 396 162 323 424 264 95 346 106
140 324 96 38 22 396 137 222 140 414
275 118 221 336 205 19 298 156 363 298
282 64 121 380 61 116 323 136 136 268
403 324 416 85 379 156 359 364 321 31
200 335 98 296 75 320 122 283 367 367
377 417 377 78 356 96 79 235 23 428
376 80 44 400 242 423 321 358 375 428
118 40 84 61 78 421 37 59 226 312
358 58 319 380 24 222 80 382 50 388
258 81 117 278 320 257 82 320 85 27
84 380 22 103 157 245 219 119 388 396
275 402 162 381 345 199 417 245 48 328
79 60 244 417 45 121 182 379 143 382
384 320 317 120 319 185 425 102 107 112
57 242 260 256 58 77 223 282 330 336
224 60 280 202 44 184 320 18 228 300
155 316 223 61 235 375 183 155 364 167
44 305 124 119 98 299 144 283 305 308
182 263 38 216 258 220 212 45
238 305 95 165 152 21
84 423 263 352
176 258
都市の数は 19 から 23 までで、配置は乱数で決めたものです。実行結果は次のようになりました。
julia> @time solver(read_data("r19.txt"))
結果: OPTIMAL
値 1444.0588618791194
[0, 7, 15, 16, 1, 5, 3, 18, 12, 17, 14, 8, 4, 11, 2, 10, 13, 6, 9]
0.384641 seconds (32.09 k allocations: 1.774 MiB)
julia> @time solver(read_data("r20.txt"))
結果: OPTIMAL
値 1672.311517728647
[0, 15, 3, 7, 19, 9, 18, 8, 14, 4, 10, 13, 6, 5, 17, 2, 1, 16, 12, 11]
3.033542 seconds (45.56 k allocations: 2.710 MiB)
julia> @time solver(read_data("r21.txt"))
結果: OPTIMAL
値 1623.610871744113
[0, 3, 13, 11, 15, 4, 14, 1, 2, 8, 18, 6, 16, 10, 17, 12, 9, 20, 5, 7, 19]
1.588736 seconds (41.19 k allocations: 2.288 MiB)
julia> @time solver(read_data("r22.txt"))
結果: OPTIMAL
値 1873.8983412008597
[0, 16, 1, 2, 6, 14, 17, 7, 19, 12, 13, 21, 5, 10, 4, 9, 3, 20, 15, 8, 18, 11]
0.665769 seconds (42.51 k allocations: 2.295 MiB)
julia> @time solver(read_data("r23.txt"))
結果: OPTIMAL
値 1694.4037091208925
[0, 11, 21, 7, 1, 20, 10, 18, 13, 9, 4, 19, 8, 12, 5, 17, 14, 22, 16, 2, 3, 15, 6]
2.723834 seconds (51.23 k allocations: 2.943 MiB)
経路図は Julia のライブラリ Plots.jl を使って作成しました。プログラムは次のようになります。
リスト : 経路の描画
function draw(ps, tour)
idx = sortperm(tour)
zss = ps[idx]
push!(zss, ps[1])
xss = [x[1] for x = zss]
yss = [x[2] for x = zss]
plot(xss, yss, linecolor = :blue, linewidth = 2, label = :none)
scatter!([x[1] for x = ps], [x[2] for x = ps],
color = :red, markersize = 6, label = :none)
end
ご参考までに、ソルバーを CBC に変更した場合の実行時間も示します。
表 : 実行結果
秒
: 距離 : HiGHS : CBC
--------+--------+-------+-------
r19.txt : 1444.1 : 0.38 : 0.88
r20.txt : 1672.3 : 3.03 : 6.47
r21.txt : 1623.6 : 1.59 : 8.70
r22.txt : 1873.9 : 0.67 : 6.07
r23.txt : 1694.4 : 2.72 : 7.99
- 実行環境 : Julia 1.12.6, Ubunts 24.04 LTS (WSL2, Windows 11), intel CORE i5-1235U 1.30GHz
- ソルバー : HiGHS v1.24.0, CBC v1.3.0
MTZ 方式の制約はまだ弱いのですが、HiGHS は CBC よりも高速に解くことができました。CBC を使う場合は、制約の強化が必要なようです。そこで、参考 URL 『Python を用いた最適化ソルバー Gurobi 入門』を参考に、MTZ 方式の制約を強化してみましょう。
●MTZ 方式の強化
MTZ 方式は 辺 \(X_{ij}\) (i -> j) を使った制約式ですが、 ここに辺 \(X_{ji}\) (j -> i) の条件を加えます。次の式を見てください。
\(X_{ji}\) が 0 の場合、式 2 は (a) または (b) になります。これは式 1 と同じ制約ですね。次に \(X_{ji}\) が 1 の場合を考えます。\(X_{ij}\) が 0 であれば、式 2 は \(us[i] \leq us[j] + 1\) (c) となります。j -> i に行く辺を選ぶのですから、\(us[i]\) の値は \(us[j] + 1\) になるので式 (c) を満たします。
最後に、\(X_{ij}\) と \(X_{ji}\) が 1 の場合を考えます (式 (d))。\(us[i]\) が最小値の 1 のとき、式 (d) は \(size - 1 \leq us[j]\) になるので条件を満たしますが、それ以外の値では条件を満たさないことがわかります。つまり、\(X_{ij}\) と \(X_{ji}\) は同時に選ぶことはできない、という制約を表していると考えることができます。
プログラムの修正は簡単です
# 頂点に対応する変数の制約
for i = 2:size
for j = 2:size
if i != j
@constraint(model, us[i] + 1 - size * (1 - xs[i, j]) + (size - 3) * xs[j ,i] <= us[j])
end
end
end
MTZ 方式の制約を設定するとき、for ループの変数 i も 1 を除外することをお忘れなく。
実行結果は次のようになりました。
表 : 実行結果 (秒)
: : MTZ0 : MTZ1
: 距離 : HiGHS : CBC : HiGHS : CBC
--------+--------+-------+-------+-------+--------
r19.txt : 1444.1 : 0.38 : 0.88 : 0.14 : 0.25
r20.txt : 1672.3 : 3.03 : 6.47 : 1.29 : 2.04
r21.txt : 1623.6 : 1.59 : 8.70 : 2.08 : 2.54
r22.txt : 1873.9 : 0.67 : 6.07 : 1.07 : 0.54
r23.txt : 1694.4 : 2.72 : 7.99 : 6.67 : 8.53
CBC の場合、MTZ 方式の制約強化はとても効果的で、r23.txt 以外のデータは実行時間を短縮することができました。HiGHS の場合、データによって効果はまちまちで、大幅に遅くなる場合もあります。制約強化の効果は CBC よりも少ないようです。
●MTZ 方式の強化 (2)
制約式の強化はこれだけではありません。us[i] (i != 1) の範囲は 1 から size - 1 までですが、ここにも辺 \(X_{ij}\) と \(X_{ji}\) を使った制約を追加することができます。次の式を見てください。
式 2a は \(U_{i}\) の下界を表し、式 2b は上界を表しています。まず最初に、式 2a から説明しましょう。次の式を見てください。
辺 \(X_{1i}\) は最初に選択する辺、\(X_{i1}\) は最後に選択する辺を表します。a の場合、\(X_{1i}\) は 0 なので、\(U_{i}\) は 2 以上であることがわかります。b の場合、\(X_{1i}\) は 1 なので、\(U_{i}\) は 1 になります。c の場合、\(X_{i1}\) は 1 なので、\(U_{i}\) の値は size - 1 となります。d の場合、最後に選択する辺は条件を満たしますが、最初に選択する辺は条件を満たしません。つまり、2 つの辺を同時に選ぶことはできないわけです。
次に式 2b を説明します。
a の場合、\(X_{i1}\) は 0 なので、\(U_{i}\) の値は size - 2 以下になります。b の場合、\(X_{1i}\) は 1 なので、\(U_{i}\) の値は 1 になります。c の場合、\(X_{i1}\) が 1 なので、\(U_{i}\) の値は size - 1 となります。d の場合、最初に選択する辺は条件を満たしていますが、最後に選択する辺は条件を満たしません。つまり、2 つの辺を同時に選ぶことはできないわけです。
以上の制約条件を設定するプログラムは次のようになります。
リスト : 制約条件の設定
# 頂点に対応する変数の制約
for i = 2:size
for j = 2:size
if i != j
@constraint(model, us[i] + 1 - size * (1 - xs[i, j]) + (size - 3) * xs[j ,i] <= us[j])
end
end
end
# 上界と下界の制約
for i = 2:size
@constraint(model, 1 + (1 - xs[1, i]) + (size - 3) * xs[i, 1] <= us[i])
@constraint(model, us[i] <= (size - 1) - (1 - xs[i, 1]) - (size - 3) * xs[1, i])
end
制約条件の式をそのままプログラムしただけなので、とくに難しいところはないと思います。
それでは実行してみましょう。
表 : 実行結果 (秒)
: : MTZ0 : MTZ1 : MTZ2
: 距離 : HiGHS : CBC : HiGHS : CBC : HiGHS : CBC
--------+--------+-------+-------+-------+-------+-------+-------
r19.txt : 1444.1 : 0.38 : 0.88 : 0.14 : 0.25 : 0.18 : 0.11
r20.txt : 1672.3 : 3.03 : 6.47 : 1.29 : 2.04 : 2.01 : 0.72
r21.txt : 1623.6 : 1.59 : 8.70 : 2.08 : 2.54 : 0.07 : 0.19
r22.txt : 1873.9 : 0.67 : 6.07 : 1.07 : 0.54 : 0.43 : 0.54
r23.txt : 1694.4 : 2.72 : 7.99 : 6.67 : 8.53 : 2.15 : 2.12
上界と下界の制約を強化した効果はとても大きいですね。最適化ソルバーを使う場合、適切な制約条件を設定することがいかに重要であるか、実感することができました。
ご参考までに、都市を 30 から 35 に増やして実行してみました。都市のデータはプログラムリストの中に記述してあります。興味のある方はいろいろ試してみてください。
表 : 実行結果 (秒)
都市 : 距離 : HiGHS : CBC
------+--------+-------+-------
30 : 2107.1 : 8.60 : 32.28
31 : 2110.0 : 2.78 : 1.36
32 : 2000.5 : 0.35 : 0.45
33 : 2407.4 : 21.20 : 54.09
34 : 2318.3 : 0.82 : 1.09
35 : 2332.2 : 11.12 : 3.04
都市数 30
都市数 31
都市数 32
都市数 33
都市数 34
都市数 35
●プログラムリスト
#
# tsp.jl : JuMP による巡回セールスマン問題の解法
#
# Copyright (C) 2026 Makoto Hiroi
#
using JuMP
using HiGHS
# using Cbc
using Plots
# ファイルよりデータを読み込む
function read_data(filename)
buff = []
open(filename, "r") do fin
for s = eachline(fin)
push!(buff, [parse(Int, x) for x = split(s)])
end
end
buff
end
# 距離を求める
function distance(p1, p2)
dx = p1[1] - p2[1]
dy = p1[2] - p2[2]
sqrt(dx * dx + dy * dy)
end
# 隣接行列の生成
function make_matrix(ps)
size = length(ps)
xs = zeros(size, size)
for i = 1:size
for j = 1:size
xs[i, j] = distance(ps[i], ps[j])
end
end
xs
end
# 経路の描画
function draw(ps, tour)
idx = sortperm(tour)
zss = ps[idx]
push!(zss, ps[1])
xss = [x[1] for x = zss]
yss = [x[2] for x = zss]
plot(xss, yss, linecolor = :blue, linewidth = 2, label = :none)
scatter!([x[1] for x = ps], [x[2] for x = ps],
color = :red, markersize = 6, label = :none)
end
# 解を求める
function solver(ps)
size = length(ps)
ws = make_matrix(ps) # 隣接行列
model = Model(HiGHS.Optimizer)
# model = Model(Cbc.Optimizer)
# 変数の生成
@variable(model, xs[1:size, 1:size], Bin)
@variable(model, 0 <= us[1:size] <= size-1, Int)
# 目的関数
@objective(model, Min, sum(ws .* xs))
# 制約条件
for i = 1:size
@constraint(model, sum(xs[j, i] for j = 1:size) == 1)
@constraint(model, sum(xs[i, j] for j = 1:size) == 1)
end
# 自分自身への辺は選択しない
for i in 1:size
@constraint(model, xs[i, i] == 0)
end
# 頂点に対応する変数の制約
for i = 2:size
for j = 2:size
if i != j
@constraint(model, us[i] + 1 - size * (1 - xs[i, j]) + (size - 3) * xs[j ,i] <= us[j])
end
end
end
# 上界と下界の制約
for i = 2:size
@constraint(model, 1 + (1 - xs[1, i]) + (size - 3) * xs[i, 1] <= us[i])
@constraint(model, us[i] <= (size - 1) - (1 - xs[i, 1]) - (size - 3) * xs[1, i])
end
@constraint(model, us[1] == 0) # スタート
for i = 2:size
@constraint(model, 1 <= us[i])
end
# 実行
# println(model)
set_silent(model)
optimize!(model)
# 結果の表示
println("結果: ", termination_status(model))
println("値 ", objective_value(model))
tour = [round(Int, x) for x = value.(us)]
println(tour)
# 描画
draw(ps, tour)
end
# データ
ps30 = [[255, 352], [431, 109], [354, 238], [141, 32], [67, 351], [417, 112], [395, 108],
[398, 237], [162, 230], [488, 255], [464, 304], [436, 28], [47, 206], [181, 480],
[424, 343], [106, 287], [429, 134], [408, 31], [63, 75], [429, 225], [158, 345],
[307, 13], [129, 339], [120, 479], [143, 133], [92, 194], [225, 349], [130, 446],
[454, 156], [210, 70]]
ps31 = [[183, 411], [279, 134], [373, 16], [221, 171], [144, 21], [341, 48], [401, 260],
[198, 159], [121, 245], [234, 411], [249, 443], [123, 53], [69, 128], [108, 439],
[34, 457], [205, 315], [106, 237], [244, 302], [349, 238], [54, 101], [256, 450],
[200, 42], [223, 219], [186, 171], [201, 414], [170, 156], [15, 216], [324, 194],
[289, 139], [405, 400], [401, 205]]
ps32 = [[221, 82], [78, 417], [235, 434], [453, 409], [220, 341], [352, 162], [190, 460],
[463, 390], [414, 428], [101, 195], [300, 327], [109, 365], [278, 68],
[143, 141], [244, 367], [173, 67], [406, 181], [49, 313], [447, 86], [146, 377],
[382, 207], [226, 140], [356, 416], [374, 354], [100, 111], [455, 159],
[398, 45], [373, 241], [133, 91], [451, 322], [50, 298], [157, 159]]
ps33 = [[17, 93], [83, 282], [52, 429], [238, 258], [81, 199], [309, 169], [207, 48],
[179, 132], [169, 65], [390, 366], [108, 217], [37, 160], [266, 324], [16, 219],
[68, 222], [120, 430], [11, 36], [260, 454], [489, 51], [44, 378], [370, 229],
[362, 272], [283, 337], [191, 335], [266, 313], [475, 211], [377, 97],
[263, 199], [479, 240], [44, 418], [75, 395], [122, 300], [428, 323]]
ps34 = [[271, 18], [448, 428], [176, 326], [432, 207], [405, 367], [116, 294],
[459, 466], [103, 294], [336, 20], [417, 458], [130, 351], [307, 165],
[86, 39], [118, 161], [213, 199], [169, 43], [361, 228], [211, 485], [422, 446],
[69, 362], [169, 238], [356, 310], [344, 107], [99, 37], [22, 38], [325, 484],
[486, 248], [41, 370], [478, 442], [210, 480], [473, 126], [66, 282],
[129, 265], [130, 26]]
ps35 = [[341, 207], [196, 401], [15, 399], [367, 51], [30, 136], [425, 172], [278, 188],
[307, 164], [263, 54], [220, 122], [294, 165], [257, 11], [243, 121],
[144, 390], [156, 263], [49, 230], [293, 444], [477, 299], [377, 424],
[337, 239], [197, 47], [204, 109], [435, 100], [468, 308], [55, 446], [386, 325],
[483, 282], [38, 315], [147, 28], [398, 131], [100, 457], [437, 108], [265, 317],
[366, 37], [205, 178]]
●フロー型部分巡回路除去制約
今回は flow formulation (単品種フロー) と呼ばれる部分巡回路除去制約を紹介します。セールスマンは「もの」を所持していて、都市を訪問するたびに「もの」を一つ置いていきます。n 個の都市を巡回するのであれば、セールスマンは n 個の「もの」を持っていて、都市を移動するたびに持っている「もの」は一つずつ減っていきます。セールスマンが運ぶ「もの」の個数をフローとみなして定式化するのが単品種フローの考え方です。
簡単な例を示しましょう。次の図を見てください。
: A B C D
A -- 3 -- B ---+-------------
| | A : X 3 0 0
0 2 B : 0 X 2 0
| | C : 0 0 X 1
D -- 1 -- C D : 0 0 0 X
都市 A から出発して B, C, D の順番で巡回します。都市 i から j にを移動するときサラリーマンが持っている「もの」の個数を \(f_{ij}\) とすると、\(f_{AB} = 3, f_{BC} = 2, f_{CD} = 1, f_{DA} = 0\) になります。これを表にすると上図 (右) のようになり、次式に示すような関係が成り立ちます。
出発点 A を除いて、頂点から流出するフローから流入するフローを引き算すると -1 になります。出発点は特別で、都市の個数を n とすると、フローの差は n - 1 になります。
フローを表す変数 \(f_{ij}\) を使って制約条件を記述すると次のようになります。
このほかにも、変数 \(f_{ij}\) には制約条件があります。\(f_{ij}\) の値は、辺 \(X_{ij}\) が 1 であれば n - 1 以下になり、\(X_{ij}\) が 0 であれば 0 になります。したがって、制約条件は次式のようになります。
あとはこれを JuMP でプログラムするだけです。次のリストを見てください。
リスト : 単品種フロー型部分巡回路除去制約
function solver(ps)
size = length(ps)
ws = make_matrix(ps) # 隣接行列
model = Model(HiGHS.Optimizer)
# 変数の生成
@variable(model, xs[1:size, 1:size], Bin)
@variable(model, 0 <= fs[1:size, 1:size], Int)
# 目的関数
@objective(model, Min, sum(ws .* xs))
# 制約条件
for i = 1:size
@constraint(model, sum(xs[j, i] for j = 1:size) == 1)
@constraint(model, sum(xs[i, j] for j = 1:size) == 1)
end
# 自分自身への辺は選択しない
for i = 1:size
@constraint(model, xs[i, i] == 0)
end
# 単品種フロー型部分巡回路除去制約
for i = 1:size
for j = 1:size
if i != j
@constraint(model, fs[i, j] <= (size - 1) * xs[i, j])
end
end
end
@constraint(model, sum(fs[1,i] for i = 2:size) - sum(fs[i,1] for i = 2:size) == size - 1)
for i = 2:size
@constraint(model, sum(fs[i, j] for j = 1:size if i != j) -
sum(fs[j, i] for j = 1:size if i != j) == -1)
end
# 実行
# println(model)
set_silent(model)
optimize!(model)
# 結果の表示
println("結果: ", termination_status(model))
println("値 ", objective_value(model))
end
フロー \(f_{ij}\) を表す変数は配列 fs に格納します。あとは制約条件をそのままプログラムしただけなので、とくに難しいところはないと思います。
実行結果を示します。
表 : 実行結果 (HiGHS, 単位 秒)
: 距離 : MTZ0 : MTZ1 : MTZ2 : FLOW
--------+--------+-------+--------+--------+-------
r19.txt : 1444.1 : 0.38 : 0.14 : 0.18 : 0.37
r20.txt : 1672.3 : 3.03 : 1.29 : 2.01 : 0.69
r21.txt : 1623.6 : 1.59 : 2.08 : 0.07 : 0.71
r22.txt : 1873.9 : 0.67 : 1.07 : 0.43 : 5.56
r23.txt : 1694.4 : 2.72 : 6.67 : 2.15 : 2.78
単品種フロー (FLOW) はデータによって実行速度が大きく変わりますが、おおむね MTZ0 よりも速くて MTZ2 よりも遅くなる場合が多いようです。そこで、単品種フローの制約をもう少しだけ強化してみましょう。
●制約の強化
制約の強化は簡単です。変数 \(f_{ij}\) の値が n - 1 になるのは i が 1 のときだけです。それ以外の変数の値は n - 2 以下になります。これを制約として追加しましょう。プログラムの修正も簡単です。次のリストを見てください。
リスト : 制約の強化
# 単品種フロー型部分巡回路除去制約 (強化版)
for i = 2:size
@constraint(model, fs[1, i] <= (size - 1) * xs[1, i])
end
for i = 2:size
for j = 1:size
if i != j
@constraint(model, fs[i, j] <= (size - 2) * xs[i, j])
end
end
end
fs[1, i] の上限値を (size - 1) * xs[1, i] とし、それ以外の変数 fs[i, j] の上限値は (size - 2) * xs[i, j] とするだけです。
実行結果は次のようになりました。
表 : 実行結果 (HiGHS, 単位 秒)
: 距離 : MTZ0 : MTZ1 : MTZ2 : FLOW : FLOW1
--------+--------+-------+--------+--------+--------+--------
r19.txt : 1444.1 : 0.38 : 0.14 : 0.18 : 0.37 : 0.92
r20.txt : 1672.3 : 3.03 : 1.29 : 2.01 : 0.69 : 1.09
r21.txt : 1623.6 : 1.59 : 2.08 : 0.07 : 0.71 : 0.16
r22.txt : 1873.9 : 0.67 : 1.07 : 0.43 : 5.56 : 3.88
r23.txt : 1694.4 : 2.72 : 6.67 : 2.15 : 2.78 : 4.35
データによって効果はまちまちです。単品種フローの場合、制約を強化したからといって、高速になるとは限らないようです。制約式によって得手不得手があるのかもしれませんね。ご参考までに都市の数を 30 から 35 に増やしたときの実行結果を示します。
表 : 実行結果 (HiGHS, 秒)
都市 : 距離 : MTZ2 : FLOW : FLOW1
------+--------+-------+-------+-------
30 : 2107.1 : 8.60 : 5.22 : 3.59
31 : 2110.0 : 2.78 : 5.71 : 43.88
32 : 2000.5 : 0.35 : 2.34 : 3.06
33 : 2407.4 : 21.20 : 21.34 : 50.69
34 : 2318.3 : 0.82 : 6.56 : 6.60
35 : 2332.2 : 11.12 : 26.42 : 10.12
興味のある方はいろいろ試してみてください。
●参考 URL
- 汎用MIPソルバによる巡回セールスマン問題の求解 : 多項式オーダ本数の部分巡回路除去制約 [PDF], (沼田一道さん)
- ロジスティクス工学 07 配送計画, (宮本裕一郎さん)